From 2b6919940238fa6fa03f75237e8d30342f42e799 Mon Sep 17 00:00:00 2001 From: rhaas Date: Sat, 6 Jul 2013 18:11:18 +0000 Subject: GRHydro: re-implement prim2con routine in C * new experimental prim2con routine that is about twice as fast as the old one. * Right now tested only for simple EOS. * Right now handles only prim2con call after reconstruction and must be enabled by changing comments in GRHydro_Reconstruct.F90 * remove duplicate code in all reconstruction *_drv.F90 routines and just have the initialization of plus and minus variables in the main GRHydro_Reconstruct.F90 routine From: Christian Ott git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@556 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- src/GRHydro_Con2Prim.F90 | 2 +- src/GRHydro_ENOReconstruct_drv.F90 | 64 +------ src/GRHydro_MP5Reconstruct_drv.F90 | 59 +------ src/GRHydro_PPMMReconstruct_drv.F90 | 56 +----- src/GRHydro_PPMReconstruct_drv.F90 | 42 ----- src/GRHydro_Prim2Con.c | 334 ++++++++++++++++++++++++++++++++++++ src/GRHydro_Reconstruct.F90 | 81 +++++++-- src/GRHydro_ReconstructPoly.F90 | 122 ++++++------- src/GRHydro_TVDReconstruct_drv.F90 | 61 +------ src/GRHydro_WENOReconstruct_drv.F90 | 61 +------ src/make.code.defn | 4 +- 11 files changed, 494 insertions(+), 392 deletions(-) create mode 100644 src/GRHydro_Prim2Con.c diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 176ba77..64785b5 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -684,7 +684,7 @@ subroutine Con2Prim_pt(cctk_iteration,ii,jj,kk,& IF_BELOW_ATMO(rho, GRHydro_rho_min, GRHydro_atmo_tolerance, r) then SET_ATMO_MIN(rho, GRHydro_rho_min, r) !GRHydro_rho_min udens = rho - dens = isdetg * rho + dens = sdetg * rho ! epsilon = (sqrt( (utau + pnew + udens)**2) - pnew - udens)/udens epsilon = epsmin ! w_lorentz=1, so the expression for utau reduces to: diff --git a/src/GRHydro_ENOReconstruct_drv.F90 b/src/GRHydro_ENOReconstruct_drv.F90 index 17a2494..0dfdde9 100644 --- a/src/GRHydro_ENOReconstruct_drv.F90 +++ b/src/GRHydro_ENOReconstruct_drv.F90 @@ -86,6 +86,8 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) call CCTK_WARN(0, "Allocation problems with trivial_rp") end if + + !!$ The dum variables are used as dummies if MHD on but divergence cleaning off allocate(dum(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& dump(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),& @@ -115,64 +117,12 @@ subroutine GRHydro_ENOReconstruct_drv(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) -!!$ Initialize variables that store reconstructed quantities: - - !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) + !initialize trivial_rp to false + !$OMP PARALLEL DO + do k=1,nz + do j=1,ny + do i=1,nx trivial_rp(i,j,k) = .false. - rhoplus(i,j,k) = 0.0d0 - rhominus(i,j,k)= 0.0d0 - epsplus(i,j,k) = 0.0d0 - epsminus(i,j,k) = 0.0d0 - velxplus(i,j,k) = 0.0d0 - velxminus(i,j,k) = 0.0d0 - velyplus(i,j,k) = 0.0d0 - velyminus(i,j,k) = 0.0d0 - velzplus(i,j,k) = 0.0d0 - velzminus(i,j,k) = 0.0d0 - - if(evolve_mhd.ne.0) then - Bvecxplus(i,j,k) = 0.0d0 - Bvecxminus(i,j,k) = 0.0d0 - Bvecyplus(i,j,k) = 0.0d0 - Bvecyminus(i,j,k) = 0.0d0 - Bveczplus(i,j,k) = 0.0d0 - Bveczminus(i,j,k) = 0.0d0 - if(clean_divergence.ne.0) then - psidcplus(i,j,k) = 0.0d0 - psidcminus(i,j,k) = 0.0d0 - endif - endif - - if (evolve_entropy .ne. 0) then - entropyplus(i,j,k) = 0.0d0 - entropyminus(i,j,k) = 0.0d0 - endif - - if (evolve_tracer .ne. 0) then - tracerplus(i,j,k,:) = 0.0d0 - tracerminus(i,j,k,:) = 0.0d0 - endif - - if (evolve_Y_e .ne. 0) then - ! set this to the cell center values - ! to make sure we have good Y_e even in - ! the boundary region (for full GF EOS calls) - Y_e_plus(i,j,k) = Y_e(i,j,k) - Y_e_minus(i,j,k) = Y_e(i,j,k) - endif - - if (evolve_temper .ne. 0) then - ! initialize with cell-center values - ! to make sure we have safe initial - ! guesses at interfaces even if - ! temperature is not reconstructed - tempplus(i,j,k) = temperature(i,j,k) - tempminus(i,j,k) = temperature(i,j,k) - endif - enddo enddo enddo diff --git a/src/GRHydro_MP5Reconstruct_drv.F90 b/src/GRHydro_MP5Reconstruct_drv.F90 index fc839c0..bb222dd 100644 --- a/src/GRHydro_MP5Reconstruct_drv.F90 +++ b/src/GRHydro_MP5Reconstruct_drv.F90 @@ -133,65 +133,18 @@ subroutine GRHydro_MP5Reconstruct_drv(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) -!!$ Initialize variables that store reconstructed quantities: - - !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) + !initialize trivial_rp to false + !$OMP PARALLEL DO + do k=1,nz + do j=1,ny + do i=1,nx trivial_rp(i,j,k) = .false. - rhoplus(i,j,k) = 0.0d0 - rhominus(i,j,k)= 0.0d0 - epsplus(i,j,k) = 0.0d0 - epsminus(i,j,k) = 0.0d0 - velxplus(i,j,k) = 0.0d0 - velxminus(i,j,k) = 0.0d0 - velyplus(i,j,k) = 0.0d0 - velyminus(i,j,k) = 0.0d0 - velzplus(i,j,k) = 0.0d0 - velzminus(i,j,k) = 0.0d0 - - if(evolve_mhd.ne.0) then - Bvecxplus(i,j,k) = 0.0d0 - Bvecxminus(i,j,k) = 0.0d0 - Bvecyplus(i,j,k) = 0.0d0 - Bvecyminus(i,j,k) = 0.0d0 - Bveczplus(i,j,k) = 0.0d0 - Bveczminus(i,j,k) = 0.0d0 - if(clean_divergence.ne.0) then - psidcplus(i,j,k) = 0.0d0 - psidcminus(i,j,k) = 0.0d0 - endif - endif - - if (evolve_entropy .ne. 0) then - entropyplus(i,j,k) = 0.0d0 - entropyminus(i,j,k) = 0.0d0 - endif - - if (evolve_tracer .ne. 0) then - tracerplus(i,j,k,:) = 0.0d0 - tracerminus(i,j,k,:) = 0.0d0 - endif - - if (evolve_Y_e .ne. 0) then - Y_e_plus(i,j,k) = 0.0d0 - Y_e_minus(i,j,k) = 0.0d0 - endif - - if (evolve_temper .ne. 0) then - ! set this to cell center value to have - ! good initial guesses at interfaces - ! in case we don't reconstruct temp - tempplus(i,j,k) = temperature(i,j,k) - tempminus(i,j,k) = temperature(i,j,k) - endif - enddo enddo enddo !$OMP END PARALLEL DO + !!$ MP5 starts: if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers diff --git a/src/GRHydro_PPMMReconstruct_drv.F90 b/src/GRHydro_PPMMReconstruct_drv.F90 index 370d764..f9ceaf7 100644 --- a/src/GRHydro_PPMMReconstruct_drv.F90 +++ b/src/GRHydro_PPMMReconstruct_drv.F90 @@ -135,60 +135,16 @@ subroutine GRHydro_PPMMReconstruct_drv(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) -!!$ Initialize variables that store reconstructed quantities: - - !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) + !initialize trivial_rp to false + !$OMP PARALLEL DO + do k=1,nz + do j=1,ny + do i=1,nx trivial_rp(i,j,k) = .false. - rhoplus(i,j,k) = 0.0d0 - rhominus(i,j,k)= 0.0d0 - epsplus(i,j,k) = 0.0d0 - epsminus(i,j,k) = 0.0d0 - velxplus(i,j,k) = 0.0d0 - velxminus(i,j,k) = 0.0d0 - velyplus(i,j,k) = 0.0d0 - velyminus(i,j,k) = 0.0d0 - velzplus(i,j,k) = 0.0d0 - velzminus(i,j,k) = 0.0d0 - - Bvecxplus(i,j,k) = 0.0d0 - Bvecxminus(i,j,k) = 0.0d0 - Bvecyplus(i,j,k) = 0.0d0 - Bvecyminus(i,j,k) = 0.0d0 - Bveczplus(i,j,k) = 0.0d0 - Bveczminus(i,j,k) = 0.0d0 - if(clean_divergence.ne.0) then - psidcplus(i,j,k) = 0.0d0 - psidcminus(i,j,k) = 0.0d0 - endif - - if (evolve_entropy .ne. 0) then - entropyplus(i,j,k) = 0.0d0 - entropyminus(i,j,k) = 0.0d0 - endif - - if (evolve_tracer .ne. 0) then - tracerplus(i,j,k,:) = 0.0d0 - tracerminus(i,j,k,:) = 0.0d0 - endif - - if (evolve_Y_e .ne. 0) then - Y_e_plus(i,j,k) = Y_e(i,j,k) - Y_e_minus(i,j,k) = Y_e(i,j,k) - endif - - if(evolve_temper .ne. 0) then - tempplus(i,j,k) = temperature(i,j,k) - tempminus(i,j,k) = temperature(i,j,k) - endif - - enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END PARALLEL DO !!$ PPM starts: if (flux_direction == 1) then diff --git a/src/GRHydro_PPMReconstruct_drv.F90 b/src/GRHydro_PPMReconstruct_drv.F90 index da4d1cf..f8959a2 100644 --- a/src/GRHydro_PPMReconstruct_drv.F90 +++ b/src/GRHydro_PPMReconstruct_drv.F90 @@ -128,48 +128,6 @@ subroutine GRHydro_PPMReconstruct_drv(CCTK_ARGUMENTS) apply_enhanced_ppm = .false. end if -!!$ Initialize variables that store reconstructed quantities: - - !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) - trivial_rp(i,j,k) = .false. - rhoplus(i,j,k) = 0.0d0 - rhominus(i,j,k)= 0.0d0 - epsplus(i,j,k) = 0.0d0 - epsminus(i,j,k) = 0.0d0 - velxplus(i,j,k) = 0.0d0 - velxminus(i,j,k) = 0.0d0 - velyplus(i,j,k) = 0.0d0 - velyminus(i,j,k) = 0.0d0 - velzplus(i,j,k) = 0.0d0 - velzminus(i,j,k) = 0.0d0 - - if (evolve_tracer .ne. 0) then - tracerplus(i,j,k,:) = 0.0d0 - tracerminus(i,j,k,:) = 0.0d0 - endif - - if (evolve_Y_e .ne. 0) then - ! set this to the cell center values - ! to make sure we have good Y_e even in - ! the boundary region (for full GF EOS calls) - Y_e_plus(i,j,k) = Y_e(i,j,k) - Y_e_minus(i,j,k) = Y_e(i,j,k) - endif - - if (evolve_temper .ne. 0) then - ! set this to temperature (piecewise constant), because we won't - ! set it if reconstruct_temper .eq. 0 - tempplus(i,j,k) = temperature(i,j,k) - tempminus(i,j,k) = temperature(i,j,k) - endif - - enddo - enddo - enddo - !$OMP END PARALLEL DO !!$ PPM starts: if (flux_direction == 1) then diff --git a/src/GRHydro_Prim2Con.c b/src/GRHydro_Prim2Con.c new file mode 100644 index 0000000..36578b0 --- /dev/null +++ b/src/GRHydro_Prim2Con.c @@ -0,0 +1,334 @@ +#include +#include +#include +#include + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" +#define MIN(a,b) (((a)<(b))?(a):(b)) +#define MAX(a,b) (((a)>(b))?(a):(b)) + + +// some prototypes +CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); +static inline double SpatialDeterminantC(double gxx, double gxy, + double gxz, double gyy, + double gyz, double gzz); + +static inline __attribute__((always_inline)) void prim2conC(double *w, double *dens, double *sx, + double *sy, double *sz, double *tau, + const double gxx, const double gxy, + const double gxz, const double gyy, const double gyz, + const double gzz, const double sdet, + const double rho, const double vx, + const double vy, const double vz, + const double eps, const double press); + + + + +void Primitive2ConservativeC(CCTK_ARGUMENTS); + +CCTK_FCALL void CCTK_FNAME(Primitive2ConservativeCforF)(cGH ** p_cctkGH) { + Primitive2ConservativeC(*p_cctkGH); +} + + +void Primitive2ConservativeC(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + // save memory when multipatch is not used + CCTK_REAL * restrict g11, * restrict g12, * restrict g13; + CCTK_REAL * restrict g22, * restrict g23, * restrict g33; + + if(GRHydro_UseGeneralCoordinates(cctkGH)) { + g11 = gaa; + g12 = gab; + g13 = gac; + g22 = gbb; + g23 = gbc; + g33 = gcc; + } else { + g11 = gxx; + g12 = gxy; + g13 = gxz; + g22 = gyy; + g23 = gyz; + g33 = gzz; + } + + // 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) { + int n = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; + int keyerr[n]; int anyerr = 0; + int keytemp = 0; + + // don't need special error handling for analytic EOS + EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, + rhominus,epsminus,NULL,NULL,pressminus,keyerr,&anyerr); + + EOS_Omni_pressOMP(*GRHydro_eos_handle,keytemp,GRHydro_eos_rf_prec,n, + rhoplus,epsplus,NULL,NULL,pressplus,keyerr,&anyerr); + + } else { + if(reconstruct_temper) { + int n = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]; + // Linux usually offers ~8MB of stack which is good for ~1M points. + // We should consider malloc() if we ever smash the stack. + int keyerr[n]; int anyerr = 0; + int keytemp = 1; + + // ensure Y_e and temperature within bounds +#pragma omp parallel for + for(int i=0;i gzz vup => vel end if + + ! Initialize plus and minus states + !$OMP PARALLEL DO + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + ! must initialize rho and eps plus minus + ! to cell average in order to have sane values + ! in the boundary zones for EOS call + rhoplus(i,j,k) = rho(i,j,k) + rhominus(i,j,k)= rho(i,j,k) + epsplus(i,j,k) = eps(i,j,k) + epsminus(i,j,k) = eps(i,j,k) + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 + + if(evolve_mhd.ne.0) then + Bvecxplus(i,j,k) = 0.0d0 + Bvecxminus(i,j,k) = 0.0d0 + Bvecyplus(i,j,k) = 0.0d0 + Bvecyminus(i,j,k) = 0.0d0 + Bveczplus(i,j,k) = 0.0d0 + Bveczminus(i,j,k) = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus(i,j,k) = 0.0d0 + psidcminus(i,j,k) = 0.0d0 + endif + endif + + if (evolve_entropy .ne. 0) then + entropyplus(i,j,k) = 0.0d0 + entropyminus(i,j,k) = 0.0d0 + endif + + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + ! set this to the cell center values + ! to make sure we have good Y_e even in + ! the boundary region (for full GF EOS calls) + Y_e_plus(i,j,k) = Y_e(i,j,k) + Y_e_minus(i,j,k) = Y_e(i,j,k) + endif + + if(evolve_temper .ne. 0) then + ! set this to cell center value to have + ! good initial guesses at interfaces + ! in case we don't reconstruct temp + tempplus(i,j,k) = temperature(i,j,k) + tempminus(i,j,k) = temperature(i,j,k) + endif + enddo + enddo + enddo + !$OMP END PARALLEL DO if (CCTK_EQUALS(recon_method,"tvd")) then ! this handles MHD and non-MHD @@ -225,14 +271,13 @@ subroutine Reconstruction(CCTK_ARGUMENTS) enddo !$OMP END PARALLEL DO - - if (CCTK_EQUALS(recon_vars,"primitive").or.& CCTK_EQUALS(recon_method,"ppm")) then if(evolve_mhd.ne.0) then call primitive2conservativeM(CCTK_PASS_FTOF) else call primitive2conservative(CCTK_PASS_FTOF) +! call Primitive2ConservativeCforF(cctkGH) endif else if (CCTK_EQUALS(recon_vars,"conservative")) then if(evolve_mhd.ne.0) then diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90 index 1e3a5b8..826a49e 100644 --- a/src/GRHydro_ReconstructPoly.F90 +++ b/src/GRHydro_ReconstructPoly.F90 @@ -123,69 +123,75 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) - !$OMP PARALLEL - !$OMP WORKSHARE - trivial_rp = .false. - !$OMP END WORKSHARE NOWAIT - !!$ Currently only option is reconstruction on primitive variables. !!$ Should change this. - !$OMP WORKSHARE - psi4 = 1.0d0 - !$OMP END WORKSHARE NOWAIT - - !$OMP WORKSHARE - lbetax = betax - lbetay = betay - lbetaz = betaz - !$OMP END WORKSHARE NOWAIT - -!!$ Initialize variables that store reconstructed quantities - - !$OMP WORKSHARE - rhoplus = 0.0d0 - rhominus = 0.0d0 - epsplus = 0.0d0 - epsminus = 0.0d0 - velxplus = 0.0d0 - velxminus = 0.0d0 - velyplus = 0.0d0 - velyminus = 0.0d0 - velzplus = 0.0d0 - velzminus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - if(evolve_mhd.ne.0) then - !$OMP WORKSHARE - Bvecxplus = 0.0d0 - Bvecxminus = 0.0d0 - Bvecyplus = 0.0d0 - Bvecyminus = 0.0d0 - Bveczplus = 0.0d0 - Bveczminus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - if(clean_divergence.ne.0) then - !$OMP WORKSHARE - psidcplus = 0.0d0 - psidcminus=0.0d0 - !$OMP END WORKSHARE NOWAIT - endif - endif + ! Initialize plus and minus states + !$OMP PARALLEL DO + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + psi4(i,j,k) = 1.0d0 + lbetax(i,j,k) = betax(i,j,k) + lbetay(i,j,k) = betay(i,j,k) + lbetaz(i,j,k) = betaz(i,j,k) + trivial_rp(i,j,k) = .false. + ! must initialize rho and eps plus minus + ! to cell average in order to have sane values + ! in the boundary zones for EOS call + rhoplus(i,j,k) = rho(i,j,k) + rhominus(i,j,k)= rho(i,j,k) + epsplus(i,j,k) = eps(i,j,k) + epsminus(i,j,k) = eps(i,j,k) + velxplus(i,j,k) = 0.0d0 + velxminus(i,j,k) = 0.0d0 + velyplus(i,j,k) = 0.0d0 + velyminus(i,j,k) = 0.0d0 + velzplus(i,j,k) = 0.0d0 + velzminus(i,j,k) = 0.0d0 + + if(evolve_mhd.ne.0) then + Bvecxplus(i,j,k) = 0.0d0 + Bvecxminus(i,j,k) = 0.0d0 + Bvecyplus(i,j,k) = 0.0d0 + Bvecyminus(i,j,k) = 0.0d0 + Bveczplus(i,j,k) = 0.0d0 + Bveczminus(i,j,k) = 0.0d0 + if(clean_divergence.ne.0) then + psidcplus(i,j,k) = 0.0d0 + psidcminus(i,j,k) = 0.0d0 + endif + endif + if (evolve_entropy .ne. 0) then + entropyplus(i,j,k) = 0.0d0 + entropyminus(i,j,k) = 0.0d0 + endif + + if (evolve_tracer .ne. 0) then + tracerplus(i,j,k,:) = 0.0d0 + tracerminus(i,j,k,:) = 0.0d0 + endif + + if (evolve_Y_e .ne. 0) then + ! set this to the cell center values + ! to make sure we have good Y_e even in + ! the boundary region (for full GF EOS calls) + Y_e_plus(i,j,k) = Y_e(i,j,k) + Y_e_minus(i,j,k) = Y_e(i,j,k) + endif + if(evolve_temper .ne. 0) then + ! set this to cell center value to have + ! good initial guesses at interfaces + ! in case we don't reconstruct temp + tempplus(i,j,k) = temperature(i,j,k) + tempminus(i,j,k) = temperature(i,j,k) + endif + enddo + enddo + enddo + !$OMP END PARALLEL DO - if (evolve_tracer .ne. 0) then - !$OMP WORKSHARE - tracerplus = 0.0d0 - tracerminus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - endif - if (evolve_Y_e .ne. 0) then - !$OMP WORKSHARE - Y_e_plus = 0.0d0 - Y_e_minus = 0.0d0 - !$OMP END WORKSHARE NOWAIT - endif - !$OMP END PARALLEL if (CCTK_EQUALS(recon_method,"tvd")) then diff --git a/src/GRHydro_TVDReconstruct_drv.F90 b/src/GRHydro_TVDReconstruct_drv.F90 index eb633b2..6925132 100644 --- a/src/GRHydro_TVDReconstruct_drv.F90 +++ b/src/GRHydro_TVDReconstruct_drv.F90 @@ -147,63 +147,12 @@ subroutine GRHydro_TVDReconstruct_drv(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) -!!$ Initialize variables that store reconstructed quantities: - - !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) + !initialize trivial_rp to false + !$OMP PARALLEL DO + do k=1,nz + do j=1,ny + do i=1,nx trivial_rp(i,j,k) = .false. - rhoplus(i,j,k) = 0.0d0 - rhominus(i,j,k)= 0.0d0 - epsplus(i,j,k) = 0.0d0 - epsminus(i,j,k) = 0.0d0 - velxplus(i,j,k) = 0.0d0 - velxminus(i,j,k) = 0.0d0 - velyplus(i,j,k) = 0.0d0 - velyminus(i,j,k) = 0.0d0 - velzplus(i,j,k) = 0.0d0 - velzminus(i,j,k) = 0.0d0 - - if(evolve_mhd.ne.0) then - Bvecxplus(i,j,k) = 0.0d0 - Bvecxminus(i,j,k) = 0.0d0 - Bvecyplus(i,j,k) = 0.0d0 - Bvecyminus(i,j,k) = 0.0d0 - Bveczplus(i,j,k) = 0.0d0 - Bveczminus(i,j,k) = 0.0d0 - if(clean_divergence.ne.0) then - psidcplus(i,j,k) = 0.0d0 - psidcminus(i,j,k) = 0.0d0 - endif - endif - - if (evolve_entropy .ne. 0) then - entropyplus(i,j,k) = 0.0d0 - entropyminus(i,j,k) = 0.0d0 - endif - - if (evolve_tracer .ne. 0) then - tracerplus(i,j,k,:) = 0.0d0 - tracerminus(i,j,k,:) = 0.0d0 - endif - - if (evolve_Y_e .ne. 0) then - ! set this to the cell center values - ! to make sure we have good Y_e even in - ! the boundary region (for full GF EOS calls) - Y_e_plus(i,j,k) = Y_e(i,j,k) - Y_e_minus(i,j,k) = Y_e(i,j,k) - endif - - if(evolve_temper .ne. 0) then - ! set this to cell center value to have - ! good initial guesses at interfaces - ! in case we don't reconstruct temp - tempplus(i,j,k) = temperature(i,j,k) - tempminus(i,j,k) = temperature(i,j,k) - endif - enddo enddo enddo diff --git a/src/GRHydro_WENOReconstruct_drv.F90 b/src/GRHydro_WENOReconstruct_drv.F90 index ae3bd01..396c54d 100644 --- a/src/GRHydro_WENOReconstruct_drv.F90 +++ b/src/GRHydro_WENOReconstruct_drv.F90 @@ -192,63 +192,12 @@ subroutine GRHydro_WENOReconstruct_drv(CCTK_ARGUMENTS) ny = cctk_lsh(2) nz = cctk_lsh(3) -!!$ Initialize variables that store reconstructed quantities: - - !$OMP PARALLEL DO PRIVATE(i,j,k) - do k=1,cctk_lsh(3) - do j=1,cctk_lsh(2) - do i=1,cctk_lsh(1) + !initialize trivial_rp to false + !$OMP PARALLEL DO + do k=1,nz + do j=1,ny + do i=1,nx trivial_rp(i,j,k) = .false. - rhoplus(i,j,k) = 0.0d0 - rhominus(i,j,k)= 0.0d0 - epsplus(i,j,k) = 0.0d0 - epsminus(i,j,k) = 0.0d0 - velxplus(i,j,k) = 0.0d0 - velxminus(i,j,k) = 0.0d0 - velyplus(i,j,k) = 0.0d0 - velyminus(i,j,k) = 0.0d0 - velzplus(i,j,k) = 0.0d0 - velzminus(i,j,k) = 0.0d0 - - if(evolve_mhd.ne.0) then - Bvecxplus(i,j,k) = 0.0d0 - Bvecxminus(i,j,k) = 0.0d0 - Bvecyplus(i,j,k) = 0.0d0 - Bvecyminus(i,j,k) = 0.0d0 - Bveczplus(i,j,k) = 0.0d0 - Bveczminus(i,j,k) = 0.0d0 - if(clean_divergence.ne.0) then - psidcplus(i,j,k) = 0.0d0 - psidcminus(i,j,k) = 0.0d0 - endif - endif - - if (evolve_entropy .ne. 0) then - entropyplus(i,j,k) = 0.0d0 - entropyminus(i,j,k) = 0.0d0 - endif - - if (evolve_tracer .ne. 0) then - tracerplus(i,j,k,:) = 0.0d0 - tracerminus(i,j,k,:) = 0.0d0 - endif - - if (evolve_Y_e .ne. 0) then - ! set this to the cell center values - ! to make sure we have good Y_e even in - ! the boundary region (for full GF EOS calls) - Y_e_plus(i,j,k) = Y_e(i,j,k) - Y_e_minus(i,j,k) = Y_e(i,j,k) - endif - - if (evolve_temper .ne. 0) then - ! set this to cell center value to have - ! good initial guesses at interfaces - ! in case we don't reconstruct temp - tempplus(i,j,k) = temperature(i,j,k) - tempminus(i,j,k) = temperature(i,j,k) - endif - enddo enddo enddo diff --git a/src/make.code.defn b/src/make.code.defn index b55a089..1fce82b 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -81,7 +81,9 @@ SRCS = Utils.F90 \ GRHydro_WENOScalars.F90 \ GRHydro_MP5Reconstruct.F90 \ Con2Prim_fortran_interfaces.F90 \ - Con2PrimM_fortran_interfaces.F90 + Con2PrimM_fortran_interfaces.F90 \ + GRHydro_Prim2Con.c + # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3