diff options
Diffstat (limited to 'archive/Comparison.c')
-rw-r--r-- | archive/Comparison.c | 332 |
1 files changed, 332 insertions, 0 deletions
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; 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); + +} + |