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"><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));