diff options
Diffstat (limited to 'src/GRHydro_PoloidalMagFieldM.F90')
| -rw-r--r-- | src/GRHydro_PoloidalMagFieldM.F90 | 152 |
1 files changed, 152 insertions, 0 deletions
diff --git a/src/GRHydro_PoloidalMagFieldM.F90 b/src/GRHydro_PoloidalMagFieldM.F90 new file mode 100644 index 0000000..38763ed --- /dev/null +++ b/src/GRHydro_PoloidalMagFieldM.F90 @@ -0,0 +1,152 @@ + /*@@ + @file GRHydro_PoloidalMagFieldM.F90 + @date Oct 31, 2011 + @author Bruno Mundim, Joshua Faber, Scott Noble + @desc + Poloidal Magnetic field implemented as in "General relativistic +simulations of magnetized binary neutron star mergers" - +by Yuk Tung Liu, Stuart L. Shapiro, Zachariah B. Etienne, and +Keisuke Taniguchi - Phys. Rev. D 78, 024012 (2008) + + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#include "GRHydro_Macros.h" + +#define velx(i,j,k) vel(i,j,k,1) +#define vely(i,j,k) vel(i,j,k,2) +#define velz(i,j,k) vel(i,j,k,3) +#define sx(i,j,k) scon(i,j,k,1) +#define sy(i,j,k) scon(i,j,k,2) +#define sz(i,j,k) scon(i,j,k,3) +#define Bvecx(i,j,k) Bvec(i,j,k,1) +#define Bvecy(i,j,k) Bvec(i,j,k,2) +#define Bvecz(i,j,k) Bvec(i,j,k,3) +#define Bconsx(i,j,k) Bcons(i,j,k,1) +#define Bconsy(i,j,k) Bcons(i,j,k,2) +#define Bconsz(i,j,k) Bcons(i,j,k,3) + + /*@@ + @routine GRHydro_PoloidalMagFieldM + @date Oct 31, 2011 + @author Bruno Mundim, Joshua Faber, Scott Noble + @desc + Poloidal Magnetic field implemented as in "General relativistic +simulations of magnetized binary neutron star mergers" - +by Yuk Tung Liu, Stuart L. Shapiro, Zachariah B. Etienne, and +Keisuke Taniguchi - Phys. Rev. D 78, 024012 (2008) + @enddesc + @calls + @calledby + @history + Using GRHydro_ShockTubeM.F90 as a template. + @endhistory + +@@*/ + +subroutine GRHydro_PoloidalMagFieldM(CCTK_ARGUMENTS) + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, nx, ny, nz + CCTK_REAL :: det + CCTK_REAL :: dx,dy,dz + CCTK_REAL :: rhofac, delPcut, maxP_Pcut + CCTK_REAL :: Aphi, Ax, Ay, Az + CCTK_REAL :: rho_dx, rho_dy, rho_dz + CCTK_REAL :: press_dx, press_dy, press_dz + CCTK_REAL :: Aphi_dx, Aphi_dy, Aphi_dz + CCTK_REAL :: Ax_dy, Ax_dz, Ay_dx, Ay_dz + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = CCTK_DELTA_SPACE(1) + dy = CCTK_DELTA_SPACE(2) + dz = CCTK_DELTA_SPACE(3) + + do i=1,nx + do j=1,ny + do k=1,nz + + rhofac = 1.0d0-rho(i,j,k)/poloidal_rho_max + delPcut = press(i,j,k)-poloidal_P_cut + maxP_Pcut = max(delPcut,0.0d0) + Aphi = poloidal_A_b*rhofac**poloidal_n_p*maxP_Pcut + Ax = -y(i,j,k)*Aphi + Ay = x(i,j,k)*Aphi + Az = 0.0 + +!! write(*,*)'Before accessing rho(i,k,k)' +!! write(*,*)'rho(',i,',',j,',',k,') = ', rho(i,j,k) +!! write(*,*)'Ax, Ay, Az, Aphi = ', Ax, Ay, Az,Aphi +!! write(*,*)'rhofac = ', rhofac +!! write(*,*)'delPcut = ', delPcut +!! write(*,*)'maxP_Pcut = ', maxP_Pcut + + rho_dx = 0.5d0*(rho(i+1,j,k)-rho(i-1,j,k))/dx + rho_dy = 0.5d0*(rho(i,j+1,k)-rho(i,j-1,k))/dy + rho_dz = 0.5d0*(rho(i,j,k+1)-rho(i,j,k-1))/dz + press_dx = 0.5d0*(press(i+1,j,k)-press(i-1,j,k))/dx + press_dy = 0.5d0*(press(i,j+1,k)-press(i,j-1,k))/dy + press_dz = 0.5d0*(press(i,j,k+1)-press(i,j,k-1))/dz + + Aphi_dx = poloidal_A_b*rhofac**(poloidal_n_p-1)*maxP_Pcut* & + ( rhofac*press_dx/delPcut - poloidal_n_p*rho_dx/poloidal_rho_max ) + Aphi_dy = poloidal_A_b*rhofac**(poloidal_n_p-1)*maxP_Pcut* & + ( rhofac*press_dy/delPcut - poloidal_n_p*rho_dy/poloidal_rho_max ) + Aphi_dz = poloidal_A_b*rhofac**(poloidal_n_p-1)*maxP_Pcut* & + ( rhofac*press_dz/delPcut - poloidal_n_p*rho_dz/poloidal_rho_max ) + + Ax_dy = -Aphi - y(i,j,k)*Aphi_dy + Ax_dz = y(i,j,k)*Aphi_dz + Ay_dx = Aphi + x(i,j,k)*Aphi_dx + Ay_dz = x(i,j,k)*Aphi_dz + + Bvecx(i,j,k) = alp(i,j,k)*Ay_dz + Bvecy(i,j,k) = -alp(i,j,k)*Ax_dz + Bvecz(i,j,k) = alp(i,j,k)*(Ax_dy-Ay_dx) + + det=SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + if (CCTK_EQUALS(GRHydro_eos_type,"Polytype")) then + call Prim2ConPolyM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),& + w_lorentz(i,j,k)) + else + call Prim2ConGenM(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& + gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),Bconsx(i,j,k),Bconsy(i,j,k),Bconsz(i,j,k),rho(i,j,k),& + velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),Bvecx(i,j,k),Bvecy(i,j,k),Bvecz(i,j,k),& + w_lorentz(i,j,k)) + end if + + enddo + enddo + enddo + + densrhs = 0.d0 + srhs = 0.d0 + taurhs = 0.d0 + Bconsrhs = 0.d0 + + return + +end subroutine GRHydro_PoloidalMagFieldM + + |