diff options
| author | knarf <knarf@ac85fae7-cede-4708-beff-ae01c7fa1c26> | 2013-12-06 21:42:22 +0000 |
|---|---|---|
| committer | knarf <knarf@ac85fae7-cede-4708-beff-ae01c7fa1c26> | 2013-12-06 21:42:22 +0000 |
| commit | 67de59308d16c563d15aa7a26bdedb379578c087 (patch) | |
| tree | 951a36355f1a2547762a8757f4cddf16b80c45e3 | |
| parent | 8f011ea25c87d3167ef56172186d6f37bc546032 (diff) | |
Add a parameter for an exponent of the pressure term in a poloidal magnetic
field configuration. The default is set to 1, maintaining current status. The
testsuite using this still passes after applying this patch, with identical
files.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/GRHydro_InitData/trunk@220 ac85fae7-cede-4708-beff-ae01c7fa1c26
| -rw-r--r-- | param.ccl | 5 | ||||
| -rw-r--r-- | src/GRHydro_PoloidalMagFieldM.F90 | 23 |
2 files changed, 21 insertions, 7 deletions
@@ -462,6 +462,11 @@ CCTK_REAL poloidal_P_cut "Pressure used to confine the B field inside a star" (0:* :: "Anything positive." } 1.0e-8 +CCTK_INT poloidal_P_p "Index of pressure factor" +{ + (0:* :: "Any non-negative integer" +} 1 + CCTK_REAL poloidal_rho_max "Maximum initial density" { (0:* :: "Anything positive." diff --git a/src/GRHydro_PoloidalMagFieldM.F90 b/src/GRHydro_PoloidalMagFieldM.F90 index c809001..d8c454c 100644 --- a/src/GRHydro_PoloidalMagFieldM.F90 +++ b/src/GRHydro_PoloidalMagFieldM.F90 @@ -88,6 +88,10 @@ subroutine GRHydro_PoloidalMagFieldM(CCTK_ARGUMENTS) ! Initialize to zero Bvec = 0.0d0 + Aphi_dx = 0.0 + Aphi_dy = 0.0 + Aphi_dz = 0.0 + do i=2,nx-1 do j=2,ny-1 do k=2,nz-1 @@ -95,7 +99,7 @@ subroutine GRHydro_PoloidalMagFieldM(CCTK_ARGUMENTS) rhofac = 1.0d0-rho(i,j,k)/poloidal_rho_max delPcut = press(i,j,k)-poloidal_P_cut maxP_Pcut = max(delPcut,0.0d0) - AphiL = poloidal_A_b*rhofac**poloidal_n_p*maxP_Pcut + AphiL = poloidal_A_b*rhofac**poloidal_n_p*maxP_Pcut**poloidal_P_p Ax = -y(i,j,k)*AphiL Ay = x(i,j,k)*AphiL Az = 0.0 @@ -115,12 +119,17 @@ subroutine GRHydro_PoloidalMagFieldM(CCTK_ARGUMENTS) 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 ) + if (maxP_Pcut > 0.0) then + Aphi_dx = poloidal_A_b*( & + -poloidal_n_p*rho_dx/poloidal_rho_max*maxP_Pcut**poloidal_P_p & + +rhofac**poloidal_n_p*press_dx*poloidal_P_p*maxP_Pcut**(poloidal_P_p-1)) + Aphi_dy = poloidal_A_b*( & + -poloidal_n_p*rho_dy/poloidal_rho_max*maxP_Pcut**poloidal_P_p & + +rhofac**poloidal_n_p*press_dy*poloidal_P_p*maxP_Pcut**(poloidal_P_p-1)) + Aphi_dz = poloidal_A_b*( & + -poloidal_n_p*rho_dz/poloidal_rho_max*maxP_Pcut**poloidal_P_p & + +rhofac**poloidal_n_p*press_dz*poloidal_P_p*maxP_Pcut**(poloidal_P_p-1)) + endif Ax_dy = -AphiL - y(i,j,k)*Aphi_dy Ax_dz = -y(i,j,k)*Aphi_dz |