From e9bd56216de3463d4c500faa3d8a7c053ba883e3 Mon Sep 17 00:00:00 2001 From: jthorn Date: Sun, 16 Jun 2002 18:13:32 +0000 Subject: [[This is a redo of my "cvs import" of 2002/06/11, this time using proper cvs operations (commit/delete/add) to preserve the full CVS history of this thorn.]] This is a major cleanup/revision of AEIThorns/Exact. Major user-visible changes: * major expansion of doc/documentation.tex * major expansion of documentation in param.ccl file * rename all parameters, systematize spacetime/coordinate/parameter names (there is a perl script in par/convert-pars.pl to convert old parameter files to the new names) * [from Mitica Vulcanov] many additions and fixes to cosmological solutions and Schwarzschild-Lemaitre * fix stress-energy tensor computations so they work -- before they were all disabled in CVS (INCLUDES lines were commented out in interface.ccl) due to requiring excessive friendship with evolution thorns and/or public parameters; new code copies parameters to restricted grid scalars, which Cactus automagically "pushes" to friends * added some more tests to testsuite, though these don't yet work fully Additional internal changes: * rename many Fortran subroutines (and a few C ones too) so their names start with the thorn name to reduce the chances of name collisions with other thorns * move all metrics to subdirectory so the main source directory isn't so cluttered * move two files containing subroutines which were never called (they didn't belong in this thorn, but somehow got into cvs by accident) into new archive/ directory * some (small) improvements in efficiency -- the exact_model parameter is now decoded from a keyword (string) to an integer once at INITIAL, and that integer tested by the stress-energy tensor code, rather than requiring a separate series of string tests at each grid point (!) like the old stress-energy tensor code did Added Files: (these are the two files containing subroutines which were never called) Comparison.c moved from ../src/Comparison.c ComparisonSolutions.F moved from ../src/ComparisonSolutions.F git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@99 e296648e-0e4f-0410-bd07-d597d9acff87 --- archive/Comparison.c | 332 ++++++++++++++++++++++++++++++++++++++++++ archive/ComparisonSolutions.F | 185 +++++++++++++++++++++++ 2 files changed, 517 insertions(+) create mode 100644 archive/Comparison.c create mode 100644 archive/ComparisonSolutions.F (limited to 'archive') diff --git a/archive/Comparison.c b/archive/Comparison.c new file mode 100644 index 0000000..5add331 --- /dev/null +++ b/archive/Comparison.c @@ -0,0 +1,332 @@ +/* Please god comment me ... */ +/* $Header$ */ + +#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; ingridFuncs;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;inpoints;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;inpoints;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;inpoints;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;inpoints;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/archive/ComparisonSolutions.F b/archive/ComparisonSolutions.F new file mode 100644 index 0000000..1c1fcd4 --- /dev/null +++ b/archive/ComparisonSolutions.F @@ -0,0 +1,185 @@ +C $Header$ + +#include "cctk.h" +#include "cctk_Arguments.h" + + subroutine ComparisonMetric(CCTK_ARGUMENTS, + $ gxx_ex, gxy_ex, gxz_ex, + $ gyy_ex, gyz_ex, gzz_ex) + + implicit none + + DECLARE_CCTK_ARGUMENTS + + 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_ARGUMENTS, + $ hxx_ex, hxy_ex, hxz_ex, + $ hyy_ex, hyz_ex, hzz_ex) + + implicit none + + DECLARE_CCTK_ARGUMENTS + + 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_ARGUMENTS, + $ alp_ex, betax_ex, betay_ex, betaz_ex) + + implicit none + + DECLARE_CCTK_ARGUMENTS + + 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 -- cgit v1.2.3