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;
Fluxes of the magnetic field components
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;
Electric fields (eventually from E=v \134cross B=F(B))
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));
RHS for B
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:
Check RHS at center of cell
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,%);
Check constraint transport property for edge based divergence operator (see Bruno's thesis)
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));