From eb1d76b0841fa77d7158ec8392d5dcfbfba8949e Mon Sep 17 00:00:00 2001 From: rhaas Date: Fri, 4 Nov 2011 18:34:36 +0000 Subject: second iteration of constraint transport * fix some indices * move poison loop to proper spot * compute conserved divergence in track_divB * add Maple worksheet to check constraint transport indices git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@298 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- doc/constraint_transport.mw | 268 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 268 insertions(+) create mode 100644 doc/constraint_transport.mw (limited to 'doc') diff --git a/doc/constraint_transport.mw b/doc/constraint_transport.mw new file mode 100644 index 0000000..020e75c --- /dev/null +++ b/doc/constraint_transport.mw @@ -0,0 +1,268 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +Constraint transport scheme in GRHydro +The current constraint transport scheme in GRHydro is based on WhiskyMHD (arXiv:gr-qc/0701109) and Bruno Giacomazzo's PhD thesis (http://www.sissa.it/ap/PhD/Theses/giacomazzo.ps.gz). It assumes that all magnetic field components live at the cell center and uses the Reconstruction step to move them to the cell interfaces. +This worksheet helps checking the expressions entered into the code by explicitely checking that div(\134dot B) = 0. + + + + +restart; + + + + +# we later expand taylor series in eps which helps counting orders in the stepsize +dx:=eps*dX;dy:=eps*dY;dz:=eps*dZ; + + +
+<Text-field style="Heading 1" layout="Heading 1">Fluxes of the magnetic field components</Text-field> + + +BconsxfluxX := (i,j,k) -> 0; + + + + +BconsyfluxX := (i,j,k) -> add(add(add(BconsyfluxX_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); + + + + +BconszfluxX := (i,j,k) -> add(add(add(BconszfluxX_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); + + + + +# enforce F^I(B^j) = -F^j(B^i) if desired to check that RHS computed agrees to some order with the one used in divergence cleaning +BconsxfluxY := (i,j,k) -> add(add(add(BconsxfluxY_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); #-BconsyfluxX(i,j,k); + + + + +BconsyfluxY := (i,j,k) -> 0; + + + + +BconszfluxY := (i,j,k) -> add(add(add(BconszfluxY_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); + + + + +BconsxfluxZ := (i,j,k) -> add(add(add(BconsxfluxZ_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); #-BconszfluxX(i,j,k); + + + + +BconsyfluxZ := (i,j,k) -> add(add(add(BconsyfluxZ_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); #-BconszfluxY(i,j,k); + + + + +BconszfluxZ := (i,j,k) -> 0; + + +
+
+<Text-field style="Heading 1" layout="Heading 1">Lapse</Text-field> + + +alp := (i,j,k) -> add(add(add(alp_coeff[_i,_j,_k] * (i*dx)^_i * (j*dy)^_j * (k*dz)^_k,_i=0..3),_j=0..3),_k=0..3); + + +
+
+<Text-field style="Heading 1" layout="Heading 1"><Font encoding="UTF-8">Electric fields (eventually from E=v \134cross B=F(B))</Font></Text-field> + + +Ey := (i,j,k) -> 1/4*(alp(i+1/2,j,k)*BconszfluxX(i+1/2,j,k) + alp(i+1/2,j,k+1)*BconszfluxX(i+1/2,j,k+1)) - + 1/4*(alp(i,j,k+1/2)*BconsxfluxZ(i,j,k+1/2) + alp(i+1,j,k+1/2)*BconsxfluxZ(i+1,j,k+1/2)); + + + + +Ez := (i,j,k) -> 1/4*(alp(i,j+1/2,k)*BconsxfluxY(i,j+1/2,k) + alp(i+1,j+1/2,k)*BconsxfluxY(i+1,j+1/2,k)) - + 1/4*(alp(i+1/2,j,k)*BconsyfluxX(i+1/2,j,k) + alp(i+1/2,j+1,k)*BconsyfluxX(i+1/2,j+1,k)); + + + + +Ex := (i,j,k) -> 1/4*(alp(i,j,k+1/2)*BconsyfluxZ(i,j,k+1/2) + alp(i,j+1,k+1/2)*BconsyfluxZ(i,j+1,k+1/2)) - + 1/4*(alp(i,j+1/2,k)*BconszfluxY(i,j+1/2,k) + alp(i,j+1/2,k+1)*BconszfluxY(i,j+1/2,k+1)); + + +
+
+<Text-field style="Heading 1" layout="Heading 1">RHS for B</Text-field> + + +Bconsrhsx := (i,j,k) -> -1/2*(Ey(i-1,j ,k-1)-Ey(i-1,j ,k ))/dz + -1/2*(Ez(i-1,j ,k )-Ez(i-1,j-1,k ))/dy + -1/2*(Ey(i ,j ,k-1)-Ey(i ,j ,k ))/dz + -1/2*(Ez(i ,j ,k )-Ez(i ,j-1,k ))/dy: + + + + +Bconsrhsy := (i,j,k) -> -1/2 * (Ez(i-1,j-1,k )-Ez(i ,j-1,k )) / dx + -1/2 * (Ex(i ,j-1,k )-Ex(i ,j-1,k-1)) / dz + -1/2 * (Ez(i-1,j ,k )-Ez(i ,j ,k )) / dx + -1/2 * (Ex(i ,j ,k )-Ex(i ,j ,k-1)) / dz: + + + + +Bconsrhsz := (i,j,k) -> - 1/2 * (Ex(i ,j-1,k-1)-Ex(i ,j ,k-1)) / dy + - 1/2 * (Ey(i ,j ,k-1)-Ey(i-1,j ,k-1)) / dx + - 1/2 * (Ex(i ,j-1,k )-Ex(i ,j ,k )) / dy + - 1/2 * (Ey(i ,j ,k )-Ey(i-1,j ,k )) / dx: + + +
+
+<Text-field style="Heading 1" layout="Heading 1">Check RHS at center of cell</Text-field> + + +Bconsrhsx(0,0,0): + + + + +simplify(%): + + + + +algsubs(eps=0,%); + + + + + + + + + +(alp(0,-1/2,0)*BconsxfluxY(0,-1/2,0)-alp(0,1/2,0)*BconsxfluxY(0,1/2,0))/dy + +(alp(0,0,-1/2)*BconsxfluxZ(0,0,-1/2)-alp(0,0,1/2)*BconsxfluxZ(0,0,1/2))/dz: + + + + +collect(%,eps):algsubs(eps^2=0,%); + + +
+
+<Text-field style="Heading 1" layout="Heading 1">Check constraint transport property for edge based divergence operator (see Bruno's thesis)</Text-field> + + +divB := (i,j,k) -> (Bconsrhsx(i+1,j,k)-Bconsrhsx(i,j,k)+Bconsrhsx(i+1,j+1,k)-Bconsrhsx(i,j+1,k)+Bconsrhsx(i+1,j,k+1)-Bconsrhsx(i,j,k+1)+Bconsrhsx(i+1,j+1,k+1)-Bconsrhsx(i,j+1,k+1))/4/dx + +(Bconsrhsy(i,j+1,k)-Bconsrhsy(i,j,k)+Bconsrhsy(i+1,j+1,k)-Bconsrhsy(i+1,j,k)+Bconsrhsy(i,j+1,k+1)-Bconsrhsy(i,j,k+1)+Bconsrhsy(i+1,j+1,k+1)-Bconsrhsy(i+1,j,k+1))/4/dy + +(Bconsrhsz(i,j,k+1)-Bconsrhsz(i,j,k)+Bconsrhsz(i+1,j,k+1)-Bconsrhsz(i+1,j,k)+Bconsrhsz(i,j+1,k+1)-Bconsrhsz(i,j+1,k)+Bconsrhsz(i+1,j+1,k+1)-Bconsrhsz(i+1,j+1,k))/4/dz; + + + + +simplify(divB(0,0,0)); + + +
+ + + + + + + + \ No newline at end of file -- cgit v1.2.3