diff options
| -rw-r--r-- | interface.ccl | 2 | ||||
| -rw-r--r--[-rwxr-xr-x] | param.ccl | 66 | ||||
| -rw-r--r-- | schedule.ccl | 10 | ||||
| -rw-r--r-- | src/GRHydro_AdvectedLoopM.F90 | 20 | ||||
| -rw-r--r-- | src/GRHydro_AlfvenWaveM.F90 | 16 | ||||
| -rw-r--r-- | src/GRHydro_Bondi.c | 38 | ||||
| -rw-r--r--[-rwxr-xr-x] | src/GRHydro_Bondi_new.F90 | 0 |
7 files changed, 76 insertions, 76 deletions
diff --git a/interface.ccl b/interface.ccl index 06a497c..fce7f6c 100644 --- a/interface.ccl +++ b/interface.ccl @@ -151,7 +151,7 @@ void FUNCTION Con2PrimGenM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCT void FUNCTION Con2PrimPolyM(CCTK_INT INOUT handle, CCTK_REAL INOUT gamma_eos, CCTK_REAL INOUT dens, \ CCTK_REAL INOUT sx, CCTK_REAL INOUT sy, CCTK_REAL INOUT sz, \ CCTK_REAL INOUT sc, CCTK_REAL INOUT Bconsx, CCTK_REAL INOUT Bconsy, CCTK_REAL INOUT Bconsz, \ - CCTK_REAL INOUT rho, \ + CCTK_REAL INOUT rho, \ CCTK_REAL INOUT velx, CCTK_REAL INOUT vely, CCTK_REAL INOUT velz, \ CCTK_REAL INOUT epsilon, CCTK_REAL INOUT pressure, \ CCTK_REAL INOUT Bvecx, CCTK_REAL INOUT Bvecy, CCTK_REAL INOUT Bvecz, \ diff --git a/param.ccl b/param.ccl index 466fbe0..cf0a7ba 100755..100644 --- a/param.ccl +++ b/param.ccl @@ -10,9 +10,9 @@ USES KEYWORD temperature_evolution_method EXTENDS KEYWORD initial_hydro "" { - "shocktube" :: "Shocktube type" - "shocktube_hot" :: "Shocktube with hot nuclear EOS" - "only_atmo" :: "Set only a low atmosphere" + "shocktube" :: "Shocktube type" + "shocktube_hot" :: "Shocktube with hot nuclear EOS" + "only_atmo" :: "Set only a low atmosphere" "read_conformal":: "After reading in initial alp, rho and gxx from h5 files, sets the other quantities" "simple_wave" :: "Set initial data from Anile Miller Motta, Phys.Fluids. 26, 1450 (1983)" "monopole" :: "Monopole at the center" @@ -28,38 +28,38 @@ EXTENDS KEYWORD initial_hydro "" EXTENDS KEYWORD initial_Bvec { - "shocktube" :: "Shocktube type" + "shocktube" :: "Shocktube type" "poloidalmagfield" :: "Poloidal Magnetic Field" } shares:ADMBase EXTENDS KEYWORD initial_data "" { -# "shocktube" :: "Shock tube initial data for GRHydro" - "con2primtest" :: "Testing the con -> prim conversion" - "con2prim2con_test" :: "Testing the con -> prim -> con conversion" - "prim2con2prim_test" :: "Testing the prim -> con -> prim conversion" - "prim2con2prim_polytype_test" :: "Testing the prim -> con -> prim conversion - polytype version" - "reconstruction_test" :: "Testing reconstruction" +# "shocktube" :: "Shock tube initial data for GRHydro" + "con2primtest" :: "Testing the con -> prim conversion" + "con2prim2con_test" :: "Testing the con -> prim -> con conversion" + "prim2con2prim_test" :: "Testing the prim -> con -> prim conversion" + "prim2con2prim_polytype_test" :: "Testing the prim -> con -> prim conversion - polytype version" + "reconstruction_test" :: "Testing reconstruction" } private: KEYWORD shocktube_type "Diagonal or parallel shock?" { - "diagshock" :: "Diagonal across all axes" - "diagshock2d" :: "Diagonal across x-y axes" - "xshock" :: "Parallel to x axis" - "yshock" :: "Parallel to y axis" - "zshock" :: "Parallel to z axis" + "diagshock" :: "Diagonal across all axes" + "diagshock2d" :: "Diagonal across x-y axes" + "xshock" :: "Parallel to x axis" + "yshock" :: "Parallel to y axis" + "zshock" :: "Parallel to z axis" "sphere" :: "spherically symmetric shock" } "xshock" KEYWORD shock_case "Simple, Sod's problem or other?" { - "Simple" :: "GRAstro_Hydro test case" - "Sod" :: "Sod's problem" - "Blast" :: "Strong blast wave" + "Simple" :: "GRAstro_Hydro test case" + "Sod" :: "Sod's problem" + "Blast" :: "Strong blast wave" "Balsaralike1" :: "Hydro version of Balsara Test #1" "Balsara0" :: "Balsara Test #1, but unmagnetized" "Balsara1" :: "Balsara Test #1" @@ -81,17 +81,17 @@ KEYWORD shock_case "Simple, Sod's problem or other?" REAL shock_xpos "Position of shock plane: x" { - *:* :: "Anything" + *:* :: "Anything" } 0.0 REAL shock_ypos "Position of shock plane: y" { - *:* :: "Anything" + *:* :: "Anything" } 0.0 REAL shock_zpos "Position of shock plane: z" { - *:* :: "Anything" + *:* :: "Anything" } 0.0 REAL shock_radius "Radius of sperical shock" @@ -173,20 +173,20 @@ CCTK_REAL cyl_r_outer "Inner Radius" KEYWORD advectedloop_type "2-dimensional or 3-dimensional?" { - "2D" :: "2-dimensional (B^z=0)" - "3D" :: "3-dimensional (B^3=0, where B^3 || oblique cylinder axis." + "2D" :: "2-dimensional (B^z=0)" + "3D" :: "3-dimensional (B^3=0, where B^3 || oblique cylinder axis." } "2D" KEYWORD advectedloop_case "V^z=0 or not?" { - "V^z=0" :: "Useful to evaluate divB deviations" - "V^z/=0" :: "Useful to evaluate con2prim robustness in keeping V^z const." + "V^z=0" :: "Useful to evaluate divB deviations" + "V^z/=0" :: "Useful to evaluate con2prim robustness in keeping V^z const." } "V^z=0" KEYWORD advectedloop_delA "How to calculate B^i field from the potential A^b" { - "Exact" :: "Analytic, exact closed formula applied" - "Numeric" :: "Finite difference approximation of the derivatives applied" + "Exact" :: "Analytic, exact closed formula applied" + "Numeric" :: "Finite difference approximation of the derivatives applied" } "Exact" #################################3 @@ -194,8 +194,8 @@ KEYWORD advectedloop_delA "How to calculate B^i field from the potential A^b" ################################3 KEYWORD alfvenwave_type "1-dimensional or 2-dimensional?" { - "1D" :: "1-dimensional" - "2D" :: "2-dimensional (in x-y plane)" + "1D" :: "1-dimensional" + "2D" :: "2-dimensional (in x-y plane)" } "1D" CCTK_REAL alfvenwave_pressure "P_gas for the Alfven wave" @@ -279,22 +279,22 @@ CCTK_REAL bondi_bmag "B_0 parameter for magnetized Bondi" CCTK_REAL poloidal_A_b "Vector potential strength" { - *:* :: "Anything." + *:* :: "Anything." } 0.1 CCTK_INT poloidal_n_p "Vector potential strength" { - (0:* :: "Any positive integer." + (0:* :: "Any positive integer." } 3 CCTK_REAL poloidal_P_cut "Pressure used to confine the B field inside a star" { - (0:* :: "Anything positive." + (0:* :: "Anything positive." } 1.0e-8 CCTK_REAL poloidal_rho_max "Maximum initial density" { - (0:* :: "Anything positive." + (0:* :: "Anything positive." } 1.0e-3 # for the Magnetic Rotor test: diff --git a/schedule.ccl b/schedule.ccl index 7e5ab2c..c12093a 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -53,13 +53,13 @@ if (CCTK_Equals(initial_hydro,"shocktube")) { if(clean_divergence){ schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection { - LANG: Fortran + LANG: Fortran SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons, GRHydro::psidc } "Diagonal shock boundary conditions" } else { schedule GRHydro_Diagshock_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection { - LANG: Fortran + LANG: Fortran SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons } "Diagonal shock boundary conditions" } @@ -69,7 +69,7 @@ if (CCTK_Equals(initial_hydro,"shocktube")) { { schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection { - LANG: Fortran + LANG: Fortran } "2-D Diagonal shock boundary conditions" } @@ -78,13 +78,13 @@ if (CCTK_Equals(initial_hydro,"shocktube")) { # if(clean_divergence){ # schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection # { -# LANG: Fortran +# LANG: Fortran # SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons, GRHydro::psidc # } "2-D Diagonal shock boundary conditions" # } else { # schedule GRHydro_Diagshock2D_BoundaryM IN ApplyBCs AFTER BoundaryConditions AFTER Boundary::Boundary_ClearSelection # { -# LANG: Fortran +# LANG: Fortran # SYNC: GRHydro::dens, GRHydro::tau, GRHydro::scon, GRHydro::Bcons # } "2-D Diagonal shock boundary conditions" # } diff --git a/src/GRHydro_AdvectedLoopM.F90 b/src/GRHydro_AdvectedLoopM.F90 index 922e7f5..21acdb1 100644 --- a/src/GRHydro_AdvectedLoopM.F90 +++ b/src/GRHydro_AdvectedLoopM.F90 @@ -83,7 +83,7 @@ subroutine GRHydro_AdvectedLoopM(CCTK_ARGUMENTS) rhoval = 1.0d0 pressval = 3.0d0 - if (CCTK_EQUALS(advectedloop_type,"2D")) then + if (CCTK_EQUALS(advectedloop_type,"2D")) then !!$ Vx, Vy and Vz values: @@ -105,7 +105,7 @@ subroutine GRHydro_AdvectedLoopM(CCTK_ARGUMENTS) call CCTK_WARN(0,"V^z component case not recognized!") end if - else if (CCTK_EQUALS(advectedloop_type,"3D")) then + else if (CCTK_EQUALS(advectedloop_type,"3D")) then vxval=0.2d0*sqrt(2.0d0) vyval=0.2d0 @@ -149,7 +149,7 @@ subroutine GRHydro_AdvectedLoopM(CCTK_ARGUMENTS) 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) @@ -210,16 +210,16 @@ subroutine GRHydro_AdvectedLoopM(CCTK_ARGUMENTS) !!$ need to make x_d periodic! - if(x_d.gt.1.5*diaglength) then - x_d=x_d-2.0*diaglength + if(x_d.gt.1.5*diaglength) then + x_d=x_d-2.0*diaglength else if (x_d.gt.0.5*diaglength .and. x_d.lt.1.5*diaglength) then - x_d=x_d-diaglength - else if(x_d.lt.-1.5*diaglength) then - x_d=x_d+2.0*diaglength + x_d=x_d-diaglength + else if(x_d.lt.-1.5*diaglength) then + x_d=x_d+2.0*diaglength else if (x_d.lt.(-0.5*diaglength) .and. x_d.gt.(-1.5*diaglength)) then - x_d=x_d+diaglength + x_d=x_d+diaglength endif - + radius = sqrt(x_d**2+y_d**2) if (CCTK_EQUALS(advectedloop_delA,"Exact")) then diff --git a/src/GRHydro_AlfvenWaveM.F90 b/src/GRHydro_AlfvenWaveM.F90 index 23deb2a..53eb7b0 100644 --- a/src/GRHydro_AlfvenWaveM.F90 +++ b/src/GRHydro_AlfvenWaveM.F90 @@ -125,13 +125,13 @@ subroutine GRHydro_AlfvenWaveM(CCTK_ARGUMENTS) do j=1,ny do k=1,nz - rho(i,j,k)=rhoval + rho(i,j,k)=rhoval press(i,j,k)=pressval eps(i,j,k)=epsval if (CCTK_EQUALS(alfvenwave_type,"1D")) then - wnbr = 2.0d0*pi/range_x + wnbr = 2.0d0*pi/range_x velx(i,j,k)=vxval vely(i,j,k)=-valf*AA*cos(wnbr*x(i,j,k)) @@ -142,22 +142,22 @@ subroutine GRHydro_AlfvenWaveM(CCTK_ARGUMENTS) else if (CCTK_EQUALS(alfvenwave_type,"2D")) then - diaglength=range_x*range_y/range_d + diaglength=range_x*range_y/range_d range_d = sqrt(range_x**2+range_y**2) cos_theta = range_y/range_d sin_theta = range_x/range_d wnbr = 2.0d0*pi/diaglength - xnew = cos_theta*x(i,j,k)+sin_theta*y(i,j,k) + xnew = cos_theta*x(i,j,k)+sin_theta*y(i,j,k) - vparallel=vxval - vperp=-valf*AA*cos(wnbr*xnew) - velx(i,j,k)=vparallel*cos_theta-vperp*sin_theta + vparallel=vxval + vperp=-valf*AA*cos(wnbr*xnew) + velx(i,j,k)=vparallel*cos_theta-vperp*sin_theta vely(i,j,k)=vparallel*sin_theta+vperp*cos_theta velz(i,j,k)=-valf*AA*sin(wnbr*xnew) Bparallel=Bxval Bperp=AA*Bxval*cos(wnbr*xnew) - Bvecx(i,j,k)=Bparallel*cos_theta-Bperp*sin_theta + Bvecx(i,j,k)=Bparallel*cos_theta-Bperp*sin_theta Bvecy(i,j,k)=Bparallel*sin_theta+Bperp*cos_theta Bvecz(i,j,k)=AA*Bxval*sin(wnbr*xnew) diff --git a/src/GRHydro_Bondi.c b/src/GRHydro_Bondi.c index c657eec..e1b3d4c 100644 --- a/src/GRHydro_Bondi.c +++ b/src/GRHydro_Bondi.c @@ -350,7 +350,7 @@ static void setutcon(CCTK_REAL *vcon, CCTK_REAL gcov[][NDIM]) + gcov[2][2] * vcon[2] * vcon[2] + gcov[3][3] * vcon[3] * vcon[3] + 2.*( gcov[1][2] * vcon[1] * vcon[2] - + gcov[1][3] * vcon[1] * vcon[3] + + gcov[1][3] * vcon[1] * vcon[3] + gcov[2][3] * vcon[2] * vcon[3] ); c += 1. ; /* vector is timelike */ @@ -385,7 +385,7 @@ static void bl_gcov_func( CCTK_REAL *x, CCTK_REAL gcov[NDIM][NDIM]) r2 = r*r ; DD = 1. - 2.*M/r + asq/r2 ; mu = 1. + asq*cth*cth/r2 ; - + gcov[TT][TT] = -(1. - 2.*M/(r*mu)) ; gcov[TT][3] = -2.*M*a*s2/(r*mu) ; gcov[3][TT] = gcov[TT][3] ; @@ -456,17 +456,17 @@ static void set_bondi_parameters( CCTK_REAL M_in, CCTK_REAL Mdot_in, CCTK_REAL #if( LTRACE ) CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"\n#######################################################\n"); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"Bondi Solution Parameters1: \n"); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"------------------------- \n\n"); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"M = %28.20e Mdot = %28.20e rs = %28.20e \n",M,Mdot,rs); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"vs = %28.20e cs = %28.20e rhos = %28.20e \n",vs,cs,rhos); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"hs = %28.20e K = %28.20e Qdot = %28.20e \n",hs,K,Qdot); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"gam= %28.20e r_sol= %28.20e \n",gamma_eos, r_sol); - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"rs = : %28.20e \n", rs) ; - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"urs = : %28.20e \n", vs) ; - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"rhos = : %28.20e \n", rhos) ; - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"K = : %28.20e \n", K) ; - CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"#######################################################\n\n"); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"Bondi Solution Parameters1: \n"); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"------------------------- \n\n"); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"M = %28.20e Mdot = %28.20e rs = %28.20e \n",M,Mdot,rs); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"vs = %28.20e cs = %28.20e rhos = %28.20e \n",vs,cs,rhos); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"hs = %28.20e K = %28.20e Qdot = %28.20e \n",hs,K,Qdot); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"gam= %28.20e r_sol= %28.20e \n",gamma_eos, r_sol); + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"rs = : %28.20e \n", rs) ; + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"urs = : %28.20e \n", vs) ; + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"rhos = : %28.20e \n", rhos) ; + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"K = : %28.20e \n", K) ; + CCTK_VWarn(CCTK_WARN_DEBUG, __LINE__, __FILE__, CCTK_THORNSTRING,"#######################################################\n\n"); #endif return; @@ -495,9 +495,9 @@ static void set_bondi_parameters( CCTK_REAL M_in, CCTK_REAL Mdot_in, CCTK_REAL *****************************************************************/ static int gnr_bondi( CCTK_REAL x[], int n, - void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], - CCTK_REAL [][NEWT_DIM_B], CCTK_REAL *, - CCTK_REAL *, int) ) + void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], + CCTK_REAL [][NEWT_DIM_B], CCTK_REAL *, + CCTK_REAL *, int) ) { CCTK_REAL f, df, dx[NEWT_DIM_B], x_old[NEWT_DIM_B], resid[NEWT_DIM_B], jac[NEWT_DIM_B][NEWT_DIM_B]; @@ -560,7 +560,7 @@ static int gnr_bondi( CCTK_REAL x[], int n, if( doing_extra == 1 ) i_extra++ ; if( ((fabs(errx) <= NEWT_TOL_B)&&(doing_extra == 0)) || - (i_extra > EXTRA_NEWT_ITER_B) || (n_iter >= (MAX_NEWT_ITER_B-1)) ) { + (i_extra > EXTRA_NEWT_ITER_B) || (n_iter >= (MAX_NEWT_ITER_B-1)) ) { keep_iterating = 0; } @@ -704,8 +704,8 @@ static int find_bondi_solution( CCTK_REAL r, CCTK_REAL *rho, CCTK_REAL *u, CCTK_ retval = gnr_bondi( rho, NEWT_DIM_B, bondi_resid); if( retval ) { - CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,"find_bondi_solution: Incr. guess failed, decrease dfactor, retval = %d, r = %g, itry = %d \n", retval, r, itry); - return(10); + CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,"find_bondi_solution: Incr. guess failed, decrease dfactor, retval = %d, r = %g, itry = %d \n", retval, r, itry); + return(10); } } diff --git a/src/GRHydro_Bondi_new.F90 b/src/GRHydro_Bondi_new.F90 index 7bb5548..7bb5548 100755..100644 --- a/src/GRHydro_Bondi_new.F90 +++ b/src/GRHydro_Bondi_new.F90 |