diff options
| -rw-r--r-- | README | 20 | ||||
| -rw-r--r-- | interface.ccl | 27 | ||||
| -rw-r--r-- | param.ccl | 315 | ||||
| -rw-r--r-- | schedule.ccl | 2 | ||||
| -rw-r--r-- | src/Comparison.c | 331 | ||||
| -rw-r--r-- | src/ComparisonSolutions.F | 183 | ||||
| -rw-r--r-- | src/boostrotmetric.F | 165 | ||||
| -rw-r--r-- | src/bowl.F | 251 | ||||
| -rw-r--r-- | src/exactblendbound.F | 193 | ||||
| -rw-r--r-- | src/exactboundary.F | 123 | ||||
| -rw-r--r-- | src/exactcartblendbound.F | 207 | ||||
| -rw-r--r-- | src/exactdata.F | 194 | ||||
| -rw-r--r-- | src/exactgauge.F | 135 | ||||
| -rw-r--r-- | src/exactinitialize.F | 71 | ||||
| -rw-r--r-- | src/exactmetric.F | 136 | ||||
| -rw-r--r-- | src/fakebinary.F | 173 | ||||
| -rw-r--r-- | src/finkelstein.F | 66 | ||||
| -rw-r--r-- | src/flatschwarz.F | 70 | ||||
| -rw-r--r-- | src/include/pGF_ComparisonExtensions.h | 1 | ||||
| -rw-r--r-- | src/include/slice_normal.h | 103 | ||||
| -rw-r--r-- | src/kerrschild.F | 127 | ||||
| -rw-r--r-- | src/linextraponebound.F | 140 | ||||
| -rw-r--r-- | src/make.code.defn | 15 | ||||
| -rw-r--r-- | src/minkowski.F | 46 | ||||
| -rw-r--r-- | src/multiBH.F | 133 | ||||
| -rw-r--r-- | src/novikov.F | 201 | ||||
| -rw-r--r-- | src/slice_data.F | 327 | ||||
| -rw-r--r-- | src/slice_evolve.F | 132 | ||||
| -rw-r--r-- | src/slice_initialize.F | 41 | ||||
| -rw-r--r-- | src/starschwarz.F | 120 |
30 files changed, 4048 insertions, 0 deletions
@@ -0,0 +1,20 @@ +Cactus Code Thorn Exact +Authors : Carsten Gundlach initially, plus many other authors. +CVS info : $Header$ +-------------------------------------------------------------------------- + +Purpose of the thorn: + +The initial purpose of this thorn was the comparison of numerical +spacetimes against exact solutions. One particular solution, +the boost-rotation-symmetric spacetimes, was the first application. + +Later, the thorn developed as a reservoir for a variety of exact spacetimes, +some of them in different coordinate systems, and even some non-Einstein +spcetimes. All od these exact spacetimes have been found useful for +testing different aspect of the code. Since mamy different people have +contributed to this thorn by adding new exact solutions, you shoudl expect +many different styles of coding here. + +For more details see the documentation in the directory doc. + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..64eacd0 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,27 @@ +# Interface definition for thorn Exact +# $Header$ + +implements: exact +inherits: einstein + +private: + +real Exact_slice type=GF +{ + slicex, + slicey, + slicez, + slicet +} "position of an arbitrary slice in exact solution spacetime" + +real Exact_slicetemp type=GF +{ + slicetmp1x, + slicetmp1y, + slicetmp1z, + slicetmp1t, + slicetmp2x, + slicetmp2y, + slicetmp2z, + slicetmp2t +} "Temporary grid functions"
\ No newline at end of file diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..a37da6a --- /dev/null +++ b/param.ccl @@ -0,0 +1,315 @@ +# Parameter definitions for thorn Exact +# $Header$ + +shares: grid + +shares: einstein + +EXTENDS KEYWORD slicing "" +{ + "exact" :: Use lapse from exact solution +} "" + +EXTENDS KEYWORD shift "" +{ + "exact" :: Use shift from exact solution +} "" + + +restricted: + +# General parameters + +KEYWORD exactmodel "The exact solution used in thorn exact" +{ + "minkowski" :: Minkowski spacetime + "boostrot" :: Boost rotation symmetric spacetime + "finkelstein" :: Black hole in Eddington-Finkelstein coordinates + "kerrschild" :: Kerr-Schild form of boosted rotating black hole + "flatschwarz" :: Schwazschild black-hole with flat spatial metric + "starschwarz" :: Schwazschild (constant density) star + "novikov" :: Black hole in Novikov coordinates + "multiBH" :: Maximally charged multi BH solutions + "bowl" :: Non-Einstein bowl (bag-of-gold) spacetime + "fakebinary" :: Non-Einstein fake-binary of Thorn et al +} "minkowski" + + +# Parameters for the blended boundaries. + +BOOLEAN exblend_Ks "Blend the K variables with the exact solution?" +{ +} "yes" + +BOOLEAN exblend_gs "Blend the g variables with the exact solution?" +{ +} "yes" + +BOOLEAN exblend_gauge "Blend the lapse and shift with the exact solution?" +{ +} "yes" + +REAL exblend_rout "Outer boundary of blending region" +{ +: :: "Positive means radial value, negative means use outer bound of grid" +} -1.0 + +REAL exblend_width "Width of blending zone" +{ +: :: Positive means width in radius, negative means width = exbeldn_width*dx" +} -3.0 + + +# Parameters for slice + +REAL slice_gauss_ampl "Amplitude of Gauss slice in exact" +{ +0.0: :: "Positive please" +} 0.0 + +REAL slice_gauss_width "Width of Gauss slice in exact" +{ +0.0: :: "Positive please" +} 1.0 + + +# Parameters for boostrot + +REAL boostrotscale "Length scale of boost rotation data" +{ +0.0: :: "Positive please" +} 1.0 + +REAL boostrotstrength "Dimensionless strength parameter" +{ +0.0: :: "Positive please" +} 0.1 + +REAL boostrotsafedistance "Dimensionless safety distance" +{ +0.0: :: "Positive please" +} 0.01 + + +# Parameters for finkelstein and kerrschild. + +REAL kerrschild_boostv "Boost speed of black hole" +{ +: :: "Positive or negative" +} 0.0 + +REAL kerrschild_eps "Fudge parameter" +{ +: :: "Positive or negative" +} 1.e-16 + +REAL kerrschild_m "Black hole mass" +{ +0.0: :: "Positive please" +} 1.0 + +REAL kerrschild_a "Black hole a = J/m" +{ +0.0: :: "Positive or negative" +} 0.0 + + +# Parameters for Schwarzschild star. + +REAL starschwarz_m "Mass of Schwarzschild star" +{ +0.0: :: "Positive please" +} 1.0 + + +REAL starschwarz_r "Radius of Schwarzschild star" +{ +0.0: :: "Positive please" +} 1.0 + + +# Parameters for bowl metric. + +KEYWORD bowl_type "Type of bowl metric" +{ + gauss :: "Gaussian bowl" + fermi :: "Fermi bowl" +} "gauss" + +BOOLEAN bowl_evolve "Evolving bowl metric?" +{ +} "no" + +REAL bowl_a "Bowl strength" +{ +0.0: :: "Positive please" +} 0.5 + +REAL bowl_c "Center of deformation" +{ +0.0: :: "Positive please" +} 2.5 + +REAL bowl_s "Width of deformation" +{ +0.0: :: "Positive please" +} 2.5 + +REAL bowl_dx "Scale factor in x direction" +{ +0.0: :: "Positive please" +} 1.0 + +REAL bowl_dy "Scale factor in y direction" +{ +0.0: :: "Positive please" +} 1.0 + +REAL bowl_dz "Scale factor in z direction" +{ +0.0: :: "Positive please" +} 1.0 + +REAL bowl_t0 "Center of Fermi step in time" +{ +0.0: :: "Positive please" +} 1.0 + +REAL bowl_st "Width of Fermi step in time" +{ +0.0: :: "Positive please" +} 1.0 + + +# Parameters for Thorne's fake finary. + +KEYWORD fakebinary_atype "Thorne's binary type" +{ + "constant" :: + "quadrupole :: +} "constant" + +BOOLEAN fakebinary_retarded "Use retarded time?" +{ +} "no" + +REAL fakebinary_eps "Thorne's binary: fudge parameter" +{ +0.0: :: "Positive please" +} 1.e-16 + +REAL fakebinary_a0 "Thorne's binary: initial separation" +{ +0.0: :: "Positive please" +} 5.0 + +REAL fakebinary_omega0 "Thorne's binary: initial angular frequency" +{ +: :: "Positive or negative" +} 1.0 + +REAL fakebinary_m "Thorne's binary: mass" +{ +0.0: :: "Positive please" +} 0.1 + +REAL fakebinary_bround "Thorne's binary: smoothing for Newtonian potential" +{ +: :: "Positive or negative" +} 0.0 + + +# Parameters for multiBH. + +INT kt_nBH "number of black holes 0-4" +{ +0: :: "Positive please" +} 0 + +REAL kt_hubble "Hubble constant= pm sqrt{Lambda/3}" +{ +: :: "Positive or negative" +} 0.0 + +REAL m_bh1 "mass of black hole 1" +{ +0.0: :: "Positive please" +} 0.0 + +REAL m_bh2 "mass of black hole 2" +{ +0.0: :: "Positive please" +} 0.0 + +REAL m_bh3 "mass of black hole 3" +{ +0.0: :: "Positive please" +} 0.0 + +REAL m_bh4 "mass of black hole 4" +{ +0.0: :: "Positive please" +} 0.0 + +REAL co_bh1x "x coord of black hole 1" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh1y "y coord of black hole 1" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh1z "z coord of black hole 1" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh2x "x coord of black hole 2" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh2y "y coord of black hole 2" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh2z "z coord of black hole 2" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh3x "x coord of black hole 3" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh3y "y coord of black hole 3" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh3z "z coord of black hole 3" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh4x "x coord of black hole 4" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh4y "y coord of black hole 4" +{ +: :: "Positive or negative" +} 0.0 + +REAL co_bh4z "z coord of black hole 4" +{ +: :: "Positive or negative" +} 0.0 + + + diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..6e901c9 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,2 @@ +# Schedule definitions for thorn Exact +# $Header$ diff --git a/src/Comparison.c b/src/Comparison.c new file mode 100644 index 0000000..bc5554a --- /dev/null +++ b/src/Comparison.c @@ -0,0 +1,331 @@ +/* Please god comment me ... */ + +#include "cctk.h" + +static int comp_initialized = 0; +static int comp_active = -1; + +static int comp_metric = 0; +static int comp_curv = 0; +static int comp_gauge = 0; + +#define f_comparisonmetric FORTRAN_NAME(comparisonmetric_,COMPARISONMETRIC,comparisonmetric) +void FMODIFIER f_comparisonmetric(FORT_ARGS_PROTO, PROTO_FORT_FIELDS, + double *, + double *, + double *, + double *, + double *, + double *); +#define f_comparisoncurvature FORTRAN_NAME(comparisoncurvature_,COMPARISONCURVATURE,comparisoncurvature) +void FMODIFIER f_comparisoncurvature(FORT_ARGS_PROTO, PROTO_FORT_FIELDS, + double *, + double *, + double *, + double *, + double *, + double *); +#define f_comparisongauge FORTRAN_NAME(comparisongauge_,COMPARISONGAUGE,comparisongauge) +void FMODIFIER f_comparisongauge(FORT_ARGS_PROTO, PROTO_FORT_FIELDS, + double *, + double *, + double *, + double *); + +static int comptrips=1; + +void Comparison(pGH *in) { + void DoComparison(pGH *high, pGH *med, pGH *low, + int gfi, double *he, double *me, double *le); + + int i,j,k; + double *h[6]; + double *m[6]; + double *l[6]; + + int compevery = 1; + + pGH *high, *med, *low; + + /* Grab the approriate grid heirarchies */ + if (in->convlevel != 0) return; + + high = in; + med = GetGHbyLevel(1,in->level, in->mglevel); + low = GetGHbyLevel(2,in->level, in->mglevel); + + /* Only compare when grids are aligned. */ + if (med) compevery = compevery * 2; + if (low) compevery = compevery * 2; + + if (iteration % compevery != 0) return; + + /* See if I'm active */ + if (comp_active < 0) { + comp_active = Contains("comparison","yes"); + } + if (!comp_active) return; + + if (!comp_initialized) { + pGF *tGF; + comp_initialized =1; + for (i=0; i<high->ngridFuncs;i++) { + tGF = high->gridFuncs[i]; + if (Contains_bounded("compfields",tGF->name)) { + tGF->do_comparison = 1; + if (tGF->gfno >= GFI_GXX && + tGF->gfno <= GFI_GZZ) + comp_metric = 1; + + if (tGF->gfno >= GFI_HXX && + tGF->gfno <= GFI_HZZ) + comp_curv = 1; + + if (tGF->gfno == GFI_ALP || + tGF->gfno == GFI_BETAX || + tGF->gfno == GFI_BETAY || + tGF->gfno == GFI_BETAZ) + comp_gauge = 1; + } + } + } + if (!(comp_metric || comp_gauge || comp_curv)) { + comp_active = 0; + printf ("You must compare gauge, curvature, or metric"); + return; + } + + /* OK at this point we know we are going to do something, so + allocate the memory for the comparisons */ + for (i=0;i<6;i++) { + h[i] = (double *)malloc(high->npoints*sizeof(double)); + if (med) + m[i] = (double *)malloc( med->npoints*sizeof(double)); + else + m[i] = NULL; + if (low) + l[i] = (double *)malloc( low->npoints*sizeof(double)); + else + l[i] = NULL; + } + + EnableGFDataStorage(high, high->gridFuncs[high->ngridFuncs-1]); + if (med) + EnableGFDataStorage( med, med->gridFuncs[med->ngridFuncs-1]); + if (low) + EnableGFDataStorage( low, low->gridFuncs[low->ngridFuncs-1]); + + if (comp_metric) { + SetupFortranArrays(high); + f_comparisonmetric(FORT_ARGS(high), PASS_FORT_FIELDS(high), + h[0],h[1],h[2],h[3],h[4],h[5]); + + if (med) { + SetupFortranArrays(med); + f_comparisonmetric(FORT_ARGS(med), PASS_FORT_FIELDS(med), + m[0],m[1],m[2],m[3],m[4],m[5]); + } + + if (low) { + SetupFortranArrays(low); + f_comparisonmetric(FORT_ARGS(low), PASS_FORT_FIELDS(low), + l[0],l[1],l[2],l[3],l[4],l[5]); + } + + for (i=GFI_GXX; i<=GFI_GZZ;i++) { + if (high->gridFuncs[i]->do_comparison) { + DoComparison(high,med,low,i, + h[i-GFI_GXX],m[i-GFI_GXX],l[i-GFI_GXX]); + } + } + } + + if (comp_curv) { + SetupFortranArrays(high); + f_comparisoncurvature(FORT_ARGS(high), PASS_FORT_FIELDS(high), + h[0],h[1],h[2],h[3],h[4],h[5]); + + if (med) { + SetupFortranArrays(med); + f_comparisoncurvature(FORT_ARGS(med), PASS_FORT_FIELDS(med), + m[0],m[1],m[2],m[3],m[4],m[5]); + } + + if (low) { + SetupFortranArrays(low); + f_comparisoncurvature(FORT_ARGS(low), PASS_FORT_FIELDS(low), + l[0],l[1],l[2],l[3],l[4],l[5]); + } + + for (i=GFI_HXX; i<=GFI_HZZ;i++) { + if (high->gridFuncs[i]->do_comparison) { + DoComparison(high,med,low,i, + h[i-GFI_HXX],m[i-GFI_HXX],l[i-GFI_HXX]); + } + } + } + + if (comp_gauge) { + SetupFortranArrays(high); + f_comparisongauge(FORT_ARGS(high), PASS_FORT_FIELDS(high), + h[0],h[1],h[2],h[3]); + + if (med) { + SetupFortranArrays(med); + f_comparisongauge(FORT_ARGS(med), PASS_FORT_FIELDS(med), + m[0],m[1],m[2],m[3]); + } + + if (low) { + SetupFortranArrays(low); + f_comparisongauge(FORT_ARGS(low), PASS_FORT_FIELDS(low), + l[0],l[1],l[2],l[3]); + } + + if (high->gridFuncs[GFI_ALP]->do_comparison) { + DoComparison(high,med,low,GFI_ALP,h[0],m[0],l[0]); + } + + if (high->gridFuncs[GFI_BETAX]->do_comparison) { + DoComparison(high,med,low,GFI_BETAX,h[1],m[1],l[1]); + } + + if (high->gridFuncs[GFI_BETAY]->do_comparison) { + DoComparison(high,med,low,GFI_BETAY,h[2],m[2],l[2]); + } + + if (high->gridFuncs[GFI_BETAZ]->do_comparison) { + DoComparison(high,med,low,GFI_BETAZ,h[3],m[3],l[3]); + } + } + + for (i=0;i<6;i++) { + free(h[i]); if (m[i]) free(m[i]); if (l[i]) free(l[i]); + } + DisableGFDataStorage(high, high->gridFuncs[high->ngridFuncs-1]); + if (med) + DisableGFDataStorage( med, med->gridFuncs[med->ngridFuncs-1]); + if (low) + DisableGFDataStorage( low, low->gridFuncs[low->ngridFuncs-1]); + comptrips ++ ; +} + +void DoComparison(pGH *high, pGH *med, pGH *low, int gfi, + double *he, double *me, double *le) { + + pGF *ws, *cur; + double max[3], nm1[3], nm2[3]; + double sigtop, sigbot, sig3w, sig10, sig21; + + int i; + + printf ("Comparing %s\n",high->gridFuncs[gfi]->name); + + /* Output the analytic solution if we need it */ + cur = high->gridFuncs[gfi]; ws = high->gridFuncs[high->ngridFuncs-1]; + if (cur->do_1dio) { + ws->do_1dio = comptrips; + for (i=0;i<high->npoints;i++) + ws->data[i] = he[i]; + sprintf(ws->name,"%s_exact",cur->name); + IO_Write1D(high,ws); + ws->do_1dio = 0; + ws->lastio_it[1]--; + } + + /* Diffs go into workspace */ + for (i=0;i<high->npoints;i++) + ws->data[i] = cur->data[i] - he[i]; + + if (cur->do_1dio) { + ws->do_1dio = comptrips; + sprintf(ws->name,"%s_diff",cur->name); + IO_Write1D(high,ws); + ws->do_1dio = 0; + ws->lastio_it[1]--; + } + + max[0] = pGF_MaxVal(high,ws); + nm1[0] = pGF_Norm1 (high,ws); + nm2[0] = pGF_Norm2 (high,ws); + + if (med) { + cur = med->gridFuncs[gfi]; ws = med->gridFuncs[med->ngridFuncs-1]; + for (i=0;i<med->npoints;i++) + ws->data[i] = cur->data[i] - me[i]; + if (cur->do_1dio) { + ws->do_1dio = comptrips; + sprintf(ws->name,"%s_diff",cur->name); + IO_Write1D(med,ws); + ws->do_1dio = 0; + ws->lastio_it[1]--; + } + + max[1] = pGF_MaxVal(med,ws); + nm1[1] = pGF_Norm1 (med,ws); + nm2[1] = pGF_Norm2 (med,ws); + } + + if (low) { + cur = low->gridFuncs[gfi]; ws = low->gridFuncs[low->ngridFuncs-1]; + for (i=0;i<low->npoints;i++) + ws->data[i] = cur->data[i] - le[i]; + if (cur->do_1dio) { + ws->do_1dio = comptrips; + sprintf(ws->name,"%s_diff",cur->name); + IO_Write1D(low,ws); + ws->do_1dio = 0; + ws->lastio_it[1]--; + } + + max[2] = pGF_MaxVal(low,ws); + nm1[2] = pGF_Norm1 (low,ws); + nm2[2] = pGF_Norm2 (low,ws); + } else { + max[2] = 0.0; nm1[2] = 0.0; nm2[2] = 0.0; + } + + if (low) { + printf (" ------ : High Med Low\n"); + printf ("Max Diff: %lf %lf %lf\n", + max[0],max[1],max[2]); + printf ("NM1 Diff: %lf %lf %lf\n", + nm1[0],nm1[1],nm1[2]); + printf ("NM2 Diff: %lf %lf %lf\n", + nm2[0],nm2[1],nm2[2]); + } else { + if (med) { + printf (" ------ : High Med \n"); + printf ("Max Diff: %lf %lf\n", + max[0],max[1]); + printf ("NM1 Diff: %lf %lf\n", + nm1[0],nm1[1]); + printf ("NM2 Diff: %lf %lf\n", + nm2[0],nm2[1]); + } + } + + + if (low) { + sigbot = nm2[0] - nm2[1]; + if (sigbot == 0) sig3w = 0.0; + else sig3w = log(fabs((nm2[1]-nm2[2])/(nm2[0]-nm2[1])))/log(2.0); + + if (nm2[1] == 0) sig21 = 0.0; + else sig21 = log(fabs(nm2[2]/nm2[1]))/log(2.0); + } else { + sig3w = 0.0; sig21 = 0; + } + + if (nm2[0] == 0) sig10 = 0.0; + else sig10 = log(fabs(nm2[1]/nm2[0]))/log(2.0); + + if (low) + printf ("sigma: 3way %lf hm %lf ml %lf\n", + sig3w, sig10, sig21); + else + if (med) + printf ("sigma: hm %lf\n",sig10); + +} + diff --git a/src/ComparisonSolutions.F b/src/ComparisonSolutions.F new file mode 100644 index 0000000..e37ead0 --- /dev/null +++ b/src/ComparisonSolutions.F @@ -0,0 +1,183 @@ +#include "cctk.h" +#include "cctk_arguments.h" + + subroutine ComparisonMetric(CCTK_FARGUMENTS, + $ gxx_ex, gxy_ex, gxz_ex, + $ gyy_ex, gyz_ex, gzz_ex) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + CCTK_REAL gxx_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL gxy_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL gxz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL gyy_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL gyz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL gzz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + + integer i,j,k + +C Dummy arguments of subroutine boostrotdata: these are calculated +C (at a point) but not used. + + CCTK_REAL hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + +C Call boostrotdata pointwise. Most of what it calculates is +C thrown away, variables ending in ...junk. + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ gxx_ex(i,j,k), gyy_ex(i,j,k), gzz_ex(i,j,k), + $ gxy_ex(i,j,k), gyz_ex(i,j,k), gxz_ex(i,j,k), + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + end do + end do + end do + + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine ComparisonCurvature(CCTK_FARGUMENTS, + $ hxx_ex, hxy_ex, hxz_ex, + $ hyy_ex, hyz_ex, hzz_ex) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + CCTK_REAL hxx_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL hxy_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL hxz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL hyy_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL hyz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL hzz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + + integer i,j,k + + CCTK_REAL gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxx_ex(i,j,k), hyy_ex(i,j,k), hzz_ex(i,j,k), + $ hxy_ex(i,j,k), hyz_ex(i,j,k), hxz_ex(i,j,k), + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + end do + end do + end do + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +c Note the exact shift comes in even if the shift is not +c allocated (eg if shift is "none"). In this case just +c don't use it, since it won't be compared against acctk_lsh(2)thing. + + subroutine ComparisonGauge(CCTK_FARGUMENTS, + $ alp_ex, betax_ex, betay_ex, betaz_ex) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + CCTK_REAL alp_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL betax_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL betay_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + CCTK_REAL betaz_ex(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + + integer i,j,k + + CCTK_REAL gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ axjunk, ayjunk, azjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + + do k=1,cctk_lsh(3) + do j=1,cctk_lsh(2) + do i=1,cctk_lsh(1) + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alp_ex(i,j,k), axjunk, ayjunk, azjunk, + $ betax_ex(i,j,k), betay_ex(i,j,k), betaz_ex(i,j,k), + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + end do + end do + end do + + return + end diff --git a/src/boostrotmetric.F b/src/boostrotmetric.F new file mode 100644 index 0000000..6a775d3 --- /dev/null +++ b/src/boostrotmetric.F @@ -0,0 +1,165 @@ + +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine boostrotmetric( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + CCTK_REAL a, b, mu0, mu1, lam1, mu2, lam2, + $ lam3, mu4, lam4, mu5, lam5, num, div, f, + $ elam, emu0, delta, gfunc, tmp + CCTK_REAL h, d, numlim + + external gfunc + + logical firstcall + + data firstcall /.true./ + save firstcall, h, d, numlim + +C Get parameters of the exact solution. + + if (firstcall) then + + h = boostrotscale + d = boostrotstrength + + numlim = boostrotsafedistance + + firstcall = .false. + + end if + +C Intermediate quantities. + + a = x**2 + y**2 + b = z**2 - t**2 + + num = (0.5d0 * (a + b) - h)**2 + 2.d0 * h * a + +C Make sure we are not sitting on one of the two source wordlines, +C given by x = y = 0, z = +/- sqrt(h^2 + t^2) + + if (num / h**4 .le. numlim) stop 'too close to source wordline' + + div = 1.d0 / sqrt(num**3) + f = d**2 * ((0.25d0 * (a + b)**2 - h**2)**2 + $ - 0.5d0 * h**2 * a * b) / num**4 + + mu0 = - d * div * (0.5d0 * a**2 + h * a) + mu1 = - d * div * (0.5d0 * b + a - h) + lam1 = d * div * (0.5d0 * b - h) - a * f + mu2 = gfunc(b, mu1) + lam2 = gfunc(b, lam1) + + lam3 = d * div * (0.5d0 * b**2 - h * b) + lam4 = - d * div * (0.5d0 * a + h) - b * f + mu4 = - d * div * (0.5d0 * a + b + h) + mu5 = gfunc(a, - mu4) + lam5 = gfunc(a, lam4) + + elam = exp(lam3 + a * lam4) + emu0 = exp(mu0) + delta = exp(lam3) * (mu5 - lam5) + +C All nonvanishing metric coefficients (downstairs). + + gdxx = elam + y**2 * Delta + gdyy = elam + x**2 * Delta + gdxy = - x * y * Delta + gdzz = emu0 * (1.d0 + lam2 * z**2 - mu2 * t**2) + gdtz = - emu0 * z * t * (lam2 - mu2) + gdtt = - emu0 * (1.d0 + mu2 * z**2 - lam2 * t**2) + +C Others. + + gdzx = 0.d0 + gdyz = 0.d0 + gdtx = 0.d0 + gdty = 0.d0 + +C It is clear that the 3-metric is always spacelike in the xy plane. So +C we only need to check that gdzz is positive. + + if (gdzz .le. 0.d0) then + write(*,*) 'WARNING 3-metric not spacelike in boostrot at' + write(*,*) 't =', t, 'z =', z + write(*,*) 'x =', x, 'y =', y + pause + end if + +C Calculate inverse metric. That is not too difficult as it is +c in block-diagonal form. + + tmp = gdtt * gdzz - gdtz**2 + + if (tmp .eq. 0.d0) then + write(*,*) 'boostrot metric inverse failed in tz plane' + STOP + end if + + gutt = gdzz / tmp + guzz = gdtt / tmp + gutz = - gdtz / tmp + + tmp = gdxx * gdyy - gdxy**2 + + if (tmp .eq. 0.d0) then + write(*,*) 'boostrot metric inverse failed in xy plane' + STOP + end if + + guxx = gdyy / tmp + guyy = gdxx / tmp + guxy = - gdxy / tmp + + guzx = 0.d0 + guyz = 0.d0 + gutx = 0.d0 + guty = 0.d0 + + return + end + +ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + +C Calculates g = [exp (x f) - 1] / x as a power series for small x, +C so that the expression is regular at x = 0. + + function gfunc(x, f) + + implicit none + + integer n + + CCTK_REAL x, f, gfunc + CCTK_REAL sum, tmp + + if (abs(x*f) .ge. 1.d-6) then + gfunc = (exp(x*f) - 1.d0) / x + else + sum = 0.d0 + tmp = f + do n=1,10 + tmp = tmp / dble(n) + sum = sum + tmp + tmp = tmp * x * f + end do + gfunc = sum + end if + + return + end diff --git a/src/bowl.F b/src/bowl.F new file mode 100644 index 0000000..c9ba202 --- /dev/null +++ b/src/bowl.F @@ -0,0 +1,251 @@ +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine bowl( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + +c The metric given here is not a solution of +c Einstein's equations! It is nevertheless +c useful for tests since it has a particularly +c nice geometry. It is a static, spherically +c symmetric metric with no shift. +c +c In spherical coordinates, the metric has the +c form: +c +c 2 2 2 2 +c ds = dr + R(r) d Omega +c +c Clearly, r measures radial proper distance, and R(r) +c is the areal (Schwarzschild) radius. +c +c I choose a form of R(r) such that: +c +c R --> r r<<1, r>>1 +c +c So close in, and far away we have a flat metric. +c In the middle region, I take R to be smaller than +c r, but still larger than zero. This deficit in +c areal radius produces the geometry of a `bag of gold'. +c +c The size of the deviation from a flat geometry +c is controled by the parameter "bowl_a". +c If this parameter is 0, we are in flat space. +c The width of the curved region is controled by +c the paramter "bowl_s", and the place where the +c curvature becomes significant (the center of the +c deformation) is controled by "bowl_c". +c +c The specific form of the function R(r) is: +c +c R(r) = ( r - a f(r) ) +c +c and the form of thr function f(r) depends on the +c parameter bowl_type: +c 2 2 +c bowl_type = "gauss" f(r) = exp(-(r-c) / s ) +c +c bowl_type = "fermi" f(r) = 1 / ( 1 + exp(-s(r-c)) ) +c +c There are three extra paramters (bowl_dx,bowl_dy,bowl_dz) +c that scale the (x,y,z) axis respectively. Their default +c value is 1. These parameters are useful to hide the +c spherical symmetry of the metric. + + implicit none + + DECLARE_CCTK_PARAMETERS + +c Input. + CCTK_REAL x, y, z, t + +c Output. + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + + logical firstcall,evolve + + integer type + integer CCTK_Equals + + CCTK_REAL a,c,s + CCTK_REAL dx,dy,dz + CCTK_REAL r,r2,rr2 + CCTK_REAL xr,yr,zr,xr2,yr2,zr2 + CCTK_REAL fac,det + CCTK_REAL tfac,st,t0 + CCTK_REAL zero,one,two + + data firstcall /.true./ + save firstcall,evolve,type,a,c,s,dx,dy,dz,t0,st + +c Find {zero,one,two}. + + zero = 0.0d0 + one = 1.0d0 + two = 2.0d0 + +c Get parameters of the metric. + + if (firstcall) then + + a = bowl_a + c = bowl_c + s = bowl_s + + dx = bowl_dx + dy = bowl_dy + dz = bowl_dz + + if (dx.le.zero) then + call CCTK_WARN(0,"The parameter bowl_dx must be greater than 0") + end if + + if (dy.le.zero) then + call CCTK_WARN(0,"The parameter bowl_dy must be greater than 0") + end if + + if (dz.le.zero) then + call CCTK_WARN(0,"The parameter bowl_dz must be greater than 0") + end if + + if (CCTK_Equals(bowl_type,"gauss").ne.0) then + type = 1 + else if (CCTK_Equals(bowl_type,"fermi").ne.0) then + type = 2 + else + call CCTK_WARN(0,"Unknown type of bowl metric") + end if + + if (CCTK_Equals(bowl_evolve,"yes").ne.0) then + evolve = .true. + t0 = bowl_t0 + st = bowl_st + else + evolve = .false. + end if + + firstcall = .false. + + end if + +c Multiply the bowl strength `a' with a time evolution factor. +c The time evolution factor is taken to be a Fermi step centered +c in `t0' and with a width `st'. The size of this step is always +c 1 so that far in the past we will always have flat space, and +c far in the future we will have the static bowl. + + if (evolve) then + tfac = one/(one + dexp(-st*(t-t0))) + else + tfac = one + end if + + a = a*tfac + +c Find {r2,r}. + + r2 = (x/dx)**2 + (y/dy)**2 + (z/dz)**2 + r = dsqrt(r2) + +c Find the form function rr2 +c +c 2 2 2 +c rr2 = (r - a f) / r = (1 - a f / r) + + if (type.eq.1) then + +c Gaussian bowl: +c 2 2 2 +c rr2 = [ 1 - a exp(-(r-c) /s ) / r ] +c +c Notice that this really does not go to 1 at the +c origin. To fix this, I multiply the gaussian +c with the factor: +c +c fac = 1 - sech(4r) +c +c This goes smoothly to 0 at the origin, and climbs +c fast to a limiting value of 1 (at r=1 it is already +c equal to 0.96). + + fac = one - two/(dexp(4.0d0*r) + dexp(-4.0d0*r)) + rr2 = (one - a*fac*dexp(-((r-c)/s)**2)/r)**2 + + else if (type.eq.2) then + +c Fermi bowl: +c 2 +c rr2 = [ 1 - 1 / ( 1 + exp(-s(r-c)) ) / r ] +c +c Again, this doesn't really go to 1 at the origin, so +c I use the same trick as above. + + fac = one - two/(dexp(4.0d0*r) + dexp(-4.0d0*r)) + rr2 = (one - a*fac/(one + dexp(-s*(r-c)))/r)**2 + + end if + +c Give metric components. + + gdtt = - one + gdtx = zero + gdty = zero + gdtz = zero + + if (r.ne.0) then + + xr = (x/dx)/r + yr = (y/dy)/r + zr = (z/dz)/r + + xr2 = xr**2 + yr2 = yr**2 + zr2 = zr**2 + + gdxx = (xr2 + rr2*(yr2 + zr2))/dx**2 + gdyy = (yr2 + rr2*(xr2 + zr2))/dy**2 + gdzz = (zr2 + rr2*(xr2 + yr2))/dz**2 + + gdxy = xr*yr*(one - rr2)/(dx*dy) + gdyz = yr*zr*(one - rr2)/(dy*dz) + gdzx = xr*zr*(one - rr2)/(dx*dz) + + else + + gdxx = one + gdyy = one + gdzz = one + + gdxy = zero + gdyz = zero + gdzx = zero + + end if + +c Find inverse metric. + + gutt = - one + gutx = zero + guty = zero + gutz = zero + + det = gdxx*gdyy*gdzz + two*gdxy*gdzx*gdyz + . - gdxx*gdyz**2 - gdyy*gdzx**2 - gdzz*gdxy**2 + + guxx = (gdyy*gdzz - gdyz**2)/det + guyy = (gdxx*gdzz - gdzx**2)/det + guzz = (gdxx*gdyy - gdxy**2)/det + + guxy = (gdzx*gdyz - gdzz*gdxy)/det + guyz = (gdxy*gdzx - gdxx*gdyz)/det + guzx = (gdxy*gdyz - gdyy*gdzx)/det + + return + end diff --git a/src/exactblendbound.F b/src/exactblendbound.F new file mode 100644 index 0000000..bb7f1b7 --- /dev/null +++ b/src/exactblendbound.F @@ -0,0 +1,193 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" +#include "CactusEinstein/Einstein/src/Einstein.h" + + subroutine exactblendboundary(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + logical doKij, doGij, doLapse, doShift + + integer i,j,k + integer nx,ny,nz + integer CCTK_Equals + + CCTK_REAL router, rinner, frac, onemfrac + CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze + CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze + CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze + CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze + CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze + + CCTK_REAL alpe, axe, aye, aze + CCTK_REAL betaxe,betaye,betaze + CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze + + CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL vxe,vye,vze,sav + + CCTK_REAL :: dx,dy,dz,time + CCTK_REAL :: ierr,xmx,xmn,ymx,ymn,zmx,zmn,rmx + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + + time = cctk_time + +C Other parameters. + + doKij = (exblend_Ks.eq.1) + doGij = (exblend_gs.eq.1) + + doLapse = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(slicing,"exact").ne.0)) + doShift = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(shift,"exact").ne.0)) + + call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,"x") + call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,"y") + call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,"z") + + rmx = min(xmx,ymx,zmx) + + if (exblend_rout.lt.0) then + router = rmx - 2.0d0*dx + else + router = exblend_rout + endif + + if (exblend_width.lt.0) then + rinner = router + exblend_width*dx + else + rinner = router - exblend_width + endif + + do k=1,nz + do j=1,ny + do i=1,nx + +c We only do anything if r >= rinner so only evaluate exact data +c there. + + if (r(i,j,k) .ge. rinner) then + + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), time, + $ gxxe, gyye, gzze, gxye, gyze, gxze, + $ kxxe, kyye, kzze, kxye, kyze, kxze, + $ dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze, + $ dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze, + $ dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze, + $ alpe, axe, aye, aze, betaxe, betaye, betaze, + $ bxxe, bxye, bxze, byxe, + $ byye, byze, bzxe, bzye, bzze) + +c This sucks, but we want the exact vs so we can blend them also. + + det = -(gxze**2*gyye) + & + 2.d0*gxye*gxze*gyze + & - gxxe*gyze**2 + & - gxye**2*gzze + & + gxxe*gyye*gzze + + uxx=(-gyze**2 + gyye*gzze)/det + uxy=(gxze*gyze - gxye*gzze)/det + uyy=(-gxze**2 + gxxe*gzze)/det + uxz=(-gxze*gyye + gxye*gyze)/det + uyz=(gxye*gxze - gxxe*gyze)/det + uzz=(-gxye**2 + gxxe*gyye)/det + +c Outside of router we want to place exact data on our grid + + if (r(i,j,k) .gt. router) then + +c This is one of those things I will invariably screw up if I type +c it in so let the computer do it + +#define exassign(q) q(i,j,k) = q e +#define exassign_grp(p) \ + exassign(p xx) &&\ + exassign(p xy) &&\ + exassign(p xz) &&\ + exassign(p yy) &&\ + exassign(p yz) &&\ + exassign(p zz) + +c Note this plays on the nasty trick that fortran doesn't give a +c damn about spaces so gxx e is the same as gxxe for the parser... +c Grody but effective! + + if (doGij) then + exassign_grp(g) + endif + if (doKij) then + exassign_grp(k) + endif + + if (doLapse) then + exassign(alp) + endif + + if (doShift.and.(shift_state.ne.SHIFT_INACTIVE)) then + exassign(betax) + exassign(betay) + exassign(betaz) + endif + +c OK so we don't want to blend so use a goto to jump. + else + +c Evaluate the linear weighting fraction. Obvious... + + frac = (r(i,j,k) - rinner) / (router - rinner) + onemfrac = 1.0D0 - frac + +c Once again some c-preprocessor tricks based on the whole fortran +c space thing... + +#define INTPOINT(f,v) f(i,j,k) = frac * v + onemfrac * f(i,j,k) +#define intone(f) INTPOINT(f, f e) +#define int_grp(p) \ + intone(p xx) &&\ + intone(p xy) &&\ + intone(p xz) &&\ + intone(p yy) &&\ + intone(p yz) &&\ + intone(p zz) + + if (doGij) then + int_grp(g) + endif + if (doKij) then + int_grp(k) + endif + + if (doLapse) then + intone(alp) + endif + + if (doShift.and.(shift_state.ne.SHIFT_INACTIVE)) then + intone(betax) + intone(betay) + intone(betaz) + endif + + endif ! r > router else + endif ! r > rinner + + enddo + enddo + enddo + + return + end diff --git a/src/exactboundary.F b/src/exactboundary.F new file mode 100644 index 0000000..5e079ec --- /dev/null +++ b/src/exactboundary.F @@ -0,0 +1,123 @@ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine exactboundary(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL tplusone + CCTK_REAL + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ axjunk, ayjunk, azjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Set all initial data including dijk and vi on all points which +C are on the boundary of the domain if it really is the boundary +C of the complete grid. Treat all six sides of the grid cube this way. + +c Set t = time + dt. This is necessary here because by the time +c we reach this point the geometry has been evolved one time step +c but the variable `time' still hasn't been updated. + + tplusone = cctk_time + cctk_delta_time + +C Note we also always set the lapse and shift at the boundaries at +C time t+1. This is to provide boundary conditions for testing +C elliptic gauge conditions. If they are not used, they will be +C overwritten by exactgauge. + +#define EXACTDATAPOINT \ + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), tplusone, \ + gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), \ + gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), \ + kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), \ + kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), \ + dxgxxjunk, dxgyyjunk, dxgzzjunk, \ + dxgxyjunk, dxgyzjunk, dxgxzjunk, \ + dygxxjunk, dygyyjunk, dygzzjunk, \ + dygxyjunk, dygyzjunk, dygxzjunk, \ + dzgxxjunk, dzgyyjunk, dzgzzjunk, \ + dzgxyjunk, dzgyzjunk, dzgxzjunk, \ + alp(i,j,k), axjunk, ayjunk, azjunk, \ + betax(i,j,k), betay(i,j,k), betaz(i,j,k), \ + bxxjunk, bxyjunk, bxzjunk, \ + byxjunk, byyjunk, byzjunk, \ + bzxjunk, bzyjunk, bzzjunk) &&\ + + if (cctk_bbox(1) .eq. 1) then + i=1 + do j=1,ny + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(2) .eq. 1) then + i=nx + do j=1,ny + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(3) .eq. 1) then + j=1 + do i=1,nx + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(4) .eq. 1) then + j=ny + do i=1,nx + do k=1,nz + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(5) .eq. 1) then + k=1 + do j=1,ny + do i=1,nx + EXACTDATAPOINT + end do + end do + end if + + if (cctk_bbox(6) .eq. 1) then + k=nz + do j=1,ny + do i=1,nx + EXACTDATAPOINT + end do + end do + end if + + return + end + diff --git a/src/exactcartblendbound.F b/src/exactcartblendbound.F new file mode 100644 index 0000000..91efdcc --- /dev/null +++ b/src/exactcartblendbound.F @@ -0,0 +1,207 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" +#include "CactusEinstein/Einstein/src/Einstein.h" + + subroutine exactcartblendboundary(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + logical doKij, doGij, doDs, doVs, doLapse, doShift + + integer i,j,k + integer nx,ny,nz + integer ninterps + integer CCTK_Equals + + CCTK_REAL xblend, yblend, zblend + CCTK_REAL xmin, xmax, ymin, ymax, zmin, zmax + CCTK_REAL xfrac, yfrac, zfrac, onemxfrac, onemyfrac, onemzfrac + CCTK_REAL oonints + CCTK_REAL sfrac, onemsfrac + CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze + CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze + CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze + CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze + CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze + CCTK_REAL alpe, axe, aye, aze + CCTK_REAL betaxe,betaye,betaze + CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze + CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz + CCTK_REAL vxe,vye,vze,sav + + CCTK_REAL :: dx,dy,dz,time,ierr + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + + time = cctk_time + +C Other parameters. + + doKij = (exblend_Ks.eq.1) + doGij = (exblend_gs.eq.1) + + doLapse = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals("slicing","exact").ne.0)) + doShift = ((exblend_gauge.eq.1).and. + $ (CCTK_Equals(shift,"exact").ne.0)) + + if (exblend_width.lt.0) then + xblend = - exblend_width*dx + yblend = - exblend_width*dy + zblend = - exblend_width*dz + else + xblend = exblend_width + yblend = exblend_width + zblend = exblend_width + endif + + call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,"x") + call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,"y") + call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z") + + do k=1,nz + do j=1,ny + do i=1,nx + +c We only do anything if in the blending region + + if (x(i,j,k) .ge. xmax - xblend .or. + $ x(i,j,k) .le. xmin + xblend .or. + $ y(i,j,k) .ge. ymax - yblend .or. + $ y(i,j,k) .le. ymin + yblend .or. + $ z(i,j,k) .ge. zmax - zblend .or. + $ z(i,j,k) .le. zmin + zblend) then + + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), time, + $ gxxe, gyye, gzze, gxye, gyze, gxze, + $ kxxe, kyye, kzze, kxye, kyze, kxze, + $ dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze, + $ dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze, + $ dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze, + $ alpe, axe, aye, aze, betaxe, betaye, betaze, + $ bxxe, bxye, bxze, byxe, + $ byye, byze, bzxe, bzye, bzze) + +c This sucks, but we want the exact vs so we can blend them also. + + det = -(gxze**2*gyye) + & + 2.d0*gxye*gxze*gyze + & - gxxe*gyze**2 + & - gxye**2*gzze + & + gxxe*gyye*gzze + + uxx=(-gyze**2 + gyye*gzze)/det + uxy=(gxze*gyze - gxye*gzze)/det + uyy=(-gxze**2 + gxxe*gzze)/det + uxz=(-gxze*gyye + gxye*gyze)/det + uyz=(gxye*gxze - gxxe*gyze)/det + uzz=(-gxye**2 + gxxe*gyye)/det + +c OK so 6 blending cases. If frac = 1 we get all exact + + ninterps = 0 + + xfrac = 0.0D0 + onemxfrac = 0.0D0 + yfrac = 0.0D0 + onemyfrac = 0.0D0 + zfrac = 0.0D0 + onemzfrac = 0.0D0 + + if (x(i,j,k) .le. xmin + xblend) then + xfrac = 1.0D0 - (x(i,j,k)-xmin) / xblend + onemxfrac = 1.0D0 - xfrac + ninterps = ninterps + 1 + endif + + if (x(i,j,k) .ge. xmax - xblend) then + xfrac = 1.0D0 - (xmax - x(i,j,k)) / xblend + onemxfrac = 1.0D0 - xfrac + ninterps = ninterps + 1 + endif + + if (y(i,j,k) .le. ymin + yblend) then + yfrac = 1.0D0 - (y(i,j,k)-ymin) / yblend + onemyfrac = 1.0D0 - yfrac + ninterps = ninterps + 1 + endif + + if (y(i,j,k) .ge. ymax - yblend) then + yfrac = 1.0D0 - (ymax - y(i,j,k)) / yblend + onemyfrac = 1.0D0 - yfrac + ninterps = ninterps + 1 + endif + + if (z(i,j,k) .le. zmin + zblend) then + zfrac = 1.0D0 - (z(i,j,k)-zmin) / zblend + onemzfrac = 1.0D0 - zfrac + ninterps = ninterps + 1 + endif + + if (z(i,j,k) .ge. zmax - zblend) then + zfrac = 1.0D0 - (zmax - z(i,j,k)) / zblend + onemzfrac = 1.0D0 - zfrac + ninterps = ninterps + 1 + endif + + oonints = 1.0D0 / ninterps + + if (ninterps .eq. 0 .or. ninterps .gt. 3) then + print *,"NINTERPS error", ninterps + STOP + endif + + sfrac = (xfrac + yfrac + zfrac) * oonints + onemsfrac = 1.0D0 - sfrac + +c Once again some c-preprocessor tricks based on the whole fortran +c space thing... + +#define INTPOINT(f,v) f(i,j,k) = sfrac * v + onemsfrac * f(i,j,k) +#define intone(f) INTPOINT(f, f e) +#define int_grp(p) \ + intone(p xx) &&\ + intone(p xy) &&\ + intone(p xz) &&\ + intone(p yy) &&\ + intone(p yz) &&\ + intone(p zz) + + if (doGij) then + int_grp(g) + endif + + if (doKij) then + int_grp(k) + endif + + if (doLapse) then + intone(alp) + endif + + if (doShift.and.(shift_state.ne.SHIFT_INACTIVE)) then + intone(betax) + intone(betay) + intone(betaz) + endif + + endif ! r > rinner + + enddo + enddo + enddo + + return + end diff --git a/src/exactdata.F b/src/exactdata.F new file mode 100644 index 0000000..3210715 --- /dev/null +++ b/src/exactdata.F @@ -0,0 +1,194 @@ + +C This routine calculates BM initial data from a subroutine +C exactmetric that calculates the spacetime metric and its +C inverse. + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine exactdata(x, y, z, t, + $ gxx, gyy, gzz, gxy, gyz, gxz, + $ hxx, hyy, hzz, hxy, hyz, hxz, + $ dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz, + $ dygxx, dygyy, dygzz, dygxy, dygyz, dygxz, + $ dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz, + $ alp, ax, ay, az, betax, betay, betaz, + $ bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz) + + implicit none + + CCTK_REAL x, y, z, t, + $ gxx, gyy, gzz, gxy, gyz, gxz, + $ hxx, hyy, hzz, hxy, hyz, hxz, + $ dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz, + $ dygxx, dygyy, dygzz, dygxy, dygyz, dygxz, + $ dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz, + $ alp, ax, ay, az, betax, betay, betaz, + $ bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz + +C gxx is g_xx etc. +C hxx is K_xx etc. +C dxgyy is (/2) dg_yy / dx (sic!) +C alp is N, betax is N^x etc. +C bxy is (/2) dN^y / dx (sic and sic!) +C ax is dN / dx / N (sic!) + + CCTK_REAL eps, + $ gdtt, gdtx, gdty, gdtz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m + + parameter (eps=1.d-6) + +C Get the spacetime metric and its inverse at the base point. + + call exactmetric( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gxx, gyy, gzz, + $ gxy, gyz, gxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + +C Calculate lapse and shift from the upper metric. + + alp = 1.d0 / sqrt(- gutt) + betax = - gutx / gutt + betay = - guty / gutt + betaz = - gutz / gutt + +C In order to calculate the derivatives of the lapse and shift from +C the contravariant metric, use g^tt = -1/N^2 and g^i0 = N^i / N^2 +C Note that ax is dN/dx / N and that bxy is 1/2 dN^y / dx. + +C Calculate x-derivatives. + + call exactmetric( + $ x+eps, y, z, t, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call exactmetric( + $ x-eps, y, z, t, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + dxgxx = (gdxx_p - gdxx_m) / 4.d0 / eps + dxgyy = (gdyy_p - gdyy_m) / 4.d0 / eps + dxgzz = (gdzz_p - gdzz_m) / 4.d0 / eps + dxgxy = (gdxy_p - gdxy_m) / 4.d0 / eps + dxgyz = (gdyz_p - gdyz_m) / 4.d0 / eps + dxgxz = (gdxz_p - gdxz_m) / 4.d0 / eps + ax = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt + bxx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps + bxy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps + bxz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps + +C Calculate y-derivatives. + + call exactmetric( + $ x, y+eps, z, t, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call exactmetric( + $ x, y-eps, z, t, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + dygxx = (gdxx_p - gdxx_m) / 4.d0 / eps + dygyy = (gdyy_p - gdyy_m) / 4.d0 / eps + dygzz = (gdzz_p - gdzz_m) / 4.d0 / eps + dygxy = (gdxy_p - gdxy_m) / 4.d0 / eps + dygyz = (gdyz_p - gdyz_m) / 4.d0 / eps + dygxz = (gdxz_p - gdxz_m) / 4.d0 / eps + ay = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt + byx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps + byy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps + byz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps + +C Calculate z-derivatives. + + call exactmetric( + $ x, y, z+eps, t, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call exactmetric( + $ x, y, z-eps, t, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + dzgxx = (gdxx_p - gdxx_m) / 4.d0 / eps + dzgyy = (gdyy_p - gdyy_m) / 4.d0 / eps + dzgzz = (gdzz_p - gdzz_m) / 4.d0 / eps + dzgxy = (gdxy_p - gdxy_m) / 4.d0 / eps + dzgyz = (gdyz_p - gdyz_m) / 4.d0 / eps + dzgxz = (gdxz_p - gdxz_m) / 4.d0 / eps + az = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt + bzx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps + bzy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps + bzz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps + +C Calculate t-derivatives, and extrinsic curvature. + + call exactmetric( + $ x, y, z, t+eps, + $ gdtt_p, gdtx_p, gdty_p, gdtz_p, + $ gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p, + $ gutt_p, gutx_p, guty_p, gutz_p, + $ guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p) + call exactmetric( + $ x, y, z, t-eps, + $ gdtt_m, gdtx_m, gdty_m, gdtz_m, + $ gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m, + $ gutt_m, gutx_m, guty_m, gutz_m, + $ guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m) + + hxx = - (gdxx_p - gdxx_m) / 4.d0 / eps / alp + $ + (dxgxx * betax + dygxx * betay + dzgxx * betaz + $ + 2.d0 * (bxx * gxx + bxy * gxy + bxz * gxz)) / alp + + hyy = - (gdyy_p - gdyy_m) / 4.d0 / eps / alp + $ + (dxgyy * betax + dygyy * betay + dzgyy * betaz + $ + 2.d0 * (byx * gxy + byy * gyy + byz * gyz)) / alp + + hzz = - (gdzz_p - gdzz_m) / 4.d0 / eps / alp + $ + (dxgzz * betax + dygzz * betay + dzgzz * betaz + $ + 2.d0 * (bzx * gxz + bzy * gyz + bzz * gzz)) / alp + + hxy = - (gdxy_p - gdxy_m) / 4.d0 / eps / alp + $ + (dxgxy * betax + dygxy * betay + dzgxy * betaz + $ + bxx * gxy + bxy * gyy + bxz * gyz + $ + byx * gxx + byy * gxy + byz * gxz) / alp + + hyz = - (gdyz_p - gdyz_m) / 4.d0 / eps / alp + $ + (dxgyz * betax + dygyz * betay + dzgyz * betaz + $ + byx * gxz + byy * gyz + byz * gzz + $ + bzx * gxy + bzy * gyy + bzz * gyz) / alp + + hxz = - (gdxz_p - gdxz_m) / 4.d0 / eps / alp + $ + (dxgxz * betax + dygxz * betay + dzgxz * betaz + $ + bxx * gxz + bxy * gyz + bxz * gzz + $ + bzx * gxx + bzy * gxy + bzz * gxz) / alp + + + return + end + + diff --git a/src/exactgauge.F b/src/exactgauge.F new file mode 100644 index 0000000..eae23f1 --- /dev/null +++ b/src/exactgauge.F @@ -0,0 +1,135 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + + subroutine exactgauge(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + + integer i,j,k + integer nx,ny,nz + integer CCTK_Equals + + CCTK_REAL tplushalf + CCTK_REAL gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Set t = time + dt/2. This is for consistency with a second +C order numerical scheme. + + tplushalf = cctk_time + 0.5D0*cctk_delta_time + +C Set both lapse and shift. + + if ((CCTK_Equals(slicing,"exact").ne.0).and. + . (CCTK_Equals(shift,"exact").ne.0) ) then + + do k=1,nz + do j=1,ny + do i=1,nx + + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alp(i,j,k), axjunk, ayjunk, azjunk, + $ betax(i,j,k), betay(i,j,k), betaz(i,j,k), + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + +C Set lapse only. + + else if (CCTK_Equals(slicing,"exact").ne.0) then + + do k=1,nz + do j=1,ny + do i=1,nx + + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alp(i,j,k), axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + +C Set shift only. + + else if (CCTK_Equals(shift,"exact").ne.0) then + + do k=1,nz + do j=1,ny + do i=1,nx + + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), tplushalf, + $ gxxjunk, gyyjunk, gzzjunk, + $ gxyjunk, gyzjunk, gxzjunk, + $ hxxjunk, hyyjunk, hzzjunk, + $ hxyjunk, hyzjunk, hxzjunk, + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betax(i,j,k), betay(i,j,k), betaz(i,j,k), + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + + end if + + return + end diff --git a/src/exactinitialize.F b/src/exactinitialize.F new file mode 100644 index 0000000..e345775 --- /dev/null +++ b/src/exactinitialize.F @@ -0,0 +1,71 @@ +C Wrapper for boostrotdata. Calls it and vectorini. +C Sets Cauchy data, lapse and shift, and what else is needed +C in the Bona-Masso formalism, at an initial time. + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" +#include "CactusEinstein/Einstein/src/Einstein.h" + + subroutine exactinitialize(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + integer i,j,k + integer nx,ny,nz + + CCTK_REAL alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk + CCTK_REAL + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk + +C Note I assume time has been initialized to physical time. +C Set data pointwise. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + do k=1,nz + do j=1,ny + do i=1,nx + + call exactdata(x(i,j,k), y(i,j,k), z(i,j,k), cctk_time, + $ gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), + $ gxy(i,j,k), gyz(i,j,k), gxz(i,j,k), + $ kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), + $ kxy(i,j,k), kyz(i,j,k), kxz(i,j,k), + $ dxgxxjunk, dxgyyjunk, dxgzzjunk, + $ dxgxyjunk, dxgyzjunk, dxgxzjunk, + $ dygxxjunk, dygyyjunk, dygzzjunk, + $ dygxyjunk, dygyzjunk, dygxzjunk, + $ dzgxxjunk, dzgyyjunk, dzgzzjunk, + $ dzgxyjunk, dzgyzjunk, dzgxzjunk, + $ alpjunk, axjunk, ayjunk, azjunk, + $ betaxjunk, betayjunk, betazjunk, + $ bxxjunk, bxyjunk, bxzjunk, + $ byxjunk, byyjunk, byzjunk, + $ bzxjunk, bzyjunk, bzzjunk) + + end do + end do + end do + +C Tell the code there is no need to treat the conformal factor +C as a separate field. That is, we have set the physical metric here. + + conformal_state = NOCONFORMAL_METRIC + psi = 1.d0 + + return + end diff --git a/src/exactmetric.F b/src/exactmetric.F new file mode 100644 index 0000000..48a604f --- /dev/null +++ b/src/exactmetric.F @@ -0,0 +1,136 @@ +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +C This is just a wrapper for different routines in place of exactmetric. + + subroutine exactmetric( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz + + logical firstcall + + integer CCTK_Equals + +C Call the corresponding routine. + + if (CCTK_Equals(exactmodel,"boostrot").eq.1) then + + call boostrotmetric( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"finkelstein").eq.1) then + + call finkelstein( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"kerrschild").eq.1) then + + call kerrschild( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"flatschwarz").eq.1) then + + call flatschwarz( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"starschwarz").eq.1) then + + call starschwarz( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"novikov").eq.1) then + + call novikov(x,y,z,t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"bowl").eq.1) then + + call bowl( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"minkowski").eq.1) then + + call minkowski( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"fakebinary").eq.1) then + + call fakebinary( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + elseif (CCTK_Equals(exactmodel,"multiBH").eq.1) then + + call multi_BHs( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, + $ gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + else + + call CCTK_WARN(0,"Unknown value of parameter exactmodel") + + end if + + return + end diff --git a/src/fakebinary.F b/src/fakebinary.F new file mode 100644 index 0000000..2db00a6 --- /dev/null +++ b/src/fakebinary.F @@ -0,0 +1,173 @@ +C fakebinary.F +C Bernd Bruegmann, 6/98 +C +C Compute Thorne's four metric that, although not a solution to the +C Einstein equations, has several characteristic features of a binary +C star system. See gr-qc/98...... + +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine fakebinary( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz) + + implicit none + + DECLARE_CCTK_PARAMETERS + +C input + CCTK_REAL x, y, z, t + +C output + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdxz, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guxz + +C static + + logical firstcall + + integer CCTK_Equals + + CCTK_REAL eps, m, a0, omega0, bround, atype, aretarded + + data firstcall /.true./ + save firstcall, eps, m, a0, omega0, bround, atype, aretarded + +C temps + CCTK_REAL a, omega, tau, f + CCTK_REAL c, c0, c1, c2, c3 + CCTK_REAL rho, r, sinp, cosp, phi, sint, cost, tx, ty, tz, px, py, pz + CCTK_REAL a2, b2, bx, by, bz, detgd + +C get parameters of the exact solution. + + if (firstcall) then + + firstcall = .false. + + eps = fakebinary_eps + m = fakebinary_m + a0 = fakebinary_a0 + + omega0 = fakebinary_omega0 + bround = fakebinary_bround + + bround = max(bround, eps) + + if (CCTK_Equals(fakebinary_atype,"constant").ne.0) then + atype = 0.d0 + elseif (CCTK_Equals(fakebinary_atype,"quadrupole").ne.0) then + atype = 1.d0 + else + call CCTK_Warn(0,"Unknown value of parameter fakebinary_atype") + endif + + if (fakebinary_retarded.ne.0) then + aretarded = 1.d0 + else + aretarded = 0.d0 + endif + + end if + +C spherical coordinates + + rho = max(sqrt(x**2 + y**2), eps) + r = sqrt(rho**2 + z**2) + sinp = y / rho + cosp = x / rho + phi = acos(cosp) + sint = rho / r + cost = z / r + tx = cost*cosp + ty = cost*sinp + tz = sint + px = - sinp + py = cosp + pz = 0 + +C distance function a(T-R) + + tau = 5.d0/128.d0 * a0**4 / m**3 + a = a0 * (1.d0 - atype * 4.d0*(t - aretarded*r)/tau)**(0.25d0) + +C orbital frequency omega(T-R) + + omega = 0.5d0*(m/a**3)**2 + +C 1/r type potential f + + c = y**2 + z**2 + bround**2; + f = ((x-a)**2 + c)**(-0.5d0) + ((x+a)**2 + c)**(-0.5d0) + +C the three metric, tt part + + c3 = 2.d0*(phi + omega*r) + c0 = - 4.d0 * m * a**2 * omega**3 * (omega*r)**4 + . / (1 + (omega*r)**2)**(2.5d0) + c1 = (1 + cost**2) * cos(c3) * c0 + c2 = - 2.d0 * cost * sin(c3) * c0 + gdxx = c1 * (tx*tx - px*px) + c2 * (tx*px + px*tx) + gdxy = c1 * (tx*ty - px*py) + c2 * (tx*py + px*ty) + gdxz = c1 * (tx*tz - px*pz) + c2 * (tx*pz + px*tz) + gdyy = c1 * (ty*ty - py*py) + c2 * (ty*py + py*ty) + gdyz = c1 * (ty*tz - py*pz) + c2 * (ty*pz + py*tz) + gdzz = c1 * (tz*tz - pz*pz) + c2 * (tz*pz + pz*tz) + +C the three metric, add conformally flat part + + c = (1.d0 + m * f)**2 + gdxx = gdxx + c + gdyy = gdyy + c + gdzz = gdzz + c + +C the shift vector and covector + + c = (1.d0 - 2*m*a**2/(r**2+a**2) * f) * omega * rho + bx = c * px + by = c * py + bz = c * pz + gdtx = gdxx*bx + gdxy*by + gdxz*bz + gdty = gdxy*bx + gdyy*by + gdyz*bz + gdtz = gdxz*bx + gdyz*by + gdzz*bz + b2 = gdtx*bx + gdty*by + gdtz*bz + +C lapse squard and time-time component of the four metric + + a2 = (1.d0 - m * f)**2 + gdtt = b2 - a2 + +C done with metric, find its inverse +C inverse three metric + + detgd = -(gdxz**2*gdyy) + 2*gdxy*gdxz*gdyz - gdxx*gdyz**2 + . - gdxy**2*gdzz - gdxx*gdyy*gdzz + guxx = (-gdyz**2 + gdyy*gdzz) / detgd + guxy = (gdxz*gdyz - gdxy*gdzz) / detgd + guxz = (-(gdxz*gdyy) + gdxy*gdyz) / detgd + guyy = (-gdxz**2 + gdxx*gdzz) / detgd + guyz = (gdxy*gdxz - gdxx*gdyz) / detgd + guzz = (-gdxy**2 + gdxx*gdyy) / detgd + +C inverse four metric + + gutt = - 1.d0/a2 + gutx = bx/a2 + guty = by/a2 + gutz = bz/a2 + guxx = guxx - bx*bx/a2 + guxy = guxy - bx*by/a2 + guxz = guxz - bx*bz/a2 + guyy = guyy - by*by/a2 + guyz = guyz - by*bz/a2 + guzz = guzz - bz*bz/a2 + +C done! + return + end diff --git a/src/finkelstein.F b/src/finkelstein.F new file mode 100644 index 0000000..1e96fe9 --- /dev/null +++ b/src/finkelstein.F @@ -0,0 +1,66 @@ +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine finkelstein( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL x, y, z, t + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + + logical firstcall + + CCTK_REAL eps, m + + data firstcall /.true./ + save firstcall, eps, m + + CCTK_REAL r + +C Get parameters of the exact solution. + + if (firstcall) then + + eps = kerrschild_eps + m = kerrschild_m + + firstcall = .false. + + end if + + r = max(sqrt(x**2 + y**2 + z**2), eps) + + gdtt = - (1.d0 - 2.d0 * m / r) + gdtx = 2.d0 * m * x / r**2 + gdty = 2.d0 * m * y / r**2 + gdtz = 2.d0 * m * z / r**2 + gdxx = 1.d0 + 2.d0 * m * x**2 / r**3 + gdyy = 1.d0 + 2.d0 * m * y**2 / r**3 + gdzz = 1.d0 + 2.d0 * m * z**2 / r**3 + gdxy = 2.d0 * m * x * y / r**3 + gdyz = 2.d0 * m * y * z / r**3 + gdzx = 2.d0 * m * z * x / r**3 + + gutt = - (1.d0 + 2.d0 * m / r) + gutx = 2.d0 * m * x / r**2 + guty = 2.d0 * m * y / r**2 + gutz = 2.d0 * m * z / r**2 + guxx = 1.d0 - 2.d0 * m * x**2 / r**3 + guyy = 1.d0 - 2.d0 * m * y**2 / r**3 + guzz = 1.d0 - 2.d0 * m * z**2 / r**3 + guxy = - 2.d0 * m * x * y / r**3 + guyz = - 2.d0 * m * y * z / r**3 + guzx = - 2.d0 * m * z * x / r**3 + + return + end diff --git a/src/flatschwarz.F b/src/flatschwarz.F new file mode 100644 index 0000000..9566517 --- /dev/null +++ b/src/flatschwarz.F @@ -0,0 +1,70 @@ +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine flatschwarz( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL x, y, z, t + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + + logical firstcall + + CCTK_REAL eps, m + + data firstcall /.true./ + save firstcall, eps, m + + CCTK_REAL r, bx, by, bz, b2 + +C Get parameters of the exact solution. + + if (firstcall) then + + eps = kerrschild_eps + m = kerrschild_m + + firstcall = .false. + + end if + + r = max(sqrt(x**2 + y**2 + z**2), eps) + bx = sqrt(2.d0 * m / r) * x / r + by = sqrt(2.d0 * m / r) * y / r + bz = sqrt(2.d0 * m / r) * z / r + b2 = 2.d0 * m / r + + gdtt = - 1.d0 + b2 + gdtx = bx + gdty = by + gdtz = bz + gdxx = 1.d0 + gdyy = 1.d0 + gdzz = 1.d0 + gdxy = 0.d0 + gdyz = 0.d0 + gdzx = 0.d0 + + gutt = - 1.d0 + gutx = bx + guty = by + gutz = bz + guxx = 1.d0 - bx**2 + guyy = 1.d0 - by**2 + guzz = 1.d0 - bz**2 + guxy = - bx * by + guyz = - by * bz + guzx = - bz * bx + + return + end diff --git a/src/include/pGF_ComparisonExtensions.h b/src/include/pGF_ComparisonExtensions.h new file mode 100644 index 0000000..0cca003 --- /dev/null +++ b/src/include/pGF_ComparisonExtensions.h @@ -0,0 +1 @@ + int do_comparison; diff --git a/src/include/slice_normal.h b/src/include/slice_normal.h new file mode 100644 index 0000000..4b41d85 --- /dev/null +++ b/src/include/slice_normal.h @@ -0,0 +1,103 @@ +C Fill in remaining components of metric. + do m=1,4 + do n=1,m-1 + gu(m,n) = gu(n,m) + gd(m,n) = gd(n,m) + end do + end do + +C Calculate normal vector with respect to the slice, with indices down, +C and not yet normalized. +C n_A = epsilon_ABCD X^A_,1 X^B_,2 X^C_,3. +C Easiest to hardwired this with an editor. + nd(1) = s1d(2,1)*s1d(3,2)*s1d(4,3) + $ + s1d(3,1)*s1d(4,2)*s1d(2,3) + $ + s1d(4,1)*s1d(2,2)*s1d(3,3) + $ - s1d(2,1)*s1d(4,2)*s1d(3,3) + $ - s1d(4,1)*s1d(3,2)*s1d(2,3) + $ - s1d(3,1)*s1d(2,2)*s1d(4,3) + + nd(2) = - s1d(3,1)*s1d(4,2)*s1d(1,3) + $ - s1d(4,1)*s1d(1,2)*s1d(3,3) + $ -s1d(1,1)*s1d(3,2)*s1d(4,3) + $ + s1d(3,1)*s1d(1,2)*s1d(4,3) + $ + s1d(1,1)*s1d(4,2)*s1d(3,3) + $ + s1d(4,1)*s1d(3,2)*s1d(1,3) + + nd(3) = s1d(4,1)*s1d(1,2)*s1d(2,3) + $ + s1d(1,1)*s1d(2,2)*s1d(4,3) + $ + s1d(2,1)*s1d(4,2)*s1d(1,3) + $ - s1d(4,1)*s1d(2,2)*s1d(1,3) + $ - s1d(2,1)*s1d(1,2)*s1d(4,3) + $ - s1d(1,1)*s1d(4,2)*s1d(2,3) + + nd(4) = s1d(1,1)*s1d(2,2)*s1d(3,3) + $ + s1d(2,1)*s1d(3,2)*s1d(1,3) + $ + s1d(3,1)*s1d(1,2)*s1d(2,3) + $ - s1d(1,1)*s1d(3,2)*s1d(2,3) + $ - s1d(3,1)*s1d(2,2)*s1d(1,3) + $ - s1d(2,1)*s1d(1,2)*s1d(3,3) + +C Normalize normal vector and raise its index, both with the exact +C solution metric. +C - N^2 = g^AB n_A n_B). + norm = 0.d0 + do m=1,4 + do n=1,4 + norm = norm + nd(m) * nd(n) * gu(m,n) + end do + end do +C Debugging test. + if (norm .ge. 0.d0) then + write(6,*) 'slice normal no longer timelike' + write(6,*) 'grid point', i, j, k + write(6,*) 'x^i', x(i,j,k), y(i,j,k), z(i,j,k) + write(6,*) 'x^A', slicex(i,j,k), slicey(i,j,k), + $ slicez(i,j,k), slicet(i,j,k) + write(6,*) 'n_A', nd(1), nd(2), nd(3), nd(4) + write(6,*) 'n_A n^A', norm + write(6,*) 'dx^A/dx^i' + do m=1,4 + do n=1,3 + write(6,*) m, n, s1d(m,n) + end do + end do + write(6,*) 'g^AB' + do m=1,4 + do n=1,4 + write(6,*) m, n, gu(m,n) + end do + end do + STOP + end if +C And now n^A = - 1/N g^AB n_B. (Sign so that it is future-pointing.) + do m=1,4 + nu(m) = 0.d0 + do n=1,4 + nu(m) = nu(m) + gu(m,n) * nd(n) + end do + nu(m) = - nu(m) / sqrt(-norm) + end do + +C dX^A/dt = alpha n^A + beta^i X^A_,i. Store as a GF. + slicetmp2x(i,j,k) = alp(i,j,k) * nu(1) + $ + betax(i,j,k) * s1d(1,1) + $ + betay(i,j,k) * s1d(1,2) + $ + betaz(i,j,k) * s1d(1,3) + + slicetmp2y(i,j,k) = alp(i,j,k) * nu(2) + $ + betax(i,j,k) * s1d(2,1) + $ + betay(i,j,k) * s1d(2,2) + $ + betaz(i,j,k) * s1d(2,3) + + slicetmp2z(i,j,k) = alp(i,j,k) * nu(3) + $ + betax(i,j,k) * s1d(3,1) + $ + betay(i,j,k) * s1d(3,2) + $ + betaz(i,j,k) * s1d(3,3) + + slicetmp2t(i,j,k) = alp(i,j,k) * nu(4) + $ + betax(i,j,k) * s1d(4,1) + $ + betay(i,j,k) * s1d(4,2) + $ + betaz(i,j,k) * s1d(4,3) + + diff --git a/src/kerrschild.F b/src/kerrschild.F new file mode 100644 index 0000000..2c37b37 --- /dev/null +++ b/src/kerrschild.F @@ -0,0 +1,127 @@ +#include "cctk.h" +#include "cctk_parameters.h" + +C Kerr-Schild form of boosted rotating black hole. +C Program g_ab = eta_ab + H l_a l_b, g^ab = eta^ab - H l^a l^b. +C Here eta_ab is Minkowski in Cartesian coordinates, H is a scalar, +C and l is a null vector. + + subroutine kerrschild( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL x, y, z, t + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + + CCTK_REAL gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, + $ lt0, lx0, ly0, lz0, hh, lt, lx, ly, lz + + logical firstcall + CCTK_REAL boostv, eps, m, a + data firstcall /.true./ + save firstcall, boostv, eps, m, a + + +C Get parameters of the exact solution. + + if (firstcall) then + + boostv = kerrschild_boostv + eps = kerrschild_eps + + m = kerrschild_m + a = kerrschild_a + + firstcall = .false. + + end if + +C Boost factor. + + gamma = 1.d0 / sqrt(1.d0 - boostv**2) + +C Lorentz transform t,x,y,z -> t0,x0,y0,z0. +C t0 is never used, but is here for illustration, and we introduce +C x0 and y0 also only for clarity. +C Note that z0 = 0 means z = vt for the BH. + + t0 = gamma * (t - boostv * z) + z0 = gamma * (z - boostv * t) + x0 = x + y0 = y + +C Coordinate distance to center of black hole. Note it moves! + + rho02 = x0**2 + y0**2 + z0**2 + +C Spherical auxiliary coordinate r and angle theta in BH rest frame. + + r02 = 0.5d0 * (rho02 - a**2) + $ + sqrt(0.25d0 * (rho02 - a**2)**2 + a**2 * z0**2) + if (r02 .lt. eps) then + write(6,*) 'r02 too small in kerrschild' + write(6,*) x0, y0, z0 + write(6,*) rho02, a**2, r02 + end if + r02 = max(r02,eps) + r0 = sqrt(r02) + costheta0 = z0 / r0 + +C Coefficient H. Note this transforms as a scalar, so does not carry +C the suffix 0. + hh = m * r0 / (r0**2 + a**2 * costheta0**2) + +C Components of l_a in rest frame. Note indices down. + lt0 = 1.d0 + lx0 = (r0 * x0 + a * y0) / (r0**2 + a**2) + ly0 = (r0 * y0 - a * x0) / (r0**2 + a**2) + lz0 = z0 / r0 + +C Now boost it to coordinates x, y, z, t. +C This is the reverse Lorentz transformation, but applied +C to a one-form, so the sign of boostv is the same as the forward +C Lorentz transformation applied to the coordinates. + + lt = gamma * (lt0 - boostv * lz0) + lz = gamma * (lz0 - boostv * lt0) + lx = lx0 + ly = ly0 + +C Down metric. g_ab = flat_ab + H l_a l_b + + gdtt = - 1.d0 + 2.d0 * hh * lt * lt + gdtx = 2.d0 * hh * lt * lx + gdty = 2.d0 * hh * lt * ly + gdtz = 2.d0 * hh * lt * lz + gdxx = 1.d0 + 2.d0 * hh * lx * lx + gdyy = 1.d0 + 2.d0 * hh * ly * ly + gdzz = 1.d0 + 2.d0 * hh * lz * lz + gdxy = 2.d0 * hh * lx * ly + gdyz = 2.d0 * hh * ly * lz + gdzx = 2.d0 * hh * lz * lx + +C Up metric. g^ab = flat^ab - H l^a l^b. +C Notice that g^ab = g_ab and l^i = l_i and l^0 = - l_0 in flat spacetime. + gutt = - 1.d0 - 2.d0 * hh * lt * lt + gutx = 2.d0 * hh * lt * lx + guty = 2.d0 * hh * lt * ly + gutz = 2.d0 * hh * lt * lz + guxx = 1.d0 - 2.d0 * hh * lx * lx + guyy = 1.d0 - 2.d0 * hh * ly * ly + guzz = 1.d0 - 2.d0 * hh * lz * lz + guxy = - 2.d0 * hh * lx * ly + guyz = - 2.d0 * hh * ly * lz + guzx = - 2.d0 * hh * lz * lx + + return + end diff --git a/src/linextraponebound.F b/src/linextraponebound.F new file mode 100644 index 0000000..09a2ea7 --- /dev/null +++ b/src/linextraponebound.F @@ -0,0 +1,140 @@ +#include "cctk.h" +#include "cctk_arguments.h" + + subroutine linextraponebound(CCTK_FARGUMENTS,var) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + integer nx,ny,nz + + CCTK_REAL var(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)) + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +C Linear extrapolation from the interiors to the boundaries. +C Does not support octant or quadrant. + +C 6 faces. + + if (cctk_bbox(1).eq.1 .and. nx.ge.4 ) then + var(1,:,:) = 2.d0 * var(2,:,:) - var(3,:,:) + endif + if (cctk_bbox(3).eq.1 .and. ny.ge.4 ) then + var(:,1,:) = 2.d0 * var(:,2,:) - var(:,3,:) + endif + if (cctk_bbox(5).eq.1 .and. nz.ge.4 ) then + var(:,:,1) = 2.d0 * var(:,:,2) - var(:,:,3) + endif + if (cctk_bbox(2).eq.1 .and. nx.ge.4 ) then + var(nx,:,:) = 2.d0 * var(nx-1,:,:) - var(nx-2,:,:) + endif + if (cctk_bbox(4).eq.1 .and. ny.ge.4 ) then + var(:,ny,:) = 2.d0 * var(:,ny-1,:) - var(:,ny-2,:) + endif + if (cctk_bbox(6).eq.1 .and. nz.ge.4 ) then + var(:,:,nz) = 2.d0 * var(:,:,nz-1) - var(:,:,nz-2) + endif + +C 12 edges. +C 4 round face x=min. + if ( cctk_bbox(1).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(3).eq.1 .and. ny.ge.4 ) then + var(1,1,:) = 2.d0 * var(2,2,:) - var(3,3,:) + end if + + if ( cctk_bbox(1).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(4).eq.1 .and. ny.ge.4 ) then + var(1,ny,:) = 2.d0 * var(2,ny-1,:) - var(2,ny-2,:) + end if + + if ( cctk_bbox(1).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then + var(1,:,1) = 2.d0 * var(2,:,2) - var(3,:,3) + end if + + if ( cctk_bbox(1).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then + var(1,:,nz) = 2.d0 * var(2,:,nz-1) - var(3,:,nz-2) + end if + +C 4 around face x=max. + if ( cctk_bbox(2).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(3).eq.1 .and. ny.ge.4 ) then + var(nx,1,:) = 2.d0 * var(nx-1,2,:) - var(nx-2,3,:) + end if + + if ( cctk_bbox(2).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(4).eq.1 .and. ny.ge.4 ) then + var(nx,ny,:) = 2.d0 * var(nx-1,ny-1,:) - var(nx-2,ny-2,:) + end if + + if ( cctk_bbox(2).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then + var(nx,:,1) = 2.d0 * var(nx-1,:,2) - var(nx-2,:,3) + end if + + if ( cctk_bbox(2).eq.1 .and. nx.ge.4 + $ .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then + var(nx,:,nz) = 2.d0 * var(nx-1,:,nz-1) - var(nx-2,:,nz-2) + end if + +C Remaining 2 in y=min. + if ( cctk_bbox(3).eq.1 .and. ny.ge.4 + $ .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then + var(:,1,1) = 2.d0 * var(:,2,2) - var(:,3,3) + end if + + if ( cctk_bbox(3).eq.1 .and. ny.ge.4 + $ .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then + var(:,1,nz) = 2.d0 * var(:,2,nz-1) - var(:,2,nz-2) + end if + +C Remaining 2 in y=ymax. + if ( cctk_bbox(4).eq.1 .and. ny.ge.4 + $ .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then + var(:,ny,1) = 2.d0 * var(:,ny-1,2) - var(:,ny-2,3) + end if + + if ( cctk_bbox(4).eq.1 .and. ny.ge.4 + $ .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then + var(:,ny,nz) = 2.d0 * var(:,ny-1,nz-1) - var(:,ny-2,nz-2) + end if + +C 8 corners. + + if (nx.ge.4 .and. ny.ge.4 .and. nz.ge.4) then + if (cctk_bbox(1).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(5).eq.1) then + var(1,1,1) = 2.d0*var(2,2,2) - var(3,3,3) + end if + if (cctk_bbox(1).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(6).eq.1) then + var(1,1,nz) = 2.d0*var(2,2,nz-1) - var(3,3,nz-2) + end if + if (cctk_bbox(1).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(5).eq.1) then + var(1,ny,1) = 2.d0*var(2,ny-1,2) - var(3,ny-2,3) + end if + if (cctk_bbox(1).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(6).eq.1) then + var(1,ny,nz) = 2.d0*var(2,ny-1,nz-1) - var(3,ny-2,nz-2) + end if + if (cctk_bbox(2).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(5).eq.1) then + var(nx,1,1) = 2.d0*var(nx-1,2,2) - var(nx-2,3,3) + end if + if (cctk_bbox(2).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(6).eq.1) then + var(nx,1,nz) = 2.d0*var(nx-1,2,nz-1) - var(nx-2,3,nz-2) + end if + if (cctk_bbox(2).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(5).eq.1) then + var(nx,ny,1) = 2.d0*var(nx-1,ny-1,2) - var(nx-2,ny-2,3) + end if + if (cctk_bbox(2).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(6).eq.1) then + var(nx,ny,nz) = 2.d0*var(nx-1,ny-1,nz-1) - var(nx-2,ny-2,nz-2) + end if + end if + + + return + end diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..c3a4561 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,15 @@ +# Main make.code.defn file for thorn Exact +# $Header$ + +# Source files in this directory +SRCS = exactinitialize.F exactgauge.F exactboundary.F \ + exactdata.F exactmetric.F \ + boostrotmetric.F finkelstein.F kerrschild.F flatschwarz.F \ + starschwarz.F bowl.F novikov.F minkowski.F fakebinary.F multiBH.F\ + exactblendbound.F exactcartblendbound.F \ + slice_initialize.F slice_evolve.F slice_data.F linextraponebound.F +# ComparisonSolutions.F Comparison.c + +# Subdirectories containing source files +SUBDIRS = + diff --git a/src/minkowski.F b/src/minkowski.F new file mode 100644 index 0000000..ed76b17 --- /dev/null +++ b/src/minkowski.F @@ -0,0 +1,46 @@ +#include "cctk.h" + + subroutine minkowski( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + + implicit none + +c Input. + + CCTK_REAL x, y, z, t + +c Output. + + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + + gdtt = -1.d0 + gdtx = 0.d0 + gdty = 0.d0 + gdtz = 0.d0 + gdxx = 1.d0 + gdyy = 1.d0 + gdzz = 1.d0 + gdxy = 0.d0 + gdyz = 0.d0 + gdzx = 0.d0 + + gutt = -1.d0 + gutx = 0.d0 + guty = 0.d0 + gutz = 0.d0 + guxx = 1.d0 + guyy = 1.d0 + guzz = 1.d0 + guxy = 0.d0 + guyz = 0.d0 + guzx = 0.d0 + + return + end diff --git a/src/multiBH.F b/src/multiBH.F new file mode 100644 index 0000000..51f9b25 --- /dev/null +++ b/src/multiBH.F @@ -0,0 +1,133 @@ +c======================================================================= +c23456789012345678901234567890123456789012345678901234567890123456789012 +c----------------------------------------------------------------------- +c by Hisaaki Shinkai shinkai@wurel.wustl.edu 19980603 +c----------------------------------------------------------------------- +c This is for maximally charged multi BH solutions such as +c Majumdar-Papapetrou (1947) or Kastor-Traschen (1993) solution. +c See also doc/KTsol.tex for brief review of this solution. +c----------------------------------------------------------------------- +c For usage: in your par file +c exactmodel = "multiBH" +c kt_hubble = 0.0 # 0.0 means MP solution +c kt_nBH = 2 # number of BHs (upto 4, currently) +c m_bh1,co_bh1x,co_bh1y,co_bh1z = 1.0, -2.0,0.0,0.0 # masses and +c m_bh2,co_bh2x,co_bh2y,co_bh2z = 1.0, 2.0,0.0,0.0 # locations +c----------------------------------------------------------------------- + +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine multi_BHs( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + + implicit none + + DECLARE_CCTK_PARAMETERS + +c Input. + CCTK_REAL x, y, z, t + +c Output. + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + + logical firstcall + + CCTK_REAL kt_r, kt_aa, kt_omega + CCTK_REAL one,zero + CCTK_REAL kt_xbh(10),kt_ybh(10),kt_zbh(10),kt_mbh(10) + + integer i + + data firstcall /.true./ + save firstcall,kt_xbh,kt_ybh,kt_zbh,kt_mbh + +c Get parameters of the exact solution. + + if (firstcall) then + + write(*,*) ' welcome to Kastor-Traschen (Majumdar-Papapetrou)' + + if(kt_nBH.ge.1) then + kt_xbh(1) = co_bh1x + kt_ybh(1) = co_bh1y + kt_zbh(1) = co_bh1z + kt_mbh(1) = m_bh1 + endif + + if(kt_nBH.ge.2) then + kt_xbh(2) = co_bh2x + kt_ybh(2) = co_bh2y + kt_zbh(2) = co_bh2z + kt_mbh(2) = m_bh2 + endif + + if(kt_nBH.ge.3) then + kt_xbh(3) = co_bh3x + kt_ybh(3) = co_bh3y + kt_zbh(3) = co_bh3z + kt_mbh(3) = m_bh3 + endif + + if(kt_nBH.ge.4) then + kt_xbh(4) = co_bh4x + kt_ybh(4) = co_bh4y + kt_zbh(4) = co_bh4z + kt_mbh(4) = m_bh4 + endif + + write(*,*) 'time=',t + write(*,*) ' mass BH(x,y,z) ' + + do i=1,kt_nBH + write(*,'(4e12.3)') kt_mbh(i),kt_xbh(i),kt_ybh(i),kt_zbh(i) + enddo + + firstcall = .false. + + end if + + one =1.0 + zero=0.0 + + kt_aa=exp(kt_hubble*t) + kt_omega=1.0 + + do i=1,kt_nBH + kt_r=sqrt((x-kt_xbh(i))**2+(y-kt_ybh(i))**2+(z-kt_zbh(i))**2) + kt_omega= kt_omega+kt_mbh(i)/kt_aa/kt_r + enddo + +c write(*,*) kt_omega,kt_aa + + gdtt = -1.0/kt_omega**2 + gdtx = zero + gdty = zero + gdtz = zero + gdxx = (kt_aa*kt_omega)**2 + gdyy = (kt_aa*kt_omega)**2 + gdzz = (kt_aa*kt_omega)**2 + gdxy = zero + gdyz = zero + gdzx = zero + + gutt = one/gdtt + gutx = zero + guty = zero + gutz = zero + guxx = one/gdxx + guyy = one/gdyy + guzz = one/gdzz + guxy = zero + guyz = zero + guzx = zero + + return + end diff --git a/src/novikov.F b/src/novikov.F new file mode 100644 index 0000000..0b5120c --- /dev/null +++ b/src/novikov.F @@ -0,0 +1,201 @@ +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine novikov( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + +c The metric given here corresponds to the novikov solution +c in isotropic coordinates, as presented first in Bruegman96 +c then in correct form in Cactus paper 1. This code is the code +c which was used for the comparisons in cactus paper 1, and is written +c by PW with input from BB. + + implicit none + + DECLARE_CCTK_PARAMETERS + +c Input. + + CCTK_REAL x, y, z, t + +c Output. + + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + +c Internal. + + logical firstcall + + CCTK_REAL mass,radius + CCTK_REAL r,c,psi4 + CCTK_REAL zero,one,two + + CCTK_REAL nov_dr_drmax, nov_rmax, nov_r + CCTK_REAL grr, gqq, detg + + CCTK_REAL psi4_o_r2 + + data firstcall /.true./ + save firstcall + +c Find {zero,one,two}. + + zero = 0.0D0 + one = 1.0D0 + two = 2.0D0 + +c Get parameters of the metric. (just the mass) + mass = 1.0D0 + +c Find r. + r = dsqrt(x**2 + y**2 + z**2) + +c Find conformal factor. + c = mass/(two*r) + psi4 = (one + c)**4 + +c Evaluate novikov stuff. Note abs(t) since the data is time +c symmetric (the metric is, at least...) + grr = nov_dr_drmax(abs(t),abs(r)) + gqq = nov_r(abs(t),abs(r)) + + grr = grr **2 + gqq = gqq**2 / nov_rmax(abs(r))**2 + + +c Find metric components. + psi4_o_r2 = psi4 / r**2 + + gdtt = - 1.0D0 + + gdtx = zero + gdty = zero + gdtz = zero + +c This is just straightforward spherical -> cartesian I hope... ;-) +c Note at t=0 (grr = gqq = 1) this gives the expected result +c (namely diagonal psi^4, since psi4_o_r2 = psi^4 / r^2) + gdxx = (grr * x**2 + gqq * (y**2 + z**2)) * psi4_o_r2 + gdyy = (grr * y**2 + gqq * (x**2 + z**2)) * psi4_o_r2 + gdzz = (grr * z**2 + gqq * (x**2 + y**2)) * psi4_o_r2 + + gdxy = (grr - gqq) * x * y * psi4_o_r2 + gdzx = (grr - gqq) * x * z * psi4_o_r2 + gdyz = (grr - gqq) * y * z * psi4_o_r2 + +c Find inverse metric. + gutt = one/gdtt + gutx = zero + guty = zero + gutz = zero + detg = gdtt*gdxx*gdyy*gdzz-gdtt*gdxx*gdyz**2/4.D0-gdtt*gdxy**2*gdz + $z/4.D0+gdtt*gdxy*gdzx*gdyz/4.D0-gdtt*gdzx**2*gdyy/4.D0 + guxx = -gdtt*(-4.D0*gdyy*gdzz+gdyz**2)/(4.D0 * detg) + guxy = -gdtt*(2.D0*gdxy*gdzz-gdzx*gdyz)/(4.D0 * detg) + guzx = gdtt*(gdxy*gdyz-2.D0*gdzx*gdyy)/(4.D0 * detg) + guyy = gdtt*(4.D0*gdxx*gdzz-gdzx**2)/(4.D0 * detg) + guyz = -gdtt*(2.D0*gdxx*gdyz-gdxy*gdzx)/(4.D0 * detg) + guzz = gdtt*(4.D0*gdxx*gdyy-gdxy**2)/(4.D0 * detg) + + guxx = one/psi4 + guyy = one/psi4 + guzz = one/psi4 + + guxy = zero + guyz = zero + guzx = zero + + return + end + +c These are functions which evaluate the novikov stuff. + +c dr/drmax + CCTK_REAL function nov_dr_drmax(tauin,rbarin) + + implicit none + + CCTK_REAL rbarin, tauin + CCTK_REAL rt, nov_r, rmaxt, nov_rmax + + rt = nov_r(tauin, rbarin) + rmaxt = nov_rmax(rbarin) + nov_dr_drmax = 1.5D0 - rt / (2.0D0 * rmaxt) + + $ 1.5D0 * sqrt(rmaxt / rt - 1.0D0) * + $ acos(sqrt(rt/rmaxt)) + + return + end + + + + +c Bisection to invert the function below. This is pretty crappy +c but it works. + CCTK_REAL function nov_r(tauin, rbarin) + implicit none +c input + CCTK_REAL tauin, rbarin + +c funtions + CCTK_REAL nov_rmax, nov_tau + +c temps + CCTK_REAL rg, drg, delt, ttmp, rmt + CCTK_REAL eps + integer nit + nit = 0 + delt = 1000.0D0 + rmt = nov_rmax(rbarin) + rg = rmt + drg = rg / 2.0D0 + eps = 1.d-6 * rmt + do while (delt .gt. eps .and. nit .lt. 100) + ttmp = nov_tau(rg, rmt) + delt = abs(tauin - ttmp) + if (delt .gt. eps) then + if (ttmp .gt. tauin .or. rg .lt. drg) then + rg = rg + drg +c Enforce upper bound + if (rg .gt. rmt) rg = rmt + drg = drg / 2.0D0 + else + rg = rg - drg + endif + endif +c write (*,*) rg, ttmp, tauin + nit = nit + 1 + enddo + if (nit .ge. 100) then + write (*,*) "Novikov: inversion did not converge" + endif + nov_r = rg + return + end + +c Evaluate tau as a function of r and rmax + CCTK_REAL function nov_tau(r, rmax) + implicit none + CCTK_REAL r, rmax + + nov_tau= rmax * dsqrt(0.5D0 * r * (1.0D0 - r / rmax)) + + $ 2.0D0 * (rmax / 2)**(3.0/2.0) * + $ acos (dsqrt(r/rmax)) + + return + end + +c Evaluate rmax as a function of rbar + CCTK_REAL function nov_rmax(rbar) + implicit none + CCTK_REAL rbar + nov_rmax = (1.0D0 + 2.0D0*rbar)**2 / (4.0D0 * rbar) + return + end diff --git a/src/slice_data.F b/src/slice_data.F new file mode 100644 index 0000000..33e0140 --- /dev/null +++ b/src/slice_data.F @@ -0,0 +1,327 @@ +C Extract Cauchy data from the slice x^A(x^i) stored in slicex, +C slicey, slicez, slicet, and calculate dx^A/dt. + +#include "cctk.h" +#include "cctk_arguments.h" + + subroutine slice_data(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + integer i, j, k, l, m, n, p, q, s + integer nx,ny,nz + + CCTK_REAL s1d(4,3), nd(4), nu(4), norm, gd(4,4), gu(4,4), g3(3,3), + $ gd_p(4,4), gd_m(4,4), gd1d(4,4,4), s2d(4,3,3), h3(3,3), + $ ex_eps, d3(3,3,3) + + CCTK_REAL dx,dy,dz + + parameter (ex_eps=1.d-6) + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + +C Sum over interior points on the slice. + + do k=2,nz-1 + do j=2,ny-1 + do i=2,nx-1 + +C Calculate first derivatives of slice coordinates. + + s1d(1,1) = 0.5d0 * (slicex(i+1,j,k) - slicex(i-1,j,k))/dx + s1d(1,2) = 0.5d0 * (slicex(i,j+1,k) - slicex(i,j-1,k))/dy + s1d(1,3) = 0.5d0 * (slicex(i,j,k+1) - slicex(i,j,k-1))/dz + + s1d(2,1) = 0.5d0 * (slicey(i+1,j,k) - slicey(i-1,j,k))/dx + s1d(2,2) = 0.5d0 * (slicey(i,j+1,k) - slicey(i,j-1,k))/dy + s1d(2,3) = 0.5d0 * (slicey(i,j,k+1) - slicey(i,j,k-1))/dz + + s1d(3,1) = 0.5d0 * (slicez(i+1,j,k) - slicez(i-1,j,k))/dx + s1d(3,2) = 0.5d0 * (slicez(i,j+1,k) - slicez(i,j-1,k))/dy + s1d(3,3) = 0.5d0 * (slicez(i,j,k+1) - slicez(i,j,k-1))/dz + + s1d(4,1) = 0.5d0 * (slicet(i+1,j,k) - slicet(i-1,j,k))/dx + s1d(4,2) = 0.5d0 * (slicet(i,j+1,k) - slicet(i,j-1,k))/dy + s1d(4,3) = 0.5d0 * (slicet(i,j,k+1) - slicet(i,j,k-1))/dz + +C Now we need the exact solution metric in the preferred coordinates x^A. + + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k), slicez(i,j,k), + $ slicet(i,j,k), + $ gd(4,4), gd(1,4), gd(2,4), gd(3,4), + $ gd(1,1), gd(2,2), gd(3,3), + $ gd(1,2), gd(2,3), gd(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + +C Calculate n^A and dx^A/dt + +#include "include/slice_normal.h" + +C Calculate g_ij = X^A,i X^B,j g_AB, at first in a temp. + + g3 = 0.d0 + + do p=1,3 + do q=p,3 + do m=1,4 + do n=1,4 + g3(p,q) = g3(p,q) + $ + s1d(m,p) * s1d(n,q) * gd(m,n) + end do + end do + end do + end do + + gxx(i,j,k) = g3(1,1) + gxy(i,j,k) = g3(1,2) + gxz(i,j,k) = g3(1,3) + gyy(i,j,k) = g3(2,2) + gyz(i,j,k) = g3(2,3) + gzz(i,j,k) = g3(3,3) + +C Calculate x^A,ij. Do this by hand, and loop over first index +C with an editor (hint for proofreading). + + s2d(1,1,1) = (slicex(i+1,j,k) + slicex(i-1,j,k) + $ - 2.d0 * slicex(i,j,k)) / dx**2 + s2d(1,2,2) = (slicex(i,j+1,k) + slicex(i,j-1,k) + $ - 2.d0 * slicex(i,j,k)) / dy**2 + s2d(1,3,3) = (slicex(i,j,k+1) + slicex(i,j,k-1) + $ - 2.d0 * slicex(i,j,k)) / dz**2 + s2d(1,1,2) = (slicex(i+1,j+1,k) - slicex(i-1,j+1,k) + $ - slicex(i+1,j-1,k) + slicex(i-1,j-1,k)) / (4.*dx*dy) + s2d(1,1,3) = (slicex(i+1,j,k+1) - slicex(i-1,j,k+1) + $ - slicex(i+1,j,k-1) + slicex(i-1,j,k-1)) / (4.*dy*dz) + s2d(1,2,3) = (slicex(i,j+1,k+1) - slicex(i,j-1,k+1) + $ - slicex(i,j+1,k-1) + slicex(i,j-1,k-1)) / (4.*dx*dz) + s2d(1,2,1) = s2d(1,1,2) + s2d(1,3,2) = s2d(1,2,3) + s2d(1,3,1) = s2d(1,1,3) + + s2d(2,1,1) = (slicey(i+1,j,k) + slicey(i-1,j,k) + $ - 2.d0 * slicey(i,j,k)) / dx**2 + s2d(2,2,2) = (slicey(i,j+1,k) + slicey(i,j-1,k) + $ - 2.d0 * slicey(i,j,k)) / dy**2 + s2d(2,3,3) = (slicey(i,j,k+1) + slicey(i,j,k-1) + $ - 2.d0 * slicey(i,j,k)) / dz**2 + s2d(2,1,2) = (slicey(i+1,j+1,k) - slicey(i-1,j+1,k) + $ - slicey(i+1,j-1,k) + slicey(i-1,j-1,k)) / (4.*dx*dy) + s2d(2,1,3) = (slicey(i+1,j,k+1) - slicey(i-1,j,k+1) + $ - slicey(i+1,j,k-1) + slicey(i-1,j,k-1)) / (4.*dy*dz) + s2d(2,2,3) = (slicey(i,j+1,k+1) - slicey(i,j-1,k+1) + $ - slicey(i,j+1,k-1) + slicey(i,j-1,k-1)) / (4.*dx*dz) + s2d(2,2,1) = s2d(2,1,2) + s2d(2,3,2) = s2d(2,2,3) + s2d(2,3,1) = s2d(2,1,3) + + s2d(3,1,1) = (slicez(i+1,j,k) + slicez(i-1,j,k) + $ - 2.d0 * slicez(i,j,k)) / dx**2 + s2d(3,2,2) = (slicez(i,j+1,k) + slicez(i,j-1,k) + $ - 2.d0 * slicez(i,j,k)) / dy**2 + s2d(3,3,3) = (slicez(i,j,k+1) + slicez(i,j,k-1) + $ - 2.d0 * slicez(i,j,k)) / dz**2 + s2d(3,1,2) = (slicez(i+1,j+1,k) - slicez(i-1,j+1,k) + $ - slicez(i+1,j-1,k) + slicez(i-1,j-1,k)) / (4.*dx*dy) + s2d(3,1,3) = (slicez(i+1,j,k+1) - slicez(i-1,j,k+1) + $ - slicez(i+1,j,k-1) + slicez(i-1,j,k-1)) / (4.*dy*dz) + s2d(3,2,3) = (slicez(i,j+1,k+1) - slicez(i,j-1,k+1) + $ - slicez(i,j+1,k-1) + slicez(i,j-1,k-1)) / (4.*dx*dz) + s2d(3,2,1) = s2d(3,1,2) + s2d(3,3,2) = s2d(3,2,3) + s2d(3,3,1) = s2d(3,1,3) + + s2d(4,1,1) = (slicet(i+1,j,k) + slicet(i-1,j,k) + $ - 2.d0 * slicet(i,j,k)) / dx**2 + s2d(4,2,2) = (slicet(i,j+1,k) + slicet(i,j-1,k) + $ - 2.d0 * slicet(i,j,k)) / dy**2 + s2d(4,3,3) = (slicet(i,j,k+1) + slicet(i,j,k-1) + $ - 2.d0 * slicet(i,j,k)) / dz**2 + s2d(4,1,2) = (slicet(i+1,j+1,k) - slicet(i-1,j+1,k) + $ - slicet(i+1,j-1,k) + slicet(i-1,j-1,k)) / (4.*dx*dy) + s2d(4,1,3) = (slicet(i+1,j,k+1) - slicet(i-1,j,k+1) + $ - slicet(i+1,j,k-1) + slicet(i-1,j,k-1)) / (4.*dy*dz) + s2d(4,2,3) = (slicet(i,j+1,k+1) - slicet(i,j-1,k+1) + $ - slicet(i,j+1,k-1) + slicet(i,j-1,k-1)) / (4.*dx*dz) + s2d(4,2,1) = s2d(4,1,2) + s2d(4,3,2) = s2d(4,2,3) + s2d(4,3,1) = s2d(4,1,3) + + +C Calculate g_AB,C. Need to sum explicitly over C. Do this with +C the editor. + + call exactmetric( + $ slicex(i,j,k)+ex_eps, slicey(i,j,k), slicez(i,j,k), + $ slicet(i,j,k), + $ gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4), + $ gd_p(1,1), gd_p(2,2), gd_p(3,3), + $ gd_p(1,2), gd_p(2,3), gd_p(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + call exactmetric( + $ slicex(i,j,k)-ex_eps, slicey(i,j,k), slicez(i,j,k), + $ slicet(i,j,k), + $ gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4), + $ gd_m(1,1), gd_m(2,2), gd_m(3,3), + $ gd_m(1,2), gd_m(2,3), gd_m(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + do m=1,4 + do n=m,4 + gd1d(m,n,1) = (gd_p(m,n) - gd_m(m,n)) / (2.*ex_eps) + end do + end do + + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k)+ex_eps, slicez(i,j,k), + $ slicet(i,j,k), + $ gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4), + $ gd_p(1,1), gd_p(2,2), gd_p(3,3), + $ gd_p(1,2), gd_p(2,3), gd_p(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k)-ex_eps, slicez(i,j,k), + $ slicet(i,j,k), + $ gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4), + $ gd_m(1,1), gd_m(2,2), gd_m(3,3), + $ gd_m(1,2), gd_m(2,3), gd_m(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + do m=1,4 + do n=m,4 + gd1d(m,n,2) = (gd_p(m,n) - gd_m(m,n)) / (2.*ex_eps) + end do + end do + + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k), slicez(i,j,k)+ex_eps, + $ slicet(i,j,k), + $ gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4), + $ gd_p(1,1), gd_p(2,2), gd_p(3,3), + $ gd_p(1,2), gd_p(2,3), gd_p(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k), slicez(i,j,k)-ex_eps, + $ slicet(i,j,k), + $ gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4), + $ gd_m(1,1), gd_m(2,2), gd_m(3,3), + $ gd_m(1,2), gd_m(2,3), gd_m(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + do m=1,4 + do n=m,4 + gd1d(m,n,3) = (gd_p(m,n) - gd_m(m,n)) / (2.*ex_eps) + end do + end do + + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k), slicez(i,j,k), + $ slicet(i,j,k)+ex_eps, + $ gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4), + $ gd_p(1,1), gd_p(2,2), gd_p(3,3), + $ gd_p(1,2), gd_p(2,3), gd_p(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + call exactmetric( + $ slicex(i,j,k), slicey(i,j,k), slicez(i,j,k), + $ slicet(i,j,k)-ex_eps, + $ gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4), + $ gd_m(1,1), gd_m(2,2), gd_m(3,3), + $ gd_m(1,2), gd_m(2,3), gd_m(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + do m=1,4 + do n=m,4 + gd1d(m,n,4) = (gd_p(m,n) - gd_m(m,n)) / (2.*ex_eps) + end do + end do + +C Fill in remaining components of metric derivative. + + do l=1,4 + do m=1,4 + do n=1,m-1 + gd1d(m,n,l) = gd1d(n,m,l) + end do + end do + end do + +C Calculate K_ij. + + h3 = 0.d0 + do p=1,3 + do q=p,3 + do l=1,4 + do m=1,4 + h3(p,q) = h3(p,q) + $ + s2d(l,p,q) * nu(m) * gd(l,m) + do n=1,4 + h3(p,q) = h3(p,q) + 0.5d0 * ( + $ s1d(l,p) * s1d(n,q) * nu(m) + $ + s1d(n,p) * s1d(m,q) * nu(l) + $ - s1d(l,p) * s1d(m,q) * nu(n) ) + $ * gd1d(l,m,n) + end do + end do + end do + end do + end do + + kxx(i,j,k) = h3(1,1) + kxy(i,j,k) = h3(1,2) + kxz(i,j,k) = h3(1,3) + kyy(i,j,k) = h3(2,2) + kyz(i,j,k) = h3(2,3) + kzz(i,j,k) = h3(3,3) + + end do + end do + end do + +C Synchronize and bound slicetmp2, which contains dx^A/dt. +C This is really just for output, and will eventually be redundant. + + call linextraponebound(CCTK_FARGUMENTS,slicetmp2x) + call synconefunc(slicetmp2x) + call linextraponebound(CCTK_FARGUMENTS,slicetmp2y) + call synconefunc(slicetmp2y) + call linextraponebound(CCTK_FARGUMENTS,slicetmp2z) + call synconefunc(slicetmp2z) + call linextraponebound(CCTK_FARGUMENTS,slicetmp2t) + call synconefunc(slicetmp2t) + + kxx = - kxx + kxy = - kxy + kxz = - kxz + kyy = - kyy + kyz = - kyz + kzz = - kzz + + return + end + + diff --git a/src/slice_evolve.F b/src/slice_evolve.F new file mode 100644 index 0000000..8a0f786 --- /dev/null +++ b/src/slice_evolve.F @@ -0,0 +1,132 @@ +C Evolve the slice in the exact spacetime. + +#include "cctk.h" +#include "cctk_arguments.h" + + subroutine slice_evolve(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + integer i,j,k,m,n + integer nx,ny,nz + + CCTK_REAL s1d(4,3), nd(4), nu(4), norm, gd(4,4), gu(4,4) + CCTK_REAL dx,dy,dz,dt + +C Grid parameters. + + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + + dx = cctk_delta_space(1) + dy = cctk_delta_space(2) + dz = cctk_delta_space(3) + + dt = cctk_delta_time + +C Lax, or forward in time, step. +C dxA/dt has previously been stored in slicetmp2 +C by slice_data. + + slicetmp1x = slicex + 0.5d0 * dt * slicetmp2x + slicetmp1y = slicey + 0.5d0 * dt * slicetmp2y + slicetmp1z = slicez + 0.5d0 * dt * slicetmp2z + slicetmp1t = slicet + 0.5d0 * dt * slicetmp2t + +C Synchronize and bound slice. + + call linextraponebound(CCTK_FARGUMENTS,slicetmp1x) + call synconefunc(slicetmp1x) + call linextraponebound(CCTK_FARGUMENTS,slicetmp1y) + call synconefunc(slicetmp1y) + call linextraponebound(CCTK_FARGUMENTS,slicetmp1z) + call synconefunc(slicetmp1z) + call linextraponebound(CCTK_FARGUMENTS,slicetmp1t) + call synconefunc(slicetmp1t) + +C Prepare leapfrog step. +C Sum over interior points on the slice, now at midpoint slicetmp1. + + do k=2,nz-1 + do j=2,ny-1 + do i=2,nx-1 + +C Calculate first derivatives of slice coordinates. + + s1d(1,1) = 0.5d0*(slicetmp1x(i+1,j,k) - slicetmp1x(i-1,j,k))/dx + s1d(1,2) = 0.5d0*(slicetmp1x(i,j+1,k) - slicetmp1x(i,j-1,k))/dy + s1d(1,3) = 0.5d0*(slicetmp1x(i,j,k+1) - slicetmp1x(i,j,k-1))/dz + + s1d(2,1) = 0.5d0*(slicetmp1y(i+1,j,k) - slicetmp1y(i-1,j,k))/dx + s1d(2,2) = 0.5d0*(slicetmp1y(i,j+1,k) - slicetmp1y(i,j-1,k))/dy + s1d(2,3) = 0.5d0*(slicetmp1y(i,j,k+1) - slicetmp1y(i,j,k-1))/dz + + s1d(3,1) = 0.5d0*(slicetmp1z(i+1,j,k) - slicetmp1z(i-1,j,k))/dx + s1d(3,2) = 0.5d0*(slicetmp1z(i,j+1,k) - slicetmp1z(i,j-1,k))/dy + s1d(3,3) = 0.5d0*(slicetmp1z(i,j,k+1) - slicetmp1z(i,j,k-1))/dz + + s1d(4,1) = 0.5d0*(slicetmp1t(i+1,j,k) - slicetmp1t(i-1,j,k))/dx + s1d(4,2) = 0.5d0*(slicetmp1t(i,j+1,k) - slicetmp1t(i,j-1,k))/dy + s1d(4,3) = 0.5d0*(slicetmp1t(i,j,k+1) - slicetmp1t(i,j,k-1))/dz + +C Now we need the exact solution metric in the preferred coordinates +C x^A. + call exactmetric( + $ slicetmp1x(i,j,k), slicetmp1y(i,j,k), slicetmp1z(i,j,k), + $ slicetmp1t(i,j,k), + $ gd(4,4), gd(1,4), gd(2,4), gd(3,4), + $ gd(1,1), gd(2,2), gd(3,3), + $ gd(1,2), gd(2,3), gd(1,3), + $ gu(4,4), gu(1,4), gu(2,4), gu(3,4), + $ gu(1,1), gu(2,2), gu(3,3), + $ gu(1,2), gu(2,3), gu(1,3)) + +C Calculate n^A and dx^A/dt +#include "include/slice_normal.h" + + end do + end do + end do + +C Synchronize and bound slicetmp2, which contains dx^A/dt. + + call linextraponebound(CCTK_FARGUMENTS,slicetmp2x) + call synconefunc(slicetmp2x) + call linextraponebound(CCTK_FARGUMENTS,slicetmp2y) + call synconefunc(slicetmp2y) + call linextraponebound(CCTK_FARGUMENTS,slicetmp2z) + call synconefunc(slicetmp2z) + call linextraponebound(CCTK_FARGUMENTS,slicetmp2t) + call synconefunc(slicetmp2t) + +C Leapfrog step. + + slicex = slicex + dt * slicetmp2x + slicey = slicey + dt * slicetmp2y + slicez = slicez + dt * slicetmp2z + slicet = slicet + dt * slicetmp2t + +C Synchronize and bound slice. + + call linextraponebound(CCTK_FARGUMENTS,slicex) + call synconefunc(slicex) + call linextraponebound(CCTK_FARGUMENTS,slicey) + call synconefunc(slicey) + call linextraponebound(CCTK_FARGUMENTS,slicez) + call synconefunc(slicez) + call linextraponebound(CCTK_FARGUMENTS,slicet) + call synconefunc(slicet) + +C Extract Cauchy data at the new position, and store dxA/dt +C for use in the next Lax step. + + call slice_data(CCTK_FARGUMENTS) + + return + end + + + diff --git a/src/slice_initialize.F b/src/slice_initialize.F new file mode 100644 index 0000000..8871387 --- /dev/null +++ b/src/slice_initialize.F @@ -0,0 +1,41 @@ + +#include "cctk.h" +#include "cctk_arguments.h" +#include "CactusEinstein/Einstein/src/Einstein.h" + + subroutine slice_initialize(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + + CCTK_REAL slice_gauss_ampl, slice_gauss_width + +C Initialize position of slice in exact solution spacetime. +C For now something simple. + + slicex = x + slicey = y + slicez = z + slicet = slice_gauss_ampl*exp(-(x**2 + y**2 + z**2) + $ /slice_gauss_width**2) + +C Calculate Cauchy data and dx^A/dt. + + call slice_data(CCTK_FARGUMENTS) + +C Extrapolate g, K, etc. to the boundary, where slice_data cannot +C calculate them, because it cannot take derivatives there. At all +C other timesteps cactus will do this automatically. + + call boundaries(CCTK_FARGUMENTS) + +C This wasted a day of work. PW found it... + + conformal_state = NOCONFORMAL_METRIC + psi = 1.0D0 + + return + end + + diff --git a/src/starschwarz.F b/src/starschwarz.F new file mode 100644 index 0000000..fa47b2f --- /dev/null +++ b/src/starschwarz.F @@ -0,0 +1,120 @@ +#include "cctk.h" +#include "cctk_parameters.h" + + subroutine starschwarz( + $ x, y, z, t, + $ gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx) + +c The metric given here corresponds to a constant +c density star, also known as a `Schwarzschild' star. +c It is obviously not a vacuum solution, so if someone +c wants to evolve it, we will need to add here the +c matter variables. For the moment, only the metric +c is given. +c +c The metric is given as a conformally flat metric. +c Turns out that in the original areal radius, the +c metric variables have a kink at the surface of the +c star, but they are smooth in the conformal form. +c +c Thanks to Philippos Papadopoulos for suggesting +c the use of this metric. + + implicit none + + DECLARE_CCTK_PARAMETERS + +c Input. + + CCTK_REAL x, y, z, t + +c Output. + + CCTK_REAL gdtt, gdtx, gdty, gdtz, + $ gdxx, gdyy, gdzz, gdxy, gdyz, gdzx, + $ gutt, gutx, guty, gutz, + $ guxx, guyy, guzz, guxy, guyz, guzx + +c Internal. + + logical firstcall + + CCTK_REAL mass,radius + CCTK_REAL r,c,psi4 + CCTK_REAL zero,one,two + + data firstcall /.true./ + save firstcall + +c Find {zero,one,two}. + + zero = 0.0D0 + one = 1.0D0 + two = 2.0D0 + +c Get parameters of the metric. + + if (firstcall) then + + mass = starschwarz_m + radius = starschwarz_r + + firstcall = .false. + + end if + +c Find r. + + r = dsqrt(x**2 + y**2 + z**2) + +c Find conformal factor. + + if (r.le.radius) then + + c = mass/(two*radius) + + psi4 = (one + c)**6/(one + c*(r/radius)**2)**2 + + else + + c = mass/(two*r) + + psi4 = (one + c)**4 + + end if + +c Find metric components. + + gdtt = - psi4 + + gdtx = zero + gdty = zero + gdtz = zero + + gdxx = psi4 + gdyy = psi4 + gdzz = psi4 + + gdxy = zero + gdyz = zero + gdzx = zero + +c Find inverse metric. + + gutx = zero + guty = zero + gutz = zero + + guxx = one/psi4 + guyy = one/psi4 + guzz = one/psi4 + + guxy = zero + guyz = zero + guzx = zero + + return + end |