aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_AdvectedLoopM.F90
diff options
context:
space:
mode:
Diffstat (limited to 'src/GRHydro_AdvectedLoopM.F90')
-rw-r--r--src/GRHydro_AdvectedLoopM.F90262
1 files changed, 262 insertions, 0 deletions
diff --git a/src/GRHydro_AdvectedLoopM.F90 b/src/GRHydro_AdvectedLoopM.F90
new file mode 100644
index 0000000..81a3ede
--- /dev/null
+++ b/src/GRHydro_AdvectedLoopM.F90
@@ -0,0 +1,262 @@
+ /*@@
+ @file GRHydro_AdvectedLoopM.F90
+ @date Aug 15, 2011
+ @author Scott Noble, Joshua Faber, Bruno Mundim
+ @desc
+ Advected loop test as implemented by Beckwith and Stone Astrophys.J.Suppl.
+ 193 (2011) 6, arXiv:1101.3573.
+
+
+ Other relevant references: Devore JCP 92, 142 (1991),
+ Toth and Odstrcil JCP 128,82 (1996),
+ Gardiner and Stone JCP 227, 4123 (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_AdvectedLoopM
+ @date Aug 11, 2011
+ @author Scott Noble, Joshua Faber, Bruno Mundim
+ @desc
+ Initial data for advected loop test
+ @enddesc
+ @calls
+ @calledby
+ @history
+ Using GRHydro_ShockTubeM.F90 as a template.
+ @endhistory
+
+@@*/
+
+subroutine GRHydro_AdvectedLoopM(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_FUNCTIONS
+
+ CCTK_INT :: i, j, k, nx, ny, nz
+ CCTK_REAL :: det,radius,vxval,vyval,vzval,gam
+ CCTK_REAL :: radius_iph, radius_imh
+ CCTK_REAL :: radius_jph, radius_jmh
+ CCTK_REAL :: radius_kph, radius_kmh
+ CCTK_REAL :: rhoval,pressval
+ CCTK_REAL :: r_loop,A_loop,pi
+ CCTK_REAL :: dx,dy,dz
+ CCTK_REAL :: range_x,range_y,range_z,range_d
+ CCTK_REAL :: cos_theta, sin_theta, tan_theta
+ CCTK_REAL :: Bvecx_d, Bvecy_d, Bvecz_d
+ CCTK_REAL :: x_d, y_d, z_d
+ CCTK_REAL :: dx_d, dy_d, dz_d, dx_x, dz_x
+
+!!$Adiabatic index for test:
+ gam = (5.0d0/3.0d0)
+
+!!$radius of the loop:
+ r_loop = 0.3d0
+
+!!$ stregth of the A-field
+ A_loop=1.0d-3
+
+!!$pressure and density:
+ rhoval = 1.0d0
+ pressval = 3.0d0
+
+!!$ Vx, Vy and Vz values:
+ if (CCTK_EQUALS(advectedloop_case,"V^z/=0")) then
+ vxval=0.2d0/sqrt(6.0d0)
+ vyval=0.5d0*vxval
+ vzval=vyval
+ else if (CCTK_EQUALS(advectedloop_case,"V^z=0")) then
+ vxval=0.2d0/sqrt(6.0d0)
+ vyval=0.5d0*vxval
+ vzval=0.0d0
+ else
+ call CCTK_WARN(0,"V^z component case not recognized!")
+ end if
+
+ 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)
+
+!!$ Note that the 3D test wasn't deviced to be used with AMR!
+ range_x = (cctk_gsh(1)-1)*dx
+ range_y = (cctk_gsh(2)-1)*dy
+ range_z = (cctk_gsh(3)-1)*dz
+
+ range_d = sqrt(range_z**2+range_x**2)
+
+ cos_theta = range_z/range_d
+ sin_theta = range_x/range_d
+ tan_theta = sin_theta/cos_theta
+
+ do i=1,nx
+ do j=1,ny
+ do k=1,nz
+
+ rho(i,j,k)=rhoval
+ press(i,j,k)=pressval
+ eps(i,j,k)=press(i,j,k)/(gam-1.0d0)/rho(i,j,k)
+
+ if (CCTK_EQUALS(advectedloop_type,"2D")) then
+ velx(i,j,k)=vxval
+ vely(i,j,k)=vyval
+ velz(i,j,k)=vzval
+ Bvecz(i,j,k)=0.0d0
+
+ radius = sqrt(x(i,j,k)**2+y(i,j,k)**2)
+
+ if (CCTK_EQUALS(advectedloop_delA,"Exact")) then
+
+ if(radius.le.r_loop) then
+ Bvecx(i,j,k)=-1.0d0*A_loop*y(i,j,k)/radius
+ Bvecy(i,j,k)=A_loop*x(i,j,k)/radius
+ else
+ Bvecx(i,j,k)=0.0d0
+ Bvecy(i,j,k)=0.0d0
+ endif
+
+ else if (CCTK_EQUALS(advectedloop_delA,"Numeric")) then
+
+ radius_iph = max(sqrt((x(i,j,k)+0.5d0*dx)**2+y(i,j,k)**2)-r_loop,0.d0)
+ radius_imh = max(sqrt((x(i,j,k)-0.5d0*dx)**2+y(i,j,k)**2)-r_loop,0.d0)
+ radius_jph = max(sqrt(x(i,j,k)**2+(y(i,j,k)+0.5d0*dy)**2)-r_loop,0.d0)
+ radius_jmh = max(sqrt(x(i,j,k)**2+(y(i,j,k)-0.5d0*dy)**2)-r_loop,0.d0)
+
+!! if(radius.le.r_loop) then
+ Bvecx(i,j,k)=-1.0d0*A_loop*(radius_jph-radius_jmh)/dy
+ Bvecy(i,j,k)=A_loop*(radius_iph-radius_imh)/dx
+!! else
+!! Bvecx(i,j,k)=0.0d0
+!! Bvecy(i,j,k)=0.0d0
+!! endif
+
+ else
+ call CCTK_WARN(0,"A^b differentiation not recognized!")
+ end if
+
+ else if (CCTK_EQUALS(advectedloop_type,"3D")) then
+
+ velx(i,j,k)=cos_theta*vxval-sin_theta*vzval
+ vely(i,j,k)=vyval
+ velz(i,j,k)=cos_theta*vzval+sin_theta*vxval
+
+ Bvecz_d=0.0d0
+
+ x_d = cos_theta*x(i,j,k)+sin_theta*z(i,j,k)
+ y_d = y(i,j,k)
+ z_d = cos_theta*z(i,j,k)-sin_theta*x(i,j,k)
+
+ radius = sqrt(x_d**2+y_d**2)
+
+ if (CCTK_EQUALS(advectedloop_delA,"Exact")) then
+
+ if(radius.le.r_loop) then
+ Bvecx_d=-1.0d0*A_loop*y_d/radius
+ Bvecy_d=A_loop*x_d/radius
+ else
+ Bvecx_d=0.0d0
+ Bvecy_d=0.0d0
+ endif
+
+ Bvecx(i,j,k)=cos_theta*Bvecx_d-sin_theta*Bvecz_d
+ Bvecy(i,j,k)=Bvecy_d
+ Bvecz(i,j,k)=cos_theta*Bvecz_d+sin_theta*Bvecx_d
+
+ else if (CCTK_EQUALS(advectedloop_delA,"Numeric")) then
+
+!! dx_d = cos_theta*dx+sin_theta*dz
+!! dy_d = dy
+!! dz_d = cos_theta*dz-sin_theta*dx
+
+!! dx_d is the change in the rotated coords induced by a step in a direction over the Cartesian grid
+ dx_x = cos_theta*dx
+ dz_x = sin_theta*dz
+
+!! These are used for exact differencing
+ radius_iph = max(sqrt((x_d+0.5d0*dx_x)**2+y_d**2)-r_loop,0.d0)
+ radius_imh = max(sqrt((x_d-0.5d0*dx_x)**2+y_d**2)-r_loop,0.d0)
+ radius_jph = max(sqrt(x_d**2+(y_d+0.5d0*dy)**2)-r_loop,0.d0)
+ radius_jmh = max(sqrt(x_d**2+(y_d-0.5d0*dy)**2)-r_loop,0.d0)
+ radius_kph = max(sqrt((x_d+0.5d0*dz_x)**2+y_d**2)-r_loop,0.d0)
+ radius_kmh = max(sqrt((x_d-0.5d0*dz_x)**2+y_d**2)-r_loop,0.d0)
+
+!! see notes
+!! if(radius.le.r_loop) then
+ Bvecx(i,j,k)=-1.0d0*A_loop*cos_theta*(radius_jph-radius_jmh)/dy
+ Bvecy(i,j,k)=A_loop*(sin_theta*(radius_kph-radius_kmh)/dz + &
+ cos_theta*(radius_iph-radius_imh)/dx)
+ Bvecz(i,j,k)=-1.0d0*A_loop*sin_theta*(radius_jph-radius_jmh)/dy
+!! else
+!! Bvecx(i,j,k)=0.0d0
+!! Bvecy(i,j,k)=0.0d0
+!! Bvecz(i,j,k)=0.0d0
+!! endif
+
+ else
+ call CCTK_WARN(0,"A^b differentiation not recognized!")
+ end if
+
+ else
+ call CCTK_WARN(0,"Advected loop type not recognized!")
+ end if
+
+ 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_AdvectedLoopM
+
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-03 12:10:19 +0000