From ac63a7a15c061746ffd49d844f034165aa098fa8 Mon Sep 17 00:00:00 2001 From: rhaas Date: Tue, 15 Apr 2014 19:49:34 +0000 Subject: GRHydro: implement sources in C++. Combined result of commits by Christian Reisswig and Roland Haas. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@627 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45 --- param.ccl | 1 - schedule.ccl | 33 ++- src/GRHydro_Source.cc | 709 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 8 +- 4 files changed, 734 insertions(+), 17 deletions(-) create mode 100644 src/GRHydro_Source.cc diff --git a/param.ccl b/param.ccl index 3d08545..9a58eca 100644 --- a/param.ccl +++ b/param.ccl @@ -810,4 +810,3 @@ BOOLEAN constrain_to_1D "Set fluid velocities to zero for non-radial motion" BOOLEAN use_cxx_code "Use experimental C++ code?" STEERABLE = RECOVER { } no - diff --git a/schedule.ccl b/schedule.ccl index 5e08cfb..b7f2d58 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -732,22 +732,29 @@ if (number_of_particles) ### over the different directions ### #################################################################### -if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) -{ - schedule SourceTermsM IN GRHydroRHS BEFORE FluxTerms - { - LANG: Fortran - } "Source term calculation - MHD version" -} else if(CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { - schedule SourceTermsAM IN GRHydroRHS BEFORE FluxTerms - { - LANG: Fortran - } "Source term calculation - Vector Potential MHD version" -} else { +if(use_cxx_code) { schedule SourceTerms IN GRHydroRHS BEFORE FluxTerms { - LANG: Fortran + LANG: C } "Source term calculation" +} else { + if(CCTK_Equals(Bvec_evolution_method,"GRHydro")) + { + schedule SourceTermsM IN GRHydroRHS BEFORE FluxTerms + { + LANG: Fortran + } "Source term calculation - MHD version" + } else if(CCTK_Equals(Bvec_evolution_method,"GRHydro_Avec")) { + schedule SourceTermsAM IN GRHydroRHS BEFORE FluxTerms + { + LANG: Fortran + } "Source term calculation - Vector Potential MHD version" + } else { + schedule SourceTerms IN GRHydroRHS BEFORE FluxTerms + { + LANG: Fortran + } "Source term calculation" + } } ################################################################# ### Initial setup for the loop over the different directions. ### diff --git a/src/GRHydro_Source.cc b/src/GRHydro_Source.cc new file mode 100644 index 0000000..1858e45 --- /dev/null +++ b/src/GRHydro_Source.cc @@ -0,0 +1,709 @@ + /*@@ + @file GRHydro_Source.cc + @date Nov 29, 2013 + @author Christian Reisswig, Joshua Faber, Scott Noble, Bruno Mundim, Ian Hawke + @desc + The geometric source terms for the matter evolution + @enddesc + @@*/ + + + + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" +#include "cctk_Functions.h" + +#include +#include + +#define velx vup +#define vely (&vup[N]) +#define velz (&vup[2*N]) +#define Bvecx Bprim +#define Bvecy (&Bprim[N]) +#define Bvecz (&Bprim[2*N]) +#define Avecx Avec +#define Avecy (&Avec[N]) +#define Avecz (&Avec[2*N]) +#define Avecrhsx Avecrhs +#define Avecrhsy (&Avecrhs[N]) +#define Avecrhsz (&Avecrhs[2*N]) +#define Bconsx Bcons +#define Bconsy (&Bcons[N]) +#define Bconsz (&Bcons[2*N]) +#define Bconsrhsx Bconsrhs +#define Bconsrhsy (&Bconsrhs[N]) +#define Bconsrhsz (&Bconsrhs[2*N]) + + +extern "C" CCTK_INT GRHydro_UseGeneralCoordinates(const cGH * cctkGH); + +template static inline T SQR (T const & x) { return x*x; } + + +struct alldiff2 { + static void apply(const cGH* const restrict cctkGH, + CCTK_REAL dvars[][3], + const CCTK_REAL* const restrict * const restrict vars, + const int i, const int j, const int k, + const CCTK_REAL* restrict const ih, + const int nvars) + { + for (int n=0; n < nvars; ++n) { + dvars[n][0] = (vars[n][CCTK_GFINDEX3D(cctkGH, i+1,j,k)] - vars[n][CCTK_GFINDEX3D(cctkGH, i-1,j,k)]) * 0.5 * ih[0]; + dvars[n][1] = (vars[n][CCTK_GFINDEX3D(cctkGH, i,j+1,k)] - vars[n][CCTK_GFINDEX3D(cctkGH, i,j-1,k)]) * 0.5 * ih[1]; + dvars[n][2] = (vars[n][CCTK_GFINDEX3D(cctkGH, i,j,k+1)] - vars[n][CCTK_GFINDEX3D(cctkGH, i,j,k-1)]) * 0.5 * ih[2]; + } + } +}; + +struct alldiff4 { + static void apply(const cGH* const restrict cctkGH, + CCTK_REAL dvars[][3], + const CCTK_REAL* const restrict * const restrict vars, + const int i, const int j, const int k, + const CCTK_REAL* restrict const ih, + const int nvars) + { + for (int n=0; n < nvars; ++n) { + dvars[n][0] = ( vars[n][CCTK_GFINDEX3D(cctkGH, i-2,j,k)] - vars[n][CCTK_GFINDEX3D(cctkGH, i+2,j,k)] + + 8.0 * (vars[n][CCTK_GFINDEX3D(cctkGH, i+1,j,k)] - vars[n][CCTK_GFINDEX3D(cctkGH, i-1,j,k)])) * (1.0 / 12.0) * ih[0]; + dvars[n][1] = ( vars[n][CCTK_GFINDEX3D(cctkGH, i,j-2,k)] - vars[n][CCTK_GFINDEX3D(cctkGH, i,j+2,k)] + + 8.0 * (vars[n][CCTK_GFINDEX3D(cctkGH, i,j+1,k)] - vars[n][CCTK_GFINDEX3D(cctkGH, i,j-1,k)])) * (1.0 / 12.0) * ih[1]; + dvars[n][2] = ( vars[n][CCTK_GFINDEX3D(cctkGH, i,j,k-2)] - vars[n][CCTK_GFINDEX3D(cctkGH, i,j,k+2)] + + 8.0 * (vars[n][CCTK_GFINDEX3D(cctkGH, i,j,k+1)] - vars[n][CCTK_GFINDEX3D(cctkGH, i,j,k-1)])) * (1.0 / 12.0) * ih[2]; + } + } +}; + + +static inline void UpperMetric(CCTK_REAL& restrict uxx, + CCTK_REAL& restrict uxy, + CCTK_REAL& restrict uxz, + CCTK_REAL& restrict uyy, + CCTK_REAL& restrict uyz, + CCTK_REAL& restrict uzz, + const CCTK_REAL det, + const CCTK_REAL gxx, + const CCTK_REAL gxy, + const CCTK_REAL gxz, + const CCTK_REAL gyy, + const CCTK_REAL gyz, + const CCTK_REAL gzz) +{ + const CCTK_REAL invdet = 1.0 / det; + uxx = (-gyz*gyz + gyy*gzz)*invdet; + uxy = (gxz*gyz - gxy*gzz)*invdet; + uxz = (-gxz*gyy + gxy*gyz)*invdet; + uyy = (-gxz*gxz + gxx*gzz)*invdet; + uyz = (gxy*gxz - gxx*gyz)*invdet; + uzz = (-gxy*gxy + gxx*gyy)*invdet; + + return; +} + + +template +static void SourceTerms_LL(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + /* + Set up multipatch stuff + */ + + const CCTK_REAL * restrict vup; + const CCTK_REAL * restrict Bprim; + const CCTK_REAL * restrict g11; + const CCTK_REAL * restrict g12; + const CCTK_REAL * restrict g13; + const CCTK_REAL * restrict g22; + const CCTK_REAL * restrict g23; + const CCTK_REAL * restrict g33; + const CCTK_REAL * restrict k11; + const CCTK_REAL * restrict k12; + const CCTK_REAL * restrict k13; + const CCTK_REAL * restrict k22; + const CCTK_REAL * restrict k23; + const CCTK_REAL * restrict k33; + const CCTK_REAL * restrict beta1; + const CCTK_REAL * restrict beta2; + const CCTK_REAL * restrict beta3; + + //Multipatch related pointers + if(GRHydro_UseGeneralCoordinates(cctkGH)) { + vup=lvel; + g11=gaa; g12=gab; g13=gac; + g22=gbb; g23=gbc; + g33=gcc; + k11=kaa; k12=kab; k13=kac; + k22=kbb; k23=kbc; + k33=kcc; + Bprim=lBvec; + beta1 = betaa; + beta2 = betab; + beta3 = betac; + } else { + vup=vel; + g11=gxx; g12=gxy; g13=gxz; + g22=gyy; g23=gyz; + g33=gzz; + k11=kxx; k12=kxy; k13=kxz; + k22=kyy; k23=kyz; + k33=kzz; + Bprim=Bvec; + beta1 = betax; + beta2 = betay; + beta3 = betaz; + } + + // TODO: for vectorization we might want to use cctk_ash (have to make sure + // sensible values exist everywhere) + const int nx=cctk_lsh[0]; + const int ny=cctk_lsh[1]; + const int nz=cctk_lsh[2]; + + const int N = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + + const CCTK_REAL one = 1.00; + const CCTK_REAL two = 2.00; + const CCTK_REAL half = 0.50; + const CCTK_REAL dx = CCTK_DELTA_SPACE(0); + const CCTK_REAL dy = CCTK_DELTA_SPACE(1); + const CCTK_REAL dz = CCTK_DELTA_SPACE(2); + + const CCTK_REAL ih[3] = { 1.0/dx, 1.0/dy, 1.0/dz }; + +#pragma omp parallel for + for (int k=GRHydro_stencil; k < nz-GRHydro_stencil; ++k) { + for (int j=GRHydro_stencil; j < ny-GRHydro_stencil; ++j) { + for (int i=GRHydro_stencil; i < nx-GRHydro_stencil; ++i) { + + const int ijk = CCTK_GFINDEX3D(cctkGH, i, j, k); + + // TODO: these are not needed. + const CCTK_REAL localgxx = g11[ijk]; + const CCTK_REAL localgxy = g12[ijk]; + const CCTK_REAL localgxz = g13[ijk]; + const CCTK_REAL localgyy = g22[ijk]; + const CCTK_REAL localgyz = g23[ijk]; + const CCTK_REAL localgzz = g33[ijk]; + + const CCTK_REAL sqrtdet = sdetg[ijk]; + const CCTK_REAL invsqrtdet = 1./sqrtdet; + + CCTK_REAL uxx, uxy, uxz, uyy, uyz, uzz; + UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, sqrtdet*sqrtdet, localgxx, + localgxy, localgxz, localgyy, localgyz, localgzz); + + const CCTK_REAL shiftx = beta1[ijk]; + const CCTK_REAL shifty = beta2[ijk]; + const CCTK_REAL shiftz = beta3[ijk]; + + // Derivatives of the lapse, metric and shift + + const int nvars = 10; + const CCTK_REAL* const restrict vars[nvars] = { beta1, beta2, beta3, alp, + g11, g12, g13, + g22, g23, + g33 }; + CCTK_REAL dvars[nvars][3]; + + alldiff::apply(cctkGH, dvars, vars, i, j, k, ih, nvars); + + const CCTK_REAL dx_betax = dvars[0][0]; + const CCTK_REAL dx_betay = dvars[1][0]; + const CCTK_REAL dx_betaz = dvars[2][0]; + + const CCTK_REAL dy_betax = dvars[0][1]; + const CCTK_REAL dy_betay = dvars[1][1]; + const CCTK_REAL dy_betaz = dvars[2][1]; + + const CCTK_REAL dz_betax = dvars[0][2]; + const CCTK_REAL dz_betay = dvars[1][2]; + const CCTK_REAL dz_betaz = dvars[2][2]; + + const CCTK_REAL dx_alp = dvars[3][0]; + const CCTK_REAL dy_alp = dvars[3][1]; + const CCTK_REAL dz_alp = dvars[3][2]; + + const CCTK_REAL dx_gxx = dvars[4][0]; + const CCTK_REAL dx_gxy = dvars[5][0]; + const CCTK_REAL dx_gxz = dvars[6][0]; + const CCTK_REAL dx_gyy = dvars[7][0]; + const CCTK_REAL dx_gyz = dvars[8][0]; + const CCTK_REAL dx_gzz = dvars[9][0]; + const CCTK_REAL dy_gxx = dvars[4][1]; + const CCTK_REAL dy_gxy = dvars[5][1]; + const CCTK_REAL dy_gxz = dvars[6][1]; + const CCTK_REAL dy_gyy = dvars[7][1]; + const CCTK_REAL dy_gyz = dvars[8][1]; + const CCTK_REAL dy_gzz = dvars[9][1]; + const CCTK_REAL dz_gxx = dvars[4][2]; + const CCTK_REAL dz_gxy = dvars[5][2]; + const CCTK_REAL dz_gxz = dvars[6][2]; + const CCTK_REAL dz_gyy = dvars[7][2]; + const CCTK_REAL dz_gyz = dvars[8][2]; + const CCTK_REAL dz_gzz = dvars[9][2]; + + + const CCTK_REAL invalp = 1.0 / alp[ijk]; + const CCTK_REAL invalp2 = SQR(invalp); + const CCTK_REAL velxshift = velx[ijk] - shiftx*invalp; + const CCTK_REAL velyshift = vely[ijk] - shifty*invalp; + const CCTK_REAL velzshift = velz[ijk] - shiftz*invalp; + + // vel_i = g_ij v^j + // B_i = g_ij B^i + + const CCTK_REAL vlowx = g11[ijk]*velx[ijk] + g12[ijk]*vely[ijk] + g13[ijk]*velz[ijk]; + const CCTK_REAL vlowy = g12[ijk]*velx[ijk] + g22[ijk]*vely[ijk] + g23[ijk]*velz[ijk]; + const CCTK_REAL vlowz = g13[ijk]*velx[ijk] + g23[ijk]*vely[ijk] + g33[ijk]*velz[ijk]; + const CCTK_REAL Bvecxlow = do_mhd ? g11[ijk]*Bvecx[ijk] + g12[ijk]*Bvecy[ijk] + g13[ijk]*Bvecz[ijk] : 0; + const CCTK_REAL Bvecylow = do_mhd ? g12[ijk]*Bvecx[ijk] + g22[ijk]*Bvecy[ijk] + g23[ijk]*Bvecz[ijk] : 0; + const CCTK_REAL Bveczlow = do_mhd ? g13[ijk]*Bvecx[ijk] + g23[ijk]*Bvecy[ijk] + g33[ijk]*Bvecz[ijk] : 0; + + // B^i v_i (= b^0/u^0) + const CCTK_REAL Bdotv = do_mhd ? vlowx*Bvecx[ijk]+vlowy*Bvecy[ijk]+vlowz*Bvecz[ijk] : 0; + + // v^2 = v_i v^i; w=(1-v^2)^{-1/2} + + const CCTK_REAL v2 = vlowx*velx[ijk] + vlowy*vely[ijk] + vlowz*velz[ijk]; + const CCTK_REAL invw = sqrt(1.0-v2); + const CCTK_REAL w = 1./invw; + + // b^2 = B^i B_i / w^2 + (b^0/u^0)^2 + + const CCTK_REAL b2 = do_mhd ? (Bvecx[ijk]*Bvecxlow+Bvecy[ijk]*Bvecylow+Bvecz[ijk]*Bveczlow)*SQR(invw)+SQR(Bdotv) : 0; + + // b_i = B_i/w +w*(B dot v)*v_i + const CCTK_REAL bxlow = do_mhd ? Bvecxlow*invw+w*Bdotv*vlowx : 0; + const CCTK_REAL bylow = do_mhd ? Bvecylow*invw+w*Bdotv*vlowy : 0; + const CCTK_REAL bzlow = do_mhd ? Bveczlow*invw+w*Bdotv*vlowz : 0; + + + // These are the contravariant components + const CCTK_REAL bt = do_mhd ? w*invalp*Bdotv : 0; + const CCTK_REAL bx = do_mhd ? Bvecx[ijk]*invw+w*Bdotv*velxshift : 0; + const CCTK_REAL by = do_mhd ? Bvecy[ijk]*invw+w*Bdotv*velyshift : 0; + const CCTK_REAL bz = do_mhd ? Bvecz[ijk]*invw+w*Bdotv*velzshift : 0; + + // TODO: all of these can be expressed much more easily in terms of the + // conservatives + const CCTK_REAL rhohstarW2 = (rho[ijk]*(one + eps[ijk]) + (press[ijk] + b2)) * + SQR(w); + const CCTK_REAL pstar = press[ijk]+0.50*b2; + + // For a change, these are T^{ij} + + // TODO: strict IEEE compliance does not let the compiler remove "+ 0" + // terms, so we have to do something else here + const CCTK_REAL t00 = (rhohstarW2 - pstar)*invalp2-SQR(bt); + const CCTK_REAL t0x = rhohstarW2*velxshift*invalp + + pstar*shiftx*invalp2-bt*bx; + const CCTK_REAL t0y = rhohstarW2*velyshift*invalp + + pstar*shifty*invalp2-bt*by; + const CCTK_REAL t0z = rhohstarW2*velzshift*invalp + + pstar*shiftz*invalp2-bt*bz; + const CCTK_REAL txx = rhohstarW2*velxshift*velxshift + + pstar*(uxx - shiftx*shiftx*invalp2)-SQR(bx); + const CCTK_REAL txy = rhohstarW2*velxshift*velyshift + + pstar*(uxy - shiftx*shifty*invalp2)-bx*by; + const CCTK_REAL txz = rhohstarW2*velxshift*velzshift + + pstar*(uxz - shiftx*shiftz*invalp2)-bx*bz; + const CCTK_REAL tyy = rhohstarW2*velyshift*velyshift + + pstar*(uyy - shifty*shifty*invalp2)-SQR(by); + const CCTK_REAL tyz = rhohstarW2*velyshift*velzshift + + pstar*(uyz - shifty*shiftz*invalp2)-by*bz; + const CCTK_REAL tzz = rhohstarW2*velzshift*velzshift + + pstar*(uzz - shiftz*shiftz*invalp2)-SQR(bz); + +// Contract the shift with the extrinsic curvature + + const CCTK_REAL shiftshiftk = shiftx*shiftx*k11[ijk] + + shifty*shifty*k22[ijk] + + shiftz*shiftz*k33[ijk] + + two*(shiftx*shifty*k12[ijk] + + shiftx*shiftz*k13[ijk] + + shifty*shiftz*k23[ijk]); + + const CCTK_REAL shiftkx = shiftx*k11[ijk] + shifty*k12[ijk] + shiftz*k13[ijk]; + const CCTK_REAL shiftky = shiftx*k12[ijk] + shifty*k22[ijk] + shiftz*k23[ijk]; + const CCTK_REAL shiftkz = shiftx*k13[ijk] + shifty*k23[ijk] + shiftz*k33[ijk]; + +// Contract the matter terms with the extrinsic curvature + + const CCTK_REAL sumTK = txx*k11[ijk] + tyy*k22[ijk] + tzz*k33[ijk] + + two*(txy*k12[ijk] + txz*k13[ijk] + tyz*k23[ijk]); + +// Update term for tau + + const CCTK_REAL tau_source = t00* + (shiftshiftk - (shiftx*dx_alp + shifty*dy_alp + shiftz*dz_alp) ) + + t0x*(-dx_alp + two*shiftkx) + + t0y*(-dy_alp + two*shiftky) + + t0z*(-dz_alp + two*shiftkz) + + sumTK; + +// The following looks very little like the terms in the +// standard papers. Take a look in the ThornGuide to see why +// it is really the same thing. + +// Contract the shift with derivatives of the metric + + const CCTK_REAL halfshiftdgx = half*(shiftx*shiftx*dx_gxx + + shifty*shifty*dx_gyy + shiftz*shiftz*dx_gzz) + + shiftx*shifty*dx_gxy + shiftx*shiftz*dx_gxz + + shifty*shiftz*dx_gyz; + const CCTK_REAL halfshiftdgy = half*(shiftx*shiftx*dy_gxx + + shifty*shifty*dy_gyy + shiftz*shiftz*dy_gzz) + + shiftx*shifty*dy_gxy + shiftx*shiftz*dy_gxz + + shifty*shiftz*dy_gyz; + const CCTK_REAL halfshiftdgz = half*(shiftx*shiftx*dz_gxx + + shifty*shifty*dz_gyy + shiftz*shiftz*dz_gzz) + + shiftx*shifty*dz_gxy + shiftx*shiftz*dz_gxz + + shifty*shiftz*dz_gyz; + +// Contract the matter with derivatives of the metric + + const CCTK_REAL halfTdgx = half*(txx*dx_gxx + tyy*dx_gyy + tzz*dx_gzz) + + txy*dx_gxy + txz*dx_gxz + tyz*dx_gyz; + const CCTK_REAL halfTdgy = half*(txx*dy_gxx + tyy*dy_gyy + tzz*dy_gzz) + + txy*dy_gxy + txz*dy_gxz + tyz*dy_gyz; + const CCTK_REAL halfTdgz = half*(txx*dz_gxx + tyy*dz_gyy + tzz*dz_gzz) + + txy*dz_gxy + txz*dz_gxz + tyz*dz_gyz; + + + const CCTK_REAL sx_source = t00* + (halfshiftdgx - alp[ijk]*dx_alp) + halfTdgx + + t0x*(shiftx*dx_gxx + shifty*dx_gxy + shiftz*dx_gxz) + + t0y*(shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) + + t0z*(shiftx*dx_gxz + shifty*dx_gyz + shiftz*dx_gzz) + + rhohstarW2*invalp*(vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz) - + bt*(bxlow*dx_betax + bylow*dx_betay + bzlow*dx_betaz); + + const CCTK_REAL sy_source = t00* + (halfshiftdgy - alp[ijk]*dy_alp) + halfTdgy + + t0x*(shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) + + t0y*(shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) + + t0z*(shiftx*dy_gxz + shifty*dy_gyz + shiftz*dy_gzz) + + rhohstarW2*invalp*(vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz) - + bt*(bxlow*dy_betax + bylow*dy_betay + bzlow*dy_betaz); + + const CCTK_REAL sz_source = t00* + (halfshiftdgz - alp[ijk]*dz_alp) + halfTdgz + + t0x*(shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) + + t0y*(shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) + + t0z*(shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) + + rhohstarW2*invalp*(vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz) - + bt*(bxlow*dz_betax + bylow*dz_betay + bzlow*dz_betaz); + + densrhs[ijk] = 0.0; + srhs[ijk] = alp[ijk]*sqrtdet*sx_source; + srhs[ijk + N] = alp[ijk]*sqrtdet*sy_source; + srhs[ijk + 2*N] = alp[ijk]*sqrtdet*sz_source; + taurhs[ijk] = alp[ijk]*sqrtdet*tau_source; + + if (do_Avec) { + + // B^i and A^i both live in cell centers currently + const CCTK_REAL Avcx_source = alp[ijk]*sqrtdet*(velyshift*Bvecz[ijk] - velzshift*Bvecy[ijk]); + const CCTK_REAL Avcy_source = alp[ijk]*sqrtdet*(velzshift*Bvecx[ijk] - velxshift*Bvecz[ijk]); + const CCTK_REAL Avcz_source = alp[ijk]*sqrtdet*(velxshift*Bvecy[ijk] - velyshift*Bvecx[ijk]); + + Avecrhsx[ijk] = Avcx_source; + Avecrhsy[ijk] = Avcy_source; + Avecrhsz[ijk] = Avcz_source; + + } + + if(do_clean_divergence) { + + // g^{jk} d_i g_{kj} = d_i (g) / det + const CCTK_REAL dx_det_bydet = uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + + two*(uxy*dx_gxy+uxz*dx_gxz+uyz*dx_gyz); + const CCTK_REAL dy_det_bydet = uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + + two*(uxy*dy_gxy+uxz*dy_gxz+uyz*dy_gyz); + const CCTK_REAL dz_det_bydet = uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + + two*(uxy*dz_gxy+uxz*dz_gxz+uyz*dz_gyz); + + // g^{ik} d_k g_{li} + const CCTK_REAL gdg_x = uxx*dx_gxx + uxy*dy_gxx + uxz*dz_gxx + + uxy*dx_gxy + uyy*dy_gxy + uyz*dz_gxy + + uxz*dx_gxz + uyz*dy_gxz + uzz*dz_gxz; + + const CCTK_REAL gdg_y = uxx*dx_gxy + uxy*dy_gxy + uxz*dz_gxy + + uxy*dx_gyy + uyy*dy_gyy + uyz*dz_gyy + + uxz*dx_gyz + uyz*dy_gyz + uzz*dz_gyz; + + const CCTK_REAL gdg_z = uxx*dx_gxz + uxy*dy_gxz + uxz*dz_gxz + + uxy*dx_gyz + uyy*dy_gyz + uyz*dz_gyz + + uxz*dx_gzz + uyz*dy_gzz + uzz*dz_gzz; + + CCTK_REAL bvcx_source, bvcy_source, bvcz_source; + psidcrhs[ijk] = -one * (kap_dc*alp[ijk] + + dx_betax + dy_betay + dz_betaz ) * psidc[ijk] + + Bconsx[ijk] * (dx_alp - half*alp[ijk] * + ( uxx*dx_gxx + uyy*dx_gyy + uzz*dx_gzz + two*uxy*dx_gxy + + two*uxz*dx_gxz + two*uyz*dx_gyz ) )*invsqrtdet + + Bconsy[ijk] * (dy_alp - half*alp[ijk] * + ( uxx*dy_gxx + uyy*dy_gyy + uzz*dy_gzz + two*uxy*dy_gxy + + two*uxz*dy_gxz + two*uyz*dy_gyz ) )*invsqrtdet + + Bconsz[ijk] * (dz_alp - half*alp[ijk] * + ( uxx*dz_gxx + uyy*dz_gyy + uzz*dz_gzz + two*uxy*dz_gxy + + two*uxz*dz_gxz + two*uyz*dz_gyz ) )*invsqrtdet; + + bvcx_source = -one * ( Bconsx[ijk]*dx_betax + + Bconsy[ijk]*dy_betax + Bconsz[ijk]*dz_betax ) + + psidc[ijk]*sqrtdet*(( uxx*dx_alp+uxy*dy_alp+uxz*dz_alp ) + + alp[ijk]*(half*( uxx*dx_det_bydet + + uxy*dy_det_bydet + uxz*dz_det_bydet) - + ( uxx*gdg_x + uxy*gdg_y + uxz*gdg_z ))); + + bvcy_source = -one * ( Bconsx[ijk]*dx_betay + + Bconsy[ijk]*dy_betay + Bconsz[ijk]*dz_betay ) + + psidc[ijk]*sqrtdet*(( uxy*dx_alp+uyy*dy_alp+uyz*dz_alp ) + + alp[ijk]*(half*( uxy*dx_det_bydet + + uyy*dy_det_bydet + uyz*dz_det_bydet ) - + ( uxy*gdg_x + uyy*gdg_y + uyz*gdg_z ))); + + bvcz_source = -one * ( Bconsx[ijk]*dx_betaz + + Bconsy[ijk]*dy_betaz + Bconsz[ijk]*dz_betaz ) + + psidc[ijk]*sqrtdet*(( uxz*dx_alp+uyz*dy_alp+uzz*dz_alp ) + + alp[ijk]*(half*( uxz*dx_det_bydet + + uyz*dy_det_bydet + uzz*dz_det_bydet ) - + ( uxz*gdg_x + uyz*gdg_y + uzz*gdg_z ))); + + Bconsrhsx[ijk] = bvcx_source; + Bconsrhsy[ijk] = bvcy_source; + Bconsrhsz[ijk] = bvcz_source; + } // if(do_clean_divergence) + + } // for(i + } // for(j + } // for(k + +#if(0) // poison edges of domain + static int last_iteration_seen = -1; + static int reflevel = -1; + static int mol_substep; + if(last_iteration_seen != cctk_iteration || reflevel != *GRHydro_reflevel) { + last_iteration_seen = cctk_iteration; + reflevel = *GRHydro_reflevel; + mol_substep = 0; + } else { + mol_substep += 1; + } + for(k = 0 ; k < GRHydro_stencil*mol_substep ; ++k) { + for(j = 0 ; j < cctk_ash[1] ; ++j) { + for(i = 0 ; i < cctk_ash[0] ; ++i) { + dens[CCTK_GFINDEX3D(i,j,k)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + tau[CCTK_GFINDEX3D(i,j,k)] = -1e100 + if(do_mhd) { + Bcons[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + } + if(do_clean_divergene) + psidc[CCTK_GFINDEX3D(i,j,k)] = -1e100 + } + } + } + for(k = cctk_ash[2]-GRHydro_stencil*mol_substep ; k < cctk_ash[2] ; ++k) { + for(j = 0 ; j < cctk_ash[1] ; ++j) { + for(i = 0 ; i < cctk_ash[0] ; ++i) { + dens[CCTK_GFINDEX3D(i,j,k)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + tau[CCTK_GFINDEX3D(i,j,k)] = -1e100 + if(do_mhd) { + Bcons[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + } + if(do_clean_divergene) + psidc[CCTK_GFINDEX3D(i,j,k)] = -1e100 + } + } + } + for(i = 0 ; i < GRHydro_stencil*mol_substep ; ++i) { + for(k = 0 ; k < cctk_ash[2] ; ++k) { + for(j = 0 ; j < cctk_ash[1] ; ++j) { + dens[CCTK_GFINDEX3D(i,j,k)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + tau[CCTK_GFINDEX3D(i,j,k)] = -1e100 + if(do_mhd) { + Bcons[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + } + if(do_clean_divergene) + psidc[CCTK_GFINDEX3D(i,j,k)] = -1e100 + } + } + } + for(i = cctk_ash[0]-GRHydro_stencil*mol_substep ; i < cctk_ash[0] ; ++i) { + for(k = 0 ; k < cctk_ash[2] ; ++k) { + for(j = 0 ; j < cctk_ash[1] ; ++j) { + dens[CCTK_GFINDEX3D(i,j,k)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + tau[CCTK_GFINDEX3D(i,j,k)] = -1e100 + if(do_mhd) { + Bcons[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + } + if(do_clean_divergene) + psidc[CCTK_GFINDEX3D(i,j,k)] = -1e100 + } + } + } + for(j = 0 ; j < GRHydro_stencil*mol_substep ; ++j) { + for(i = 0 ; i < cctk_ash[0] ; ++i) { + for(k = 0 ; k < cctk_ash[2] ; ++k) { + dens[CCTK_GFINDEX3D(i,j,k)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + tau[CCTK_GFINDEX3D(i,j,k)] = -1e100 + if(do_mhd) { + Bcons[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + } + if(do_clean_divergene) + psidc[CCTK_GFINDEX3D(i,j,k)] = -1e100 + } + } + } + for(j = cctk_ash[1]-GRHydro_stencil*mol_substep ; j < cctk_ash[1] ; ++j) { + for(i = 0 ; i < cctk_ash[0] ; ++i) { + for(k = 0 ; k < cctk_ash[2] ; ++k) { + dens[CCTK_GFINDEX3D(i,j,k)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Scon[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + tau[CCTK_GFINDEX3D(i,j,k)] = -1e100 + if(do_mhd) { + Bcons[CCTK_VECTGFINDEX3D(i,j,k,0)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,1)] = -1e100 + Bcons[CCTK_VECTGFINDEX3D(i,j,k,2)] = -1e100 + } + if(do_clean_divergene) + psidc[CCTK_GFINDEX3D(i,j,k)] = -1e100 + } + } + } +#endif +} + +extern "C" void SourceTerms(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + void (*sources)(CCTK_ARGUMENTS) = NULL; + const bool do_Avec = CCTK_EQUALS(Bvec_evolution_method, "GRHydro_Avec"); +#define T(e,t) (!!(e)==(t)) +#define T3(t1,t2,t3) \ + (T(*evolve_MHD,t1) && \ + T(clean_divergence,t2) && \ + T(do_Avec,t3)) +#define SET_SOURCES(df,t1,t2,t3) \ + if(T3(t1,t2,t3)) \ + sources=SourceTerms_LL; + enum do_t {noMHD=0,MHD=1,noDC=0,DC=1,noAVec=0,AVec=1}; + + // order of arguments: + // do_MHD do_clean_divergence do_Avec + if(spatial_order == 4) { + // no MHD + SET_SOURCES(alldiff4,noMHD,noDC,noAVec) + + // not all combinations make sense for MHD, we list the ones that do + SET_SOURCES(alldiff4, MHD,noDC,noAVec) + SET_SOURCES(alldiff4, MHD, DC,noAVec) + SET_SOURCES(alldiff4, MHD, DC,noAVec) + SET_SOURCES(alldiff4, MHD,noDC, AVec) + } else if(spatial_order == 2) { + // no MHD + SET_SOURCES(alldiff2,noMHD,noDC,noAVec) + + // not all combinations make sense for MHD, we list the ones that do + SET_SOURCES(alldiff2, MHD,noDC,noAVec) + SET_SOURCES(alldiff2, MHD, DC,noAVec) + SET_SOURCES(alldiff2, MHD, DC,noAVec) + SET_SOURCES(alldiff2, MHD,noDC, AVec) + } + + if(sources == NULL) { + CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, + "No source terms can be computed for combination " + "do_mhd = %d " + "do_clean_divergence = %d " + "do_Avec = %d spatial_order = %d", + int(*evolve_MHD), + int(clean_divergence), + int(do_Avec), int(spatial_order)); + } + + // do some non-time critial operations here to reduce the number of templates + if(transport_constraints) { + const int N = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + // not much need to OpenMP parallelize this since it is pure memory access + // and memset is highly optimized + memset(Evec, 0, sizeof(Evec[0])*3*N); + // I think IEEE requires an all zero byte double to be 0. + assert(Evec[0] == 0.); + } + if(*evolve_Y_e) { + const int N = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + // not much need to OpenMP parallelize this since it is pure memory access + // and memset is highly optimized + memset(Y_e_con_rhs, 0, sizeof(Y_e_con_rhs[0])*N); + // I think IEEE requires an all zero byte double to be 0. + assert(Y_e_con_rhs[0] == 0.); + } + if(number_of_tracers>0) { + const int N = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + // not much need to OpenMP parallelize this since it is pure memory access + // and memset is highly optimized + memset(cons_tracerrhs, 0, sizeof(cons_tracerrhs[0])*number_of_tracers*N); + // I think IEEE requires an all zero byte double to be 0. + assert(cons_tracerrhs[0] == 0.); + } + if(*evolve_entropy) { + const int N = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + // not much need to OpenMP parallelize this since it is pure memory access + // and memset is highly optimized + memset(entropyrhs, 0, sizeof(entropyrhs[0])*N); + // I think IEEE requires an all zero byte double to be 0. + assert(entropyrhs[0] == 0.); + } + if(track_divB) { + const int N = cctk_ash[0]*cctk_ash[1]*cctk_ash[2]; + // not much need to OpenMP parallelize this since it is pure memory access + // and memset is highly optimized + memset(divB, 0, sizeof(divB[0])*N); + // I think IEEE requires an all zero byte double to be 0. + assert(divB[0] == 0.); + } + + sources(CCTK_PASS_CTOC); +} + + diff --git a/src/make.code.defn b/src/make.code.defn index 8739570..009e5f0 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -37,7 +37,6 @@ SRCS = Utils.F90 \ GRHydro_RoeSolver.F90 \ GRHydro_Scalars.F90 \ GRHydro_Shift.F90 \ - GRHydro_Source.F90 \ GRHydro_Startup.F90 \ GRHydro_TVDReconstruct.F90 \ GRHydro_Marquina.F90 \ @@ -92,8 +91,11 @@ SRCS = Utils.F90 \ GRHydro_TrivialReconstruct.cc \ GRHydro_MP5Reconstruct.cc \ GRHydro_Wrappers.F90 \ - GRHydro_HLLE.cc - + GRHydro_HLLE.cc \ + GRHydro_Source.cc \ + GRHydro_Source.F90 + + # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3