User: rhaas
Date: 2014/04/15 02:49 PM
Added:
/trunk/src/
GRHydro_Source.cc
Modified:
/trunk/
param.ccl, schedule.ccl
/trunk/src/
make.code.defn
Log:
GRHydro: implement sources in C++.
Combined result of commits by Christian Reisswig and Roland Haas.
File Changes:
Directory: /trunk/src/
======================
File [added]: GRHydro_Source.cc
Delta lines: +709 -0
===================================================================
--- trunk/src/GRHydro_Source.cc (rev 0)
+++ trunk/src/GRHydro_Source.cc 2014-04-15 19:49:34 UTC (rev 627)
@@ -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 <cstring>
+#include <cassert>
+
+#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 <typename T> 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 <class alldiff,
+ bool do_mhd,
+ bool do_clean_divergence,
+ bool do_Avec>
+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<df,t1,t2,t3>;
+ 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);
+}
+
+
File [modified]: make.code.defn
Delta lines: +5 -3
===================================================================
--- trunk/src/make.code.defn 2014-04-15 19:49:31 UTC (rev 626)
+++ trunk/src/make.code.defn 2014-04-15 19:49:34 UTC (rev 627)
@@ -37,7 +37,6 @@
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 @@
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 =
Directory: /trunk/
==================
File [modified]: param.ccl
Delta lines: +0 -1
===================================================================
--- trunk/param.ccl 2014-04-15 19:49:31 UTC (rev 626)
+++ trunk/param.ccl 2014-04-15 19:49:34 UTC (rev 627)
@@ -810,4 +810,3 @@
BOOLEAN use_cxx_code "Use experimental C++ code?" STEERABLE = RECOVER
{
} no
-
File [modified]: schedule.ccl
Delta lines: +20 -13
===================================================================
--- trunk/schedule.ccl 2014-04-15 19:49:31 UTC (rev 626)
+++ trunk/schedule.ccl 2014-04-15 19:49:34 UTC (rev 627)
@@ -732,22 +732,29 @@
### over the different directions ###
####################################################################
-if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
-{
- schedule SourceTermsM IN GRHydroRHS BEFORE FluxTerms
+if(use_cxx_code) {
+ schedule SourceTerms 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"
+ LANG: C
+ } "Source term calculation"
} else {
- schedule SourceTerms IN GRHydroRHS BEFORE FluxTerms
+ if(CCTK_Equals(Bvec_evolution_method,"GRHydro"))
{
- LANG: Fortran
- } "Source term calculation"
+ 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. ###
