Commits

commits@lists.einsteintoolkit.org

2 participants
6997 discussions

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 640)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:50 PM Added: /trunk/src/ GRHydro_Prim2Con.cc Removed: /trunk/src/ GRHydro_Prim2Con.c Modified: /trunk/src/ GRHydro_Functions.h, GRHydro_Reconstruct.F90, GRHydro_Reconstruct.cc, make.code.defn Log: GRHydro: Templated Prim2Con version and related changes in various places. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Functions.h Delta lines: <Unavailable> File [removed]: GRHydro_Prim2Con.c Delta lines: <Unavailable> File [added]: GRHydro_Prim2Con.cc Delta lines: <Unavailable> File [modified]: GRHydro_Reconstruct.F90 Delta lines: <Unavailable> File [modified]: GRHydro_Reconstruct.cc Delta lines: <Unavailable> File [modified]: make.code.defn Delta lines: <Unavailable>

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 639)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:50 PM Modified: /trunk/src/ GRHydro_Con2PrimM_pt_EOSOmni.c Log: GRhydro: detect EOS errors in MHD Con2Prim pressure call this is a cut down version of 3e18d28da06966fddcaf0612afa1006ceb027d18 which also contains vector potential changes. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM_pt_EOSOmni.c Delta lines: +6 -0 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2014-04-15 19:50:03 UTC (rev 638) +++ trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2014-04-15 19:50:06 UTC (rev 639) @@ -397,6 +397,12 @@ // eosvars.eosanyerr[0]=5; p = pressure_rho0_eps_eosomni(rho0,*epsilon,&dum1,&dum2,&eosvars) ; // EOSOMNI // p = pressure_rho0_u(rho0,u,gammaeos) ; // EOS + // error handling + if (eosvars.eosanyerr[0] != 0) { + if(eosvars.eoskeyerr[0] == 668) { + *retval = 3.0; + } + } w = rho0 + u + p;

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 638)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:50 PM Modified: /trunk/src/ GRHydro_Con2PrimM.F90, GRHydro_Con2PrimM_pt_EOSOmni.c, GRHydro_Prim2ConM.F90 Log: GRHydro: disable warnings in Prim2ConM and only set mask in Con2PrimM We don't abort in ptwise C2P call with MHD so that we can check for the mask and abort if necessary. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Con2PrimM.F90 Delta lines: +1 -1 =================================================================== --- trunk/src/GRHydro_Con2PrimM.F90 2014-04-15 19:50:00 UTC (rev 637) +++ trunk/src/GRHydro_Con2PrimM.F90 2014-04-15 19:50:03 UTC (rev 638) @@ -542,7 +542,7 @@ write(warnline,"(1P10E15.6)") rho(i,j,k),dens(i,j,k),eps(i,j,k),& temperature(i,j,k),y_e(i,j,k) call CCTK_WARN(1,warnline) - call CCTK_WARN(0,"Aborting!!!") + ! call CCTK_WARN(0,"Aborting!!!") endif !$OMP END CRITICAL endif File [modified]: GRHydro_Con2PrimM_pt_EOSOmni.c Delta lines: +40 -38 =================================================================== --- trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2014-04-15 19:50:00 UTC (rev 637) +++ trunk/src/GRHydro_Con2PrimM_pt_EOSOmni.c 2014-04-15 19:50:03 UTC (rev 638) @@ -504,44 +504,46 @@ /* Problem with solver, so return denoting error before doing anything further */ if( ((*retval) != 0.) || (W == FAIL_VAL) ) { *retval = *retval*100.+1.; - fprintf(stdout," *retval = %26.16e \n", *retval ); - fprintf(stdout," *dens = %26.16e \n", *dens_in ); - fprintf(stdout," *sx = %26.16e \n", *sx_in ); - fprintf(stdout," *sy = %26.16e \n", *sy_in ); - fprintf(stdout," *sz = %26.16e \n", *sz_in ); - fprintf(stdout," *tau = %26.16e \n", *tau_in ); - fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); - fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); - fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); - fprintf(stdout," *rho = %26.16e \n", *rho ); - fprintf(stdout," *velx = %26.16e \n", *velx ); - fprintf(stdout," *vely = %26.16e \n", *vely ); - fprintf(stdout," *velz = %26.16e \n", *velz ); - fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); - fprintf(stdout," *temp_in = %26.16e \n", *temp_in ); - fprintf(stdout," *y_e_in = %26.16e \n", *y_e_in ); - fprintf(stdout," *pressure = %26.16e \n", *pressure ); - fprintf(stdout," *Bx = %26.16e \n", *Bx ); - fprintf(stdout," *By = %26.16e \n", *By ); - fprintf(stdout," *Bz = %26.16e \n", *Bz ); - fprintf(stdout," *bsq = %26.16e \n", *bsq ); - fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); - fprintf(stdout," *gxx = %26.16e \n", *gxx ); - fprintf(stdout," *gxy = %26.16e \n", *gxy ); - fprintf(stdout," *gxz = %26.16e \n", *gxz ); - fprintf(stdout," *gyy = %26.16e \n", *gyy ); - fprintf(stdout," *gyz = %26.16e \n", *gyz ); - fprintf(stdout," *gzz = %26.16e \n", *gzz ); - fprintf(stdout," *uxx = %26.16e \n", *uxx ); - fprintf(stdout," *uxy = %26.16e \n", *uxy ); - fprintf(stdout," *uxz = %26.16e \n", *uxz ); - fprintf(stdout," *uyy = %26.16e \n", *uyy ); - fprintf(stdout," *uyz = %26.16e \n", *uyz ); - fprintf(stdout," *uzz = %26.16e \n", *uzz ); - fprintf(stdout," *sdet = %26.16e \n", *sdet ); - fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); - fprintf(stdout," *retval = %26.16e \n", *retval ); - fflush(stdout); +#if(DEBUG_CON2PRIMM) + fprintf(stdout," *retval = %26.16e \n", *retval ); + fprintf(stdout," *dens = %26.16e \n", *dens_in ); + fprintf(stdout," *sx = %26.16e \n", *sx_in ); + fprintf(stdout," *sy = %26.16e \n", *sy_in ); + fprintf(stdout," *sz = %26.16e \n", *sz_in ); + fprintf(stdout," *tau = %26.16e \n", *tau_in ); + fprintf(stdout," *Bconsx = %26.16e \n", *Bconsx_in ); + fprintf(stdout," *Bconsy = %26.16e \n", *Bconsy_in ); + fprintf(stdout," *Bconsz = %26.16e \n", *Bconsz_in ); + fprintf(stdout," *rho = %26.16e \n", *rho ); + fprintf(stdout," *velx = %26.16e \n", *velx ); + fprintf(stdout," *vely = %26.16e \n", *vely ); + fprintf(stdout," *velz = %26.16e \n", *velz ); + fprintf(stdout," *epsilon = %26.16e \n", *epsilon ); + fprintf(stdout," *temp_in = %26.16e \n", *temp_in ); + fprintf(stdout," *y_e_in = %26.16e \n", *y_e_in ); + fprintf(stdout," *pressure = %26.16e \n", *pressure ); + fprintf(stdout," *Bx = %26.16e \n", *Bx ); + fprintf(stdout," *By = %26.16e \n", *By ); + fprintf(stdout," *Bz = %26.16e \n", *Bz ); + fprintf(stdout," *bsq = %26.16e \n", *bsq ); + fprintf(stdout," *w_lorentz = %26.16e \n", *w_lorentz ); + fprintf(stdout," *gxx = %26.16e \n", *gxx ); + fprintf(stdout," *gxy = %26.16e \n", *gxy ); + fprintf(stdout," *gxz = %26.16e \n", *gxz ); + fprintf(stdout," *gyy = %26.16e \n", *gyy ); + fprintf(stdout," *gyz = %26.16e \n", *gyz ); + fprintf(stdout," *gzz = %26.16e \n", *gzz ); + fprintf(stdout," *uxx = %26.16e \n", *uxx ); + fprintf(stdout," *uxy = %26.16e \n", *uxy ); + fprintf(stdout," *uxz = %26.16e \n", *uxz ); + fprintf(stdout," *uyy = %26.16e \n", *uyy ); + fprintf(stdout," *uyz = %26.16e \n", *uyz ); + fprintf(stdout," *uzz = %26.16e \n", *uzz ); + fprintf(stdout," *sdet = %26.16e \n", *sdet ); + fprintf(stdout," *epsnegative = %10d \n", *epsnegative ); + fprintf(stdout," *retval = %26.16e \n", *retval ); + fflush(stdout); +#endif return; } else{ File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +0 -24 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:50:00 UTC (rev 637) +++ trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:50:03 UTC (rev 638) @@ -289,34 +289,10 @@ call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) keytemp=reconstruct_temper - if(anyerr.ne.0) then - !$OMP CRITICAL - call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(1,warnline) - write(warnline,"(1P10E15.6)") drho,deps,temp,ye,yein - call CCTK_WARN(1,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(1,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(1,warnline) - !$OMP END CRITICAL - endif else ! This is a way of recovering even on finer refinement levels: ! Use the average temperature at the interface instead of the ! reconstructed specific internal energy. - !$OMP CRITICAL - call CCTK_WARN(1,"EOS error in prim2con_hot: NOW using averaged temp!") - write(warnline,"(3i5,1P10E15.6)") ii,jj,kk,x,y,z - call CCTK_WARN(1,warnline) - write(warnline,"(1P10E15.6)") drho,deps,temp0,ye - call CCTK_WARN(1,warnline) - write(warnline,"(A7,i8)") "code: ",keyerr(1) - call CCTK_WARN(1,warnline) - write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel - call CCTK_WARN(1,warnline) - !$OMP END CRITICAL keytemp=1-reconstruct_temper temp = temp0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,&

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 637)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:50 PM Modified: /trunk/src/ GRHydro_Boundaries.F90 Log: GRHydro: clean up sync_conservatives_only primitives bc calls Define helper routine to do all boundary registration for a single bc type. This avoids code duplication and the possibility that eg flat and none are applied to different sets of variables. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: +114 -182 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:58 UTC (rev 636) +++ trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:50:00 UTC (rev 637) @@ -296,13 +296,15 @@ -subroutine GRHydro_SelectPrimitiveInitialGuessesBoundaries(CCTK_ARGUMENTS) +subroutine GRHydro_SelectPrimitiveInitialGuessesBoundariesForMethod(CCTK_ARGUMENTS, bcname) implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + + character(len=*), INTENT(IN) :: bcname integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates @@ -314,116 +316,80 @@ general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) -!!$Flat boundaries if required - ! The commented out code are those primitives which do not require and explicit ! initial guess for Con2Prim. The guesses are either not needed or computed from other quantities! - if (verbose.eq.1) call CCTK_INFO("Selecting primitive BC") + !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + ! "HydroBase::w_lorentz", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::rho", bcname) + !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + ! "HydroBase::press", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::eps", bcname) + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lvel", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", bcname) + endif + if(evolve_mhd.ne.0) then + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", bcname) + endif + endif + !if(evolve_tracer.ne.0) then + ! ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + ! "GRHydro::GRHydro_tracers", bcname) + !endif - if (CCTK_EQUALS(bound,"flat")) then - - !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::w_lorentz", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "Flat") - !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::press", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "Flat") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "Flat") - endif - if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif - endif + !if(evolve_y_e.ne.0) then + ! ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + ! "HydroBase::Y_e", bcname) + !endif - !if(evolve_tracer.ne.0) then - ! ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "GRHydro::GRHydro_tracers", "Flat") - !endif + if(evolve_temper.ne.0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::temperature", bcname) + !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + ! "HydroBase::entropy", bcname) + endif - !if(evolve_y_e.ne.0) then - ! ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::Y_e", "Flat") - !endif + if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - if(evolve_temper.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "Flat") - !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::entropy", "Flat") - endif +end subroutine GRHydro_SelectPrimitiveInitialGuessesBoundariesForMethod +subroutine GRHydro_SelectPrimitiveInitialGuessesBoundaries(CCTK_ARGUMENTS) - endif + implicit none - if (CCTK_EQUALS(bound,"none")) then - !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::w_lorentz", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "None") - !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::press", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "None") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "None") - endif - - if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif - - endif + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS - !if(evolve_tracer.ne.0) then - ! ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "GRHydro::GRHydro_tracers", "None") - ! - !endif + if (verbose.eq.1) call CCTK_INFO("Selecting primitive BC") - !if(evolve_y_e.ne.0) then - ! ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::Y_e", "None") - ! - !endif +!!$Flat boundaries if required + if (CCTK_EQUALS(bound,"flat")) then + call GRHydro_SelectPrimitiveInitialGuessesBoundariesForMethod(CCTK_PASS_FTOF, "Flat") + endif - if(evolve_temper.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "None") - !ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - ! "HydroBase::entropy", "None") - endif - +!!$None boundaries if required + if (CCTK_EQUALS(bound,"none")) then + call GRHydro_SelectPrimitiveInitialGuessesBoundariesForMethod(CCTK_PASS_FTOF, "None") end if +!!$Scalar boundaries if required if (CCTK_EQUALS(bound,"scalar")) then call CCTK_WARN(0, "Until somebody uses this I see no reason to support it") end if - if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - end subroutine GRHydro_SelectPrimitiveInitialGuessesBoundaries @@ -433,7 +399,7 @@ -subroutine GRHydro_SelectPrimitiveBoundaries(CCTK_ARGUMENTS) +subroutine GRHydro_SelectPrimitiveBoundariesForMethod(CCTK_ARGUMENTS, bcname) implicit none @@ -441,6 +407,8 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + character(len=*), INTENT(IN) :: bcname + integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates @@ -451,113 +419,77 @@ general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) -!!$Flat boundaries if required + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::w_lorentz", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::rho", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::press", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::eps", bcname) + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lvel", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", bcname) + endif + if(evolve_mhd.ne.0) then + if (general_coordinates .ne. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::lBvec", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", bcname) + endif + endif - if (verbose.eq.1) call CCTK_INFO("Selecting primitive BC") + if(evolve_tracer.ne.0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_tracers", bcname) + endif + if(evolve_y_e.ne.0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Y_e", bcname) + endif - if (CCTK_EQUALS(bound,"flat")) then - - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "Flat") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "Flat") - endif - if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif - endif + if(evolve_temper.ne.0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::temperature", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::entropy", bcname) + endif - if(evolve_tracer.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "Flat") - endif + if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - if(evolve_y_e.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "Flat") - endif +end subroutine GRHydro_SelectPrimitiveBoundariesForMethod - if(evolve_temper.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "Flat") - endif +subroutine GRHydro_SelectPrimitiveBoundaries(CCTK_ARGUMENTS) + implicit none - endif + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS - if (CCTK_EQUALS(bound,"none")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "None") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "None") - endif - - if(evolve_mhd.ne.0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif - - endif + if (verbose.eq.1) call CCTK_INFO("Selecting primitive BC") - if(evolve_tracer.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "None") - - endif +!!$Flat boundaries if required + if (CCTK_EQUALS(bound,"flat")) then + call GRHydro_SelectPrimitiveBoundariesForMethod(CCTK_ARGUMENTS, "Flat") + endif - if(evolve_y_e.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "None") - - endif - - if(evolve_temper.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "None") - endif - +!!$None boundaries if required + if (CCTK_EQUALS(bound,"none")) then + call GRHydro_SelectPrimitiveBoundariesForMethod(CCTK_ARGUMENTS, "None") end if +!!$Scalar boundaries if required if (CCTK_EQUALS(bound,"scalar")) then call CCTK_WARN(0, "Until somebody uses this I see no reason to support it") end if - if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - end subroutine GRHydro_SelectPrimitiveBoundaries

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 636)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_Boundaries.F90 Log: GRHydro: remove unused variables in boundary routine File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Boundaries.F90 Delta lines: +0 -10 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:55 UTC (rev 635) +++ trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:58 UTC (rev 636) @@ -175,7 +175,6 @@ character(len=*), INTENT(IN) :: bcname - integer, dimension(3) :: sw integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates @@ -185,8 +184,6 @@ ierr = 0 - sw = GRHydro_stencil - general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) do_Avec = CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") do_Bcons = CCTK_EQUALS(Bvec_evolution_method,"GRHydro") @@ -307,7 +304,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer, dimension(3) :: sw integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates @@ -315,7 +311,6 @@ CCTK_INT, parameter :: ione=1 ierr = 0 - sw = GRHydro_stencil general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) @@ -446,7 +441,6 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer, dimension(3) :: sw integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates @@ -454,7 +448,6 @@ CCTK_INT, parameter :: ione=1 ierr = 0 - sw = GRHydro_stencil general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) @@ -582,16 +575,13 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer, dimension(3) :: sw integer :: ierr CCTK_INT, parameter :: faces=CCTK_ALL_FACES CCTK_INT, parameter :: ione=1 ierr = 0 - sw = GRHydro_stencil - !!$Flat boundaries if required if (verbose.eq.1) call CCTK_INFO("Selecting atmosphere mask BC")

1 0

[svn:einsteintoolkit] GRHydro/trunk/ (Rev. 635)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/ schedule.ccl /trunk/src/ GRHydro_Analysis.F90, GRHydro_Boundaries.F90, GRHydro_CalcUpdate.F90, GRHydro_UpdateMask.F90 Log: GRHydro: clean up handling of Avec evolution and boundaries Remove almost empty Analysis_Init routine, apply boundary conditions consistently. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Analysis.F90 Delta lines: +70 -55 =================================================================== --- trunk/src/GRHydro_Analysis.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_Analysis.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -5,18 +5,6 @@ #include "cctk_Functions.h" #include "SpaceMask.h" -subroutine GRHydro_Analysis_Init(CCTK_ARGUMENTS) - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS - - divB = 0.0d0 - -end subroutine GRHydro_Analysis_Init - subroutine GRHydro_CalcDivB(CCTK_ARGUMENTS) implicit none @@ -27,56 +15,83 @@ CCTK_INT :: i,j,k,itracer CCTK_REAL :: idx, idy, idz, Bcons_l1, Bcons_r1, Bcons_l2, Bcons_r2, Bcons_l3, Bcons_r3 + logical :: do_Avec, do_Bcons + idx = 1.d0 / CCTK_DELTA_SPACE(1) idy = 1.d0 / CCTK_DELTA_SPACE(2) idz = 1.d0 / CCTK_DELTA_SPACE(3) - !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,Bcons_l1,Bcons_r1,Bcons_l2,Bcons_r2,Bcons_l3,Bcons_r3) - do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined - do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil - do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + do_Avec = CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") + do_Bcons = CCTK_EQUALS(Bvec_evolution_method,"GRHydro") - if(evolve_mhd.ne.0) then - if(track_divB.ne.0) then - if(transport_constraints.ne.0) then - ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) - divB(i,j,k) = & - 0.25d0*(Bcons(i+1,j,k,1)-Bcons(i,j,k,1)+ & - Bcons(i+1,j+1,k,1)-Bcons(i,j+1,k,1)+ & - Bcons(i+1,j,k+1,1)-Bcons(i,j,k+1,1)+ & - Bcons(i+1,j+1,k+1,1)-Bcons(i,j+1,k+1,1))*idx + & - 0.25d0*(Bcons(i,j+1,k,2)-Bcons(i,j,k,2)+ & - Bcons(i+1,j+1,k,2)-Bcons(i+1,j,k,2)+ & - Bcons(i,j+1,k+1,2)-Bcons(i,j,k+1,2)+ & - Bcons(i+1,j+1,k+1,2)-Bcons(i+1,j,k+1,2))*idy + & - 0.25d0*(Bcons(i,j,k+1,3)-Bcons(i,j,k,3)+ & - Bcons(i+1,j,k+1,3)-Bcons(i+1,j,k,3)+ & - Bcons(i,j+1,k+1,3)-Bcons(i,j+1,k,3)+ & - Bcons(i+1,j+1,k+1,3)-Bcons(i+1,j+1,k,3))*idz - else - Bcons_l1 = 0.5d0 * (Bcons(i,j,k,1) + & - Bcons(i-1,j,k,1)) - Bcons_l2 = 0.5d0 * (Bcons(i,j,k,2) + & - Bcons(i,j-1,k,2)) - Bcons_l3 = 0.5d0 * (Bcons(i,j,k,3) + & - Bcons(i,j,k-1,3)) - Bcons_r1 = 0.5d0 * (Bcons(i,j,k,1) + & - Bcons(i+1,j,k,1)) - Bcons_r2 = 0.5d0 * (Bcons(i,j,k,2) + & - Bcons(i,j+1,k,2)) - Bcons_r3 = 0.5d0 * (Bcons(i,j,k,3) + & - Bcons(i,j,k+1,3)) - - divB(i,j,k) = (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + divB = 0d0 - endif - endif - endif - - enddo - enddo + !$OMP PARALLEL DO PRIVATE(i,j,k,itracer,Bcons_l1,Bcons_r1,Bcons_l2,Bcons_r2,Bcons_l3,Bcons_r3) + do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined + do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil + do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil + + if(evolve_mhd.ne.0) then + if(track_divB.ne.0) then + if(transport_constraints.ne.0) then + ! edge based divergence (see WhiskyMHD & Bruno's thesis, Eq. 7.27) + divB(i,j,k) = & + 0.25d0*(Bcons(i+1,j,k,1)-Bcons(i,j,k,1)+ & + Bcons(i+1,j+1,k,1)-Bcons(i,j+1,k,1)+ & + Bcons(i+1,j,k+1,1)-Bcons(i,j,k+1,1)+ & + Bcons(i+1,j+1,k+1,1)-Bcons(i,j+1,k+1,1))*idx + & + 0.25d0*(Bcons(i,j+1,k,2)-Bcons(i,j,k,2)+ & + Bcons(i+1,j+1,k,2)-Bcons(i+1,j,k,2)+ & + Bcons(i,j+1,k+1,2)-Bcons(i,j,k+1,2)+ & + Bcons(i+1,j+1,k+1,2)-Bcons(i+1,j,k+1,2))*idy + & + 0.25d0*(Bcons(i,j,k+1,3)-Bcons(i,j,k,3)+ & + Bcons(i+1,j,k+1,3)-Bcons(i+1,j,k,3)+ & + Bcons(i,j+1,k+1,3)-Bcons(i,j+1,k,3)+ & + Bcons(i+1,j+1,k+1,3)-Bcons(i+1,j+1,k,3))*idz + else if(do_Avec) then + Bcons_l1 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,1) + & + sdetg(i-1,j,k)*Bvec(i-1,j,k,1)) + Bcons_l2 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,2) + & + sdetg(i,j-1,k)*Bvec(i,j-1,k,2)) + Bcons_l3 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,3) + & + sdetg(i,j,k-1)*Bvec(i,j,k-1,3)) + Bcons_r1 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,1) + & + sdetg(i+1,j,k)*Bvec(i+1,j,k,1)) + Bcons_r2 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,2) + & + sdetg(i,j+1,k)*Bvec(i,j+1,k,2)) + Bcons_r3 = 0.5d0 * (sdetg(i,j,k)*Bvec(i,j,k,3) + & + sdetg(i,j,k+1)*Bvec(i,j,k+1,3)) + + divB(i,j,k) = (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + + else if(do_Bcons) then + Bcons_l1 = 0.5d0 * (Bcons(i,j,k,1) + & + Bcons(i-1,j,k,1)) + Bcons_l2 = 0.5d0 * (Bcons(i,j,k,2) + & + Bcons(i,j-1,k,2)) + Bcons_l3 = 0.5d0 * (Bcons(i,j,k,3) + & + Bcons(i,j,k-1,3)) + Bcons_r1 = 0.5d0 * (Bcons(i,j,k,1) + & + Bcons(i+1,j,k,1)) + Bcons_r2 = 0.5d0 * (Bcons(i,j,k,2) + & + Bcons(i,j+1,k,2)) + Bcons_r3 = 0.5d0 * (Bcons(i,j,k,3) + & + Bcons(i,j,k+1,3)) + + divB(i,j,k) = (Bcons_l1 - Bcons_r1) * idx + (Bcons_l2 - Bcons_r2) * idy + (Bcons_l3 - Bcons_r3) * idz + + else + call CCTK_ERROR("Internal error. Do now know how to compute divB") + STOP + + endif + endif + endif + enddo - !$OMP END PARALLEL DO + enddo + enddo + !$OMP END PARALLEL DO end subroutine GRHydro_CalcDivB File [modified]: GRHydro_Boundaries.F90 Delta lines: +86 -149 =================================================================== --- trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_Boundaries.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -165,199 +165,136 @@ @@*/ -subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) - +subroutine GRHydro_RegisterBoundariesWithMethod(CCTK_ARGUMENTS, bcname) + implicit none - + DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - + + character(len=*), INTENT(IN) :: bcname + integer, dimension(3) :: sw integer :: ierr CCTK_INT :: GRHydro_UseGeneralCoordinates, general_coordinates CCTK_INT, parameter :: faces=CCTK_ALL_FACES CCTK_INT, parameter :: ione=1 + logical :: do_Avec, do_Bcons ierr = 0 sw = GRHydro_stencil general_coordinates = GRHydro_UseGeneralCoordinates(cctkGH) + do_Avec = CCTK_EQUALS(Bvec_evolution_method,"GRHydro_Avec") + do_Bcons = CCTK_EQUALS(Bvec_evolution_method,"GRHydro") - if (verbose.eq.1) call CCTK_INFO("Selecting conserved BC (and primitive BC if selected)") - -!!$Flat boundaries if required - - if (CCTK_EQUALS(bound,"flat")) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::dens", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::tau", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::scon", bcname) + if (sync_conserved_only .eq. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::dens", "Flat") + "HydroBase::w_lorentz", bcname) ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::tau", "Flat") + "HydroBase::rho", bcname) ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::scon", "Flat") - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "Flat") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "Flat") - endif - endif - if(evolve_mhd.ne.0) then - if (sync_conserved_only .eq. 0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif - endif + "HydroBase::press", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::eps", bcname) + if (general_coordinates .ne. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Bcons", "Flat") - if(clean_divergence.ne.0) then + "GRHydro::lvel", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::vel", bcname) + endif + endif + if(evolve_mhd.ne.0) then + if (sync_conserved_only .eq. 0) then + if (general_coordinates .ne. 0) then ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::psidc", "Flat") + "GRHydro::lBvec", bcname) + else + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Bvec", bcname) endif - endif + endif + if(do_Bcons) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::Bcons", bcname) + end if + if (do_Avec) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Avec", bcname) + endif + if(clean_divergence.ne.0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::psidc", bcname) + endif + endif - if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Avec", "Flat") - endif - if(evolve_tracer.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "Flat") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_cons_tracers", "Flat") + if(evolve_tracer.ne.0) then + if (sync_conserved_only .eq. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_tracers", bcname) endif + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::GRHydro_cons_tracers", bcname) + endif - if(evolve_y_e.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "Flat") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Y_e_con", "Flat") + if(evolve_y_e.ne.0) then + if (sync_conserved_only .eq. 0) then + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::Y_e", bcname) endif - - if(evolve_temper.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "Flat") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "Flat") - endif - endif - - + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "GRHydro::Y_e_con", bcname) endif - if (CCTK_EQUALS(bound,"none")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::dens", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::tau", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::scon", "None") - + if(evolve_temper.ne.0) then if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::w_lorentz", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::rho", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::press", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::eps", "None") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lvel", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::vel", "None") - endif - endif - - if(evolve_mhd.ne.0) then - if (sync_conserved_only .eq. 0) then - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "None") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "None") - endif - endif ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Bcons", "None") - if(clean_divergence.ne.0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::psidc", "None") - endif + "HydroBase::temperature", bcname) + ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & + "HydroBase::entropy", bcname) endif + endif - if (CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec")) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Avec", "None") - if (general_coordinates .ne. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::lBvec", "Flat") - else - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Bvec", "Flat") - endif - endif + if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - if(evolve_tracer.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_tracers", "None") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::GRHydro_cons_tracers", "None") - endif +end subroutine - if(evolve_y_e.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::Y_e", "None") - endif - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "GRHydro::Y_e_con", "None") - endif +subroutine GRHydro_Boundaries(CCTK_ARGUMENTS) - if(evolve_temper.ne.0) then - if (sync_conserved_only .eq. 0) then - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::temperature", "None") - ierr = ierr + Boundary_SelectGroupForBC(cctkGH, faces, GRHydro_stencil, -ione, & - "HydroBase::entropy", "None") - endif - endif + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + if (verbose.eq.1) call CCTK_INFO("Selecting conserved BC (and primitive BC if selected)") + +!!$Flat boundaries if required + if (CCTK_EQUALS(bound,"flat")) then + call GRHydro_RegisterBoundariesWithMethod(CCTK_PASS_FTOF, "Flat") endif +!!$None boundaries if required + if (CCTK_EQUALS(bound,"none")) then + call GRHydro_RegisterBoundariesWithMethod(CCTK_PASS_FTOF, "None") + endif + +!!$Scalar boundaries if required if (CCTK_EQUALS(bound,"scalar")) then call CCTK_WARN(0, "Until somebody uses this I see no reason to support it") end if - if (ierr < 0) call CCTK_WARN(0, "problems with applying the chosen boundary condition") - end subroutine GRHydro_Boundaries File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +16 -17 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_CalcUpdate.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -108,24 +108,23 @@ taurhs(i,j,k) = taurhs(i,j,k) + & (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - & alp_r * tauflux(i,j,k)) * idx - if(evolve_mhd.ne.0) then - if(transport_constraints.ne.1) then - Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & - (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconsxflux(i,j,k)) * idx - Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & - (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconsyflux(i,j,k)) * idx - Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & - (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * Bconszflux(i,j,k)) * idx - endif - if(clean_divergence.ne.0) then - psidcrhs(i,j,k) = psidcrhs(i,j,k) + & - (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & - alp_r * psidcflux(i,j,k)) * idx - endif + if(transport_constraints.eq.0) then + Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + & + (alp_l * Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsxflux(i,j,k)) * idx + Bconsrhs(i,j,k,2) = Bconsrhs(i,j,k,2) + & + (alp_l * Bconsyflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconsyflux(i,j,k)) * idx + Bconsrhs(i,j,k,3) = Bconsrhs(i,j,k,3) + & + (alp_l * Bconszflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * Bconszflux(i,j,k)) * idx + endif + if(clean_divergence.ne.0) then + psidcrhs(i,j,k) = psidcrhs(i,j,k) + & + (alp_l * psidcflux(i-xoffset,j-yoffset,k-zoffset) - & + alp_r * psidcflux(i,j,k)) * idx + endif endif if (evolve_tracer .ne. 0) then File [modified]: GRHydro_UpdateMask.F90 Delta lines: +6 -0 =================================================================== --- trunk/src/GRHydro_UpdateMask.F90 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/src/GRHydro_UpdateMask.F90 2014-04-15 19:49:55 UTC (rev 635) @@ -65,6 +65,9 @@ if (evolve_Bvec) then Bconsrhs(i,j,k,:) = 0.0d0 endif + if (evolve_Avec) then + Avecrhs(i,j,k,:) = 0.0d0 + endif ! TODO: Need to set Avecrhs and Aphirhs to zero for centered / or staggered case if vector potential is evolved @@ -93,6 +96,9 @@ if (evolve_Bvec) then Bconsrhs(i,j,k,:) = 0.0d0 endif + if (evolve_Avec) then + Avecrhs(i,j,k,:) = 0.0d0 + endif ! TODO: Need to set Avecrhs and Aphirhs to zero for centered / or staggered case if vector potential is evolved Directory: /trunk/ ================== File [modified]: schedule.ccl Delta lines: +1 -7 =================================================================== --- trunk/schedule.ccl 2014-04-15 19:49:52 UTC (rev 634) +++ trunk/schedule.ccl 2014-04-15 19:49:55 UTC (rev 635) @@ -80,7 +80,6 @@ } else if(CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { STORAGE: HydroBase::Bvec[timelevels] STORAGE: HydroBase::Avec[timelevels] - STORAGE: GRHydro::Bcons[timelevels] if (CCTK_Equals(Avec_gauge,"Lorenz")) { STORAGE: HydroBase::Aphi[timelevels] @@ -1611,16 +1610,11 @@ { } "Calculate analysis quantities" - schedule GRHydro_Analysis_Init IN GRHydroAnalysis - { - LANG: Fortran - } "Initialize variables" - schedule GRHydro_CalcDivB IN GRHydroAnalysis AFTER GRHydro_Analysis_Init { LANG: Fortran } "Calculate divB" -} +} if (constrain_to_1D) {

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 634)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_Prim2ConM.F90 Log: GRHydro: use temperature in MHD prim2con this uses tempplus/tempminus in MHD prim2con this helps when temperature is actually reconstructed. It *does* however change the results already when temperature is not reconstructed, though it does make the results use constant values rather then averages which is lower order. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +7 -22 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:49:50 UTC (rev 633) +++ trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:49:52 UTC (rev 634) @@ -176,14 +176,6 @@ avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l,g23l,g33l)) avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r,g23r,g33r)) - ! we do not have plus/minus vars for temperature since - ! eps is reconstructed. Hence, we do not update the - ! variable 'temperature' in prim2con at the interfaces - ! We will instead use an average temperature as an initial - ! guess. - xtemp(1) = 0.5d0*(temperature(i,j,k) + & - temperature(i-xoffset,j-yoffset,k-zoffset)) - if (y_e_minus(i,j,k) .le. 0.0d0 .or. y_e_plus(i,j,k) .le. 0.0d0) then !$OMP CRITICAL write(NaN_WarnLine,'(a100,7g15.6)') '(y_e_minus,y_e_plus,x,y,z,rho)', y_e(i,j,k), y_e_minus(i,j,k), y_e_plus(i,j,k), x(i,j,k),y(i,j,k),z(i,j,k),rho(i,j,k) @@ -197,14 +189,7 @@ Bconsxminus(i,j,k),Bconsyminus(i,j,k),Bconszminus(i,j,k), rhominus(i,j,k), & velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),& epsminus(i,j,k),pressminus(i,j,k),Bvecxminus(i,j,k), & - Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k), xtemp, y_e_minus(i,j,k)) - ! we do not have plus/minus vars for temperature since - ! eps is reconstructed. Hence, we do not update the - ! variable 'temperature' in prim2con at the interfaces - ! We will instead use an average temperature as an initial - ! guess. - xtemp(1) = 0.5d0*(temperature(i,j,k) + & - temperature(i+xoffset,j+yoffset,k+zoffset)) + Bvecyminus(i,j,k), Bveczminus(i,j,k), w_lorentzminus(i, j, k), tempminus(i,j,k), y_e_minus(i,j,k)) call prim2conM_hot(GRHydro_eos_handle, GRHydro_reflevel,& i,j,k,x(i,j,k),y(i,j,k),z(i,j,k), g11r,g12r,g13r,g22r,g23r,g33r, & @@ -214,7 +199,7 @@ rhoplus(i,j,k),velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),& Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k), & - w_lorentzplus(i,j,k), xtemp, y_e_plus(i,j,k)) + w_lorentzplus(i,j,k), tempplus(i,j,k), y_e_plus(i,j,k)) end do end do @@ -265,7 +250,7 @@ integer :: n, keytemp, anyerr, keyerr(1) ! real*8 :: xpress(1),xeps(1),xtemp(1),xye(1) real*8 :: temp0(1) - n = 1; keytemp = 0; anyerr = 0; keyerr(1) = 0 + n = 1; keytemp = reconstruct_temper; anyerr = 0; keyerr(1) = 0 !xpress = 0.0d0; xeps = 0.0d0; xtemp = 0.0d0; xye = 0.0d0 ! end EOS Omni vars @@ -299,11 +284,11 @@ ! in this case (coarse grid error that is hopefully restricted ! away), we use the average temperature between cells and call ! the EOS with keytemp=1 - keytemp=1 + keytemp=1-reconstruct_temper temp = temp0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) - keytemp=0 + keytemp=reconstruct_temper if(anyerr.ne.0) then !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot: lev 2") @@ -332,11 +317,11 @@ write(warnline,"(A10,i5)") "reflevel: ", GRHydro_reflevel call CCTK_WARN(1,warnline) !$OMP END CRITICAL - keytemp=1 + keytemp=1-reconstruct_temper temp = temp0 call EOS_Omni_press(handle,keytemp,GRHydro_eos_rf_prec,n,& drho,deps,temp,ye,dpress,keyerr,anyerr) - keytemp=0 + keytemp=reconstruct_temper if(anyerr.ne.0) then !$OMP CRITICAL call CCTK_WARN(1,"EOS error in prim2con_hot")

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 633)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_Prim2Con.c, GRHydro_Prim2ConM.F90 Log: GRHydro: make routine to compute speed of sound callable from F90 MHD p2c the code is long and complicated enough that it should better not be copied and pasted. File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_Prim2Con.c Delta lines: +77 -62 =================================================================== --- trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:49:47 UTC (rev 632) +++ trunk/src/GRHydro_Prim2Con.c 2014-04-15 19:49:50 UTC (rev 633) @@ -25,17 +25,21 @@ const double rho, const double vx, const double vy, const double vz, const double eps, const double press); +static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS); - void Primitive2ConservativeC(CCTK_ARGUMENTS); CCTK_FCALL void CCTK_FNAME(Primitive2ConservativeCforF)(cGH ** p_cctkGH) { Primitive2ConservativeC(*p_cctkGH); } +CCTK_FCALL void CCTK_FNAME(GRHydro_SpeedOfSound)(cGH ** p_cctkGH) { + GRHydro_SpeedOfSound(*p_cctkGH); +} + void Primitive2ConservativeC(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; @@ -61,12 +65,73 @@ g33 = gzz; } + GRHydro_SpeedOfSound(CCTK_PASS_CTOC); + // if padding is used the the "extra" vector elements could contain junk // which we do not know how to handle assert(cctk_lsh[0] == cctk_ash[0]); assert(cctk_lsh[1] == cctk_ash[1]); assert(cctk_lsh[2] == cctk_ash[2]); +#pragma omp parallel for + for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil+1; k++) + for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil+1; j++) +#pragma ivdep // force compiler to vectorize the goddamn loop + for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil+1; i++) { + + const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); + + const int idxl = CCTK_GFINDEX3D(cctkGH,i-*xoffset,j-*yoffset,k-*zoffset); + const int idxr = CCTK_GFINDEX3D(cctkGH,i+*xoffset,j+*yoffset,k+*zoffset); + + const double g11l = 0.5 * (g11[idx] + g11[idxl]); + const double g12l = 0.5 * (g12[idx] + g12[idxl]); + const double g13l = 0.5 * (g13[idx] + g13[idxl]); + const double g22l = 0.5 * (g22[idx] + g22[idxl]); + const double g23l = 0.5 * (g23[idx] + g23[idxl]); + const double g33l = 0.5 * (g33[idx] + g33[idxl]); + + const double g11r = 0.5 * (g11[idx] + g11[idxr]); + const double g12r = 0.5 * (g12[idx] + g12[idxr]); + const double g13r = 0.5 * (g13[idx] + g13[idxr]); + const double g22r = 0.5 * (g22[idx] + g22[idxr]); + const double g23r = 0.5 * (g23[idx] + g23[idxr]); + const double g33r = 0.5 * (g33[idx] + g33[idxr]); + + const double savg_detl = + sqrt(SpatialDeterminantC(g11l,g12l,g13l,g22l,g23l,g33l)); + const double savg_detr = + sqrt(SpatialDeterminantC(g11r,g12r,g13r,g22r,g23r,g33r)); + + // minus call to p2c + prim2conC(&w_lorentzminus[idx], &densminus[idx], &sxminus[idx], + &syminus[idx], &szminus[idx], &tauminus[idx], + g11l,g12l,g13l,g22l,g23l,g33l, + savg_detl,rhominus[idx], velxminus[idx], velyminus[idx], + velzminus[idx], epsminus[idx], pressminus[idx]); + + + // plus call to p2c + prim2conC(&w_lorentzplus[idx], &densplus[idx], &sxplus[idx], + &syplus[idx], &szplus[idx], &tauplus[idx], + g11r,g12r,g13r,g22r,g23r,g33r, + savg_detr,rhoplus[idx], velxplus[idx], velyplus[idx], + velzplus[idx], epsplus[idx], pressplus[idx]); + } + +} // end function Conservative2PrimitiveC + +static inline void GRHydro_SpeedOfSound(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_PARAMETERS; + DECLARE_CCTK_ARGUMENTS; + + // if padding is used the the "extra" vector elements could contain junk + // which we do not know how to handle + assert(cctk_lsh[0] == cctk_ash[0]); + assert(cctk_lsh[1] == cctk_ash[1]); + assert(cctk_lsh[2] == cctk_ash[2]); + // EOS calls (now GF-wide) if(!*evolve_temper) { // n needs to be computed using ash since ash is used when computing the @@ -92,7 +157,7 @@ } else { if(reconstruct_temper) { int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; - int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1; + int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset); int *keyerr = malloc(sizeof(*keyerr)*n); int keytemp = 1; // ensure Y_e and temperature within bounds @@ -107,8 +172,8 @@ // call the EOS with slices #pragma omp parallel for - for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++) - for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) { + for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++) + for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) { int anyerr = 0; int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k); EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx, @@ -136,8 +201,8 @@ #pragma omp parallel for - for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++) - for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) { + for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++) + for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) { int anyerr = 0; int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k); EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx, @@ -166,7 +231,7 @@ } else { // ******************** EPS RECONSTRUCTION BRANCH ****************** int n = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; - int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1; + int nx = cctk_lsh[0] - (GRHydro_stencil - 1) - (GRHydro_stencil) + 1 + (transport_constraints && !*xoffset); int *keyerr = malloc(sizeof(*keyerr)*n); int keytemp = 0; @@ -183,8 +248,8 @@ // call the EOS with slices for minus states #pragma omp parallel for - for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++) - for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) { + for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++) + for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) { int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k); int anyerr = 0; EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx, @@ -233,8 +298,8 @@ // call the EOS with slices for plus states #pragma omp parallel for - for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1;k++) - for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1;j++) { + for(int k=GRHydro_stencil-1;k<cctk_lsh[2]-GRHydro_stencil+1 + (transport_constraints && !*zoffset);k++) + for(int j=GRHydro_stencil-1;j<cctk_lsh[1]-GRHydro_stencil+1 + (transport_constraints && !*yoffset);j++) { int i = CCTK_GFINDEX3D(cctkGH,GRHydro_stencil-1,j,k); int anyerr = 0; EOS_Omni_press_cs2(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,nx, @@ -278,64 +343,14 @@ } // omp critical pragma } // if keyerr } // loop ii - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Aborting!"); } // if (anyerr) } // big loop over plus states free(keyerr); } // end branch for no temp reconsturction } // end of evolve temper branch +} - -#pragma omp parallel for - for(int k = GRHydro_stencil-1; k < cctk_lsh[2]-GRHydro_stencil+1; k++) - for(int j = GRHydro_stencil-1; j < cctk_lsh[1]-GRHydro_stencil+1; j++) -#pragma ivdep // force compiler to vectorize the goddamn loop - for(int i = GRHydro_stencil-1; i < cctk_lsh[0]-GRHydro_stencil+1; i++) { - - const int idx = CCTK_GFINDEX3D(cctkGH,i,j,k); - - const int idxl = CCTK_GFINDEX3D(cctkGH,i-*xoffset,j-*yoffset,k-*zoffset); - const int idxr = CCTK_GFINDEX3D(cctkGH,i+*xoffset,j+*yoffset,k+*zoffset); - - const double g11l = 0.5 * (g11[idx] + g11[idxl]); - const double g12l = 0.5 * (g12[idx] + g12[idxl]); - const double g13l = 0.5 * (g13[idx] + g13[idxl]); - const double g22l = 0.5 * (g22[idx] + g22[idxl]); - const double g23l = 0.5 * (g23[idx] + g23[idxl]); - const double g33l = 0.5 * (g33[idx] + g33[idxl]); - - const double g11r = 0.5 * (g11[idx] + g11[idxr]); - const double g12r = 0.5 * (g12[idx] + g12[idxr]); - const double g13r = 0.5 * (g13[idx] + g13[idxr]); - const double g22r = 0.5 * (g22[idx] + g22[idxr]); - const double g23r = 0.5 * (g23[idx] + g23[idxr]); - const double g33r = 0.5 * (g33[idx] + g33[idxr]); - - const double savg_detl = - sqrt(SpatialDeterminantC(g11l,g12l,g13l,g22l,g23l,g33l)); - const double savg_detr = - sqrt(SpatialDeterminantC(g11r,g12r,g13r,g22r,g23r,g33r)); - - // minus call to p2c - prim2conC(&w_lorentzminus[idx], &densminus[idx], &sxminus[idx], - &syminus[idx], &szminus[idx], &tauminus[idx], - g11l,g12l,g13l,g22l,g23l,g33l, - savg_detl,rhominus[idx], velxminus[idx], velyminus[idx], - velzminus[idx], epsminus[idx], pressminus[idx]); - - - // plus call to p2c - prim2conC(&w_lorentzplus[idx], &densplus[idx], &sxplus[idx], - &syplus[idx], &szplus[idx], &tauplus[idx], - g11r,g12r,g13r,g22r,g23r,g33r, - savg_detr,rhoplus[idx], velxplus[idx], velyplus[idx], - velzplus[idx], epsplus[idx], pressplus[idx]); - } - -} // end function Conservative2PrimitiveC - static inline double SpatialDeterminantC(double gxx, double gxy, double gxz, double gyy, double gyz, double gzz) File [modified]: GRHydro_Prim2ConM.F90 Delta lines: +11 -1 =================================================================== --- trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:49:47 UTC (rev 632) +++ trunk/src/GRHydro_Prim2ConM.F90 2014-04-15 19:49:50 UTC (rev 633) @@ -85,6 +85,17 @@ pBprim = loc(Bvec) end if + ! currently only the C++ code uses the speed-of-sound grid function, however + ! in the MHD code, both C++ and F90 code share the same prim2con routine (this + ! one). GRHydro_SpeedOfSound is a misnomer since it computes both cs2 and + ! pressure so we call it first so that the actual prim2con below overwrites + ! the pressur again. This is certainly not the "right" way of doing this but + ! it lets us have a working version for now. + ! TODO: extend prim2con in GRHydro_Prim2Con.c to handle MHD. + if (use_cxx_code .ne. 0) then + call GRHydro_SpeedOfSound(cctkGH) + end if + ! constraint transport needs to be able to average fluxes in the directions ! other that flux_direction, which in turn need the primitives on interfaces @@ -211,7 +222,6 @@ !$OMP END PARALLEL DO endif - end subroutine primitive2conservativeM /*@@

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 632)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_HLLEM.F90 Log: GRHydro: remove unused variables File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_HLLEM.F90 Delta lines: +3 -9 =================================================================== --- trunk/src/GRHydro_HLLEM.F90 2014-04-15 19:49:45 UTC (rev 631) +++ trunk/src/GRHydro_HLLEM.F90 2014-04-15 19:49:47 UTC (rev 632) @@ -52,7 +52,7 @@ CCTK_REAL, dimension(3) :: vtp,vtm,blowp,blowm,Bveclowp,Bveclowm CCTK_REAL, dimension(3) :: vellowp,vellowm CCTK_REAL :: ab0p,ab0m,b2p,b2m,bdotvp,bdotvm - CCTK_REAL :: wp,wm,v2p,v2m,vA2m,vA2p + CCTK_REAL :: wp,wm,v2p,v2m CCTK_REAL :: pressstarp,pressstarm CCTK_REAL :: entropyconsp,entropyconsm,entropyp,entropym,entropyf,entropydiff,entropyfp,entropyfm @@ -133,7 +133,7 @@ !$OMP vtp,vtm,vellowp,vellowm,Bveclowp,Bveclowm,& !$OMP bdotvp,bdotvm,b2p,b2m,v2p,v2m,wp,wm,& !$OMP blowp,blowm,& - !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,vA2p,vA2m,pressstarp,pressstarm,& + !$OMP rhoenth_p,rhoenth_m,ab0p,ab0m,pressstarp,pressstarm,& !$OMP usendh,psidcdiff,psidcfp,psidcfm,charmin,charmax,chartop,charpm,& !$OMP charmin_dc,charmax_dc,charpm_dc,m,xtemp,& !$OMP entropyconsp,entropyconsm,entropyp,entropym,entropyf,entropydiff,entropyfp,entropyfm) @@ -271,9 +271,6 @@ ab0p = wp*bdotvp ab0m = wm*bdotvm - vA2p = b2p/(rhoenth_p+b2p) - vA2m = b2m/(rhoenth_m+b2m) - !!$ p^* = p+b^2/2 in Anton et al. pressstarp = prim_p(6)+0.5d0*b2p pressstarm = prim_m(6)+0.5d0*b2m @@ -555,7 +552,7 @@ CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh - CCTK_REAL :: b2p,b2m,vA2m,vA2p + CCTK_REAL :: b2p,b2m CCTK_INT :: type_bits, trivial @@ -707,9 +704,6 @@ b2m=DOTP2(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,mag_m(1),mag_m(2),mag_m(3))/prim_m(7)**2 + & (DOTP(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,prim_m(2),prim_m(3),prim_m(4),mag_m(1),mag_m(2),mag_m(3)))**2 - vA2p = b2p/(prim_p(1)*(1.0d0+prim_p(5))+prim_p(6)+b2p) - vA2m = b2m/(prim_m(1)*(1.0d0+prim_m(5))+prim_m(6)+b2m) - !!$ Calculate the jumps in the conserved variables qdiff = cons_m - cons_p

1 0

[svn:einsteintoolkit] GRHydro/trunk/src/ (Rev. 631)
by rhaas＠tapir.caltech.edu 15 Apr '14

15 Apr '14

User: rhaas Date: 2014/04/15 02:49 PM Modified: /trunk/src/ GRHydro_CalcUpdate.F90 Log: GRHydro: low level optimizations File Changes: Directory: /trunk/src/ ====================== File [modified]: GRHydro_CalcUpdate.F90 Delta lines: +11 -7 =================================================================== --- trunk/src/GRHydro_CalcUpdate.F90 2014-04-15 19:49:43 UTC (rev 630) +++ trunk/src/GRHydro_CalcUpdate.F90 2014-04-15 19:49:45 UTC (rev 631) @@ -39,7 +39,7 @@ DECLARE_CCTK_FUNCTIONS CCTK_INT :: i,j,k,itracer - CCTK_REAL :: idx, alp_l, alp_r, Bcons_l, Bcons_r, alp_tmp + CCTK_REAL :: idx, idy, idz, alp_l, alp_r, Bcons_l, Bcons_r, alp_tmp idx = 1.d0 / CCTK_DELTA_SPACE(flux_direction) @@ -300,22 +300,26 @@ flux_direction.eq.1) then ! HACK: x direction is last ! FIXME: I think one could wrap all of this into a single do loop and remove the ! Evec storage + ! idx differs from idx which was 1d0/CCTK_DELTA_SPACE(flux_direction) + idx = 1d0 / CCTK_DELTA_SPACE(1) + idy = 1d0 / CCTK_DELTA_SPACE(2) + idz = 1d0 / CCTK_DELTA_SPACE(3) !$OMP PARALLEL DO PRIVATE(i,j,k) do k = GRHydro_stencil + 1, cctk_lsh(3) - GRHydro_stencil do j = GRHydro_stencil + 1, cctk_lsh(2) - GRHydro_stencil do i = GRHydro_stencil + 1, cctk_lsh(1) - GRHydro_stencil Bconsrhs(i,j,k,1) = - 0.5d0 * ((Evec(i-1,j ,k-1,2)-Evec(i-1,j ,k ,2)) & - + (Evec(i ,j ,k-1,2)-Evec(i ,j ,k ,2))) / CCTK_DELTA_SPACE(3) & + + (Evec(i ,j ,k-1,2)-Evec(i ,j ,k ,2))) * idz & - 0.5d0 * ((Evec(i-1,j ,k ,3)-Evec(i-1,j-1,k ,3)) & - + (Evec(i ,j ,k ,3)-Evec(i ,j-1,k ,3))) / CCTK_DELTA_SPACE(2) + + (Evec(i ,j ,k ,3)-Evec(i ,j-1,k ,3))) * idy Bconsrhs(i,j,k,2) = - 0.5d0 * ((Evec(i-1,j-1,k ,3)-Evec(i ,j-1,k ,3)) & - + (Evec(i-1,j ,k ,3)-Evec(i ,j ,k ,3))) / CCTK_DELTA_SPACE(1) & + + (Evec(i-1,j ,k ,3)-Evec(i ,j ,k ,3))) * idx & - 0.5d0 * ((Evec(i ,j-1,k ,1)-Evec(i ,j-1,k-1,1)) & - + (Evec(i ,j ,k ,1)-Evec(i ,j ,k-1,1))) / CCTK_DELTA_SPACE(3) + + (Evec(i ,j ,k ,1)-Evec(i ,j ,k-1,1))) * idz Bconsrhs(i,j,k,3) = - 0.5d0 * ((Evec(i ,j-1,k-1,1)-Evec(i ,j ,k-1,1)) & - + (Evec(i ,j-1,k ,1)-Evec(i ,j ,k ,1))) / CCTK_DELTA_SPACE(2) & + + (Evec(i ,j-1,k ,1)-Evec(i ,j ,k ,1))) * idy & - 0.5d0 * ((Evec(i ,j ,k-1,2)-Evec(i-1,j ,k-1,2)) & - + (Evec(i ,j ,k ,2)-Evec(i-1,j ,k ,2))) / CCTK_DELTA_SPACE(1) + + (Evec(i ,j ,k ,2)-Evec(i-1,j ,k ,2))) * idx enddo enddo enddo

1 0

Jump to page:

2026

2025

2024

2023

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

Commits