1037 lines
34 KiB
C
1037 lines
34 KiB
C
|
|
#include <iostream>
|
|
#include <iomanip>
|
|
#include <fstream>
|
|
#include <cstdlib>
|
|
#include <cstdio>
|
|
#include <string>
|
|
#include <cmath>
|
|
#include <new>
|
|
using namespace std;
|
|
|
|
#include "NullShellPatch2.h"
|
|
#include "Parallel.h"
|
|
#include "fmisc.h"
|
|
#include "misc.h"
|
|
#include "shellfunctions.h"
|
|
#include "NullEvol.h"
|
|
#include "NullNews.h"
|
|
#include "initial_null2.h"
|
|
#include "rungekutta4_rout.h"
|
|
#include "kodiss.h"
|
|
|
|
#define PI M_PI
|
|
|
|
#if 0
|
|
// for RT
|
|
void NullShellPatch2::Setup_Initial_Data(bool checkrun,double PhysTime)
|
|
{
|
|
if(checkrun)
|
|
{
|
|
}
|
|
else
|
|
{
|
|
MyList<ss_patch> *Pp=PatL;
|
|
while(Pp)
|
|
{
|
|
MyList<Block> *BL=Pp->data->blb;
|
|
while(BL)
|
|
{
|
|
Block *cg=BL->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
f_get_initial_null2(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
Pp->data->sst,Rmin);
|
|
// for Theta_AB
|
|
f_get_gauge_g00_K(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[g00->sgfn],Rmin);
|
|
}
|
|
if(BL == Pp->data->ble) break;
|
|
BL=BL->next;
|
|
}
|
|
Pp=Pp->next;
|
|
}
|
|
//Synchronize K
|
|
Synch(g00List,Symmetry,g00wt,1);
|
|
Pp=PatL;
|
|
int IONE=1;
|
|
while(Pp)
|
|
{
|
|
MyList<Block> *BP=Pp->data->blb;
|
|
int fngfs = Pp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
f_get_gauge_g00(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[g00->sgfn],Rmin,IONE);
|
|
}
|
|
if(BP==Pp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
Pp=Pp->next;
|
|
}
|
|
Synch(ThetaList,Symmetry,Thetawt,3);
|
|
}
|
|
}
|
|
#else
|
|
void NullShellPatch2::Setup_Initial_Data(bool checkrun, double PhysTime)
|
|
{
|
|
if (checkrun)
|
|
{
|
|
}
|
|
else
|
|
{
|
|
MyList<ss_patch> *Pp = PatL;
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BL = Pp->data->blb;
|
|
while (BL)
|
|
{
|
|
Block *cg = BL->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_get_initial_null3(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g220->sgfn], cg->fgfs[g230->sgfn], cg->fgfs[g330->sgfn],
|
|
Pp->data->sst, Rmin);
|
|
}
|
|
if (BL == Pp->data->ble)
|
|
break;
|
|
BL = BL->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
void NullShellPatch2::Step(double dT, double PhysTime, monitor *ErrorMonitor)
|
|
{
|
|
int iter_count = 0; // count RK4 substeps
|
|
int pre = 0, cor = 1;
|
|
int ERROR = 0;
|
|
double TT = PhysTime;
|
|
double neps = -0.05;
|
|
MyList<ss_patch> *sPp;
|
|
|
|
// Predictor
|
|
HyperSlice(dT, TT, ErrorMonitor, iter_count);
|
|
{
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
// rhs calculation
|
|
f_array_copy(cg->shape, cg->fgfs[g22_rhs->sgfn], cg->fgfs[Theta22->sgfn]);
|
|
f_array_copy(cg->shape, cg->fgfs[g23_rhs->sgfn], cg->fgfs[Theta23->sgfn]);
|
|
f_array_copy(cg->shape, cg->fgfs[g33_rhs->sgfn], cg->fgfs[Theta33->sgfn]);
|
|
f_kodis_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[g220->sgfn], cg->fgfs[g22_rhs->sgfn],
|
|
Thetawt[0], Symmetry, neps, sPp->data->sst);
|
|
f_kodis_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[g230->sgfn], cg->fgfs[g23_rhs->sgfn],
|
|
Thetawt[1], Symmetry, neps, sPp->data->sst);
|
|
f_kodis_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[g330->sgfn], cg->fgfs[g33_rhs->sgfn],
|
|
Thetawt[2], Symmetry, neps, sPp->data->sst);
|
|
f_rungekutta4_rout(cg->shape, dT, cg->fgfs[g220->sgfn], cg->fgfs[g22->sgfn], cg->fgfs[g22_rhs->sgfn],
|
|
iter_count);
|
|
f_rungekutta4_rout(cg->shape, dT, cg->fgfs[g230->sgfn], cg->fgfs[g23->sgfn], cg->fgfs[g23_rhs->sgfn],
|
|
iter_count);
|
|
f_rungekutta4_rout(cg->shape, dT, cg->fgfs[g330->sgfn], cg->fgfs[g33->sgfn], cg->fgfs[g33_rhs->sgfn],
|
|
iter_count);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
Synch(SynchList_pre, Symmetry, Thetawt, 3);
|
|
// Synch(SynchList_pre,Symmetry,g00wt,1);
|
|
|
|
// corrector
|
|
for (iter_count = 1; iter_count < 4; iter_count++)
|
|
{
|
|
// for RK4: t0, t0+dt/2, t0+dt/2, t0+dt;
|
|
if (iter_count == 1 || iter_count == 3)
|
|
TT += dT / 2;
|
|
HyperSlice(dT, TT, ErrorMonitor, iter_count);
|
|
{
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
// rhs calculation
|
|
f_array_copy(cg->shape, cg->fgfs[g221->sgfn], cg->fgfs[Theta22->sgfn]);
|
|
f_array_copy(cg->shape, cg->fgfs[g231->sgfn], cg->fgfs[Theta23->sgfn]);
|
|
f_array_copy(cg->shape, cg->fgfs[g331->sgfn], cg->fgfs[Theta33->sgfn]);
|
|
f_kodis_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[g22->sgfn], cg->fgfs[g221->sgfn],
|
|
Thetawt[0], Symmetry, neps, sPp->data->sst);
|
|
f_kodis_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[g23->sgfn], cg->fgfs[g231->sgfn],
|
|
Thetawt[1], Symmetry, neps, sPp->data->sst);
|
|
f_kodis_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2], cg->fgfs[g33->sgfn], cg->fgfs[g331->sgfn],
|
|
Thetawt[2], Symmetry, neps, sPp->data->sst);
|
|
f_rungekutta4_rout(cg->shape, dT, cg->fgfs[g220->sgfn], cg->fgfs[g221->sgfn], cg->fgfs[g22_rhs->sgfn],
|
|
iter_count);
|
|
f_rungekutta4_rout(cg->shape, dT, cg->fgfs[g230->sgfn], cg->fgfs[g231->sgfn], cg->fgfs[g23_rhs->sgfn],
|
|
iter_count);
|
|
f_rungekutta4_rout(cg->shape, dT, cg->fgfs[g330->sgfn], cg->fgfs[g331->sgfn], cg->fgfs[g33_rhs->sgfn],
|
|
iter_count);
|
|
}
|
|
if (iter_count < 3)
|
|
cg->swapList(SynchList_cor, SynchList_pre, myrank);
|
|
else
|
|
{
|
|
cg->swapList(StateList, SynchList_cor, myrank);
|
|
cg->swapList(OldStateList, SynchList_cor, myrank);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
if (iter_count < 3)
|
|
Synch(SynchList_pre, Symmetry, Thetawt, 3);
|
|
else
|
|
Synch(StateList, Symmetry, Thetawt, 3);
|
|
// if( iter_count < 3 ) Synch(SynchList_pre,Symmetry,g00wt,1);
|
|
// else Synch(StateList,Symmetry,g00wt,1);
|
|
}
|
|
}
|
|
// really ODEs, so we do not need Synch in this routine at all
|
|
#if 0
|
|
void NullShellPatch2::HyperSlice(double dT,double PhysTime,monitor *ErrorMonitor,int RK_count)
|
|
{
|
|
int ERROR=0;
|
|
Null_Boundary(PhysTime);
|
|
#if 1
|
|
MyList<ss_patch> *sPp;
|
|
|
|
// evolve g01
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
if(RK_count==0)
|
|
{
|
|
if(f_NullEvol_g01(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g01->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of g01 in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if(f_NullEvol_g01(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g22->sgfn],cg->fgfs[g23->sgfn],cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g01->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of g01 in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
Synch(g01List,Symmetry,g01wt,1);
|
|
if(RK_count==3) Dump_Data(g01List,0,PhysTime,dT);
|
|
//check error information
|
|
{int erh=ERROR;MPI_Allreduce(&erh,&ERROR,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); }
|
|
if(ERROR)
|
|
{
|
|
if(RK_count==0) Dump_Data(StateList,0,PhysTime,dT);
|
|
else Dump_Data(SynchList_pre,0,PhysTime,dT);
|
|
if(myrank == 0)
|
|
{
|
|
if(ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
else
|
|
cout<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
|
|
// evolve p02, p03, g02 and g03
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
if(RK_count==0)
|
|
{
|
|
if(f_NullEvol_pg0A(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of pg0A in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if(f_NullEvol_pg0A(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g22->sgfn],cg->fgfs[g23->sgfn],cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of pg0A in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
Synch(pg0AList,Symmetry,pg0Awt,2);
|
|
if(RK_count==3) Dump_Data(pg0AList,0,PhysTime,dT);
|
|
//check error information
|
|
{int erh=ERROR;MPI_Allreduce(&erh,&ERROR,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); }
|
|
if(ERROR)
|
|
{
|
|
Dump_Data(g01List,0,PhysTime,dT);
|
|
if(myrank == 0)
|
|
{
|
|
if(ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
else
|
|
cout<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#if 0
|
|
// for gauge variable g00
|
|
{
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
if(RK_count==0)
|
|
{
|
|
f_get_gauge_g00_real(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[g00->sgfn],Rmin);
|
|
}
|
|
else
|
|
{
|
|
f_get_gauge_g00_real(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g22->sgfn],cg->fgfs[g23->sgfn],cg->fgfs[g33->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[g00->sgfn],Rmin);
|
|
}
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
Synch(g00List,Symmetry,g00wt,1);
|
|
if(RK_count==3) Dump_Data(g00List,0,PhysTime,dT);
|
|
}
|
|
// evolve ThetaAB
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
if(RK_count==0)
|
|
{
|
|
if(f_NullEvol_Theta2(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g00->sgfn],cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of ThetaAB in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if(f_NullEvol_Theta2(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g22->sgfn],cg->fgfs[g23->sgfn],cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g00->sgfn],cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of ThetaAB in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
Synch(ThetaList,Symmetry,Thetawt,3);
|
|
if(RK_count==3) Dump_Data(ThetaList,0,PhysTime,dT);
|
|
//check error information
|
|
{int erh=ERROR;MPI_Allreduce(&erh,&ERROR,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); }
|
|
if(ERROR)
|
|
{
|
|
Dump_Data(pg0AList,0,PhysTime,dT);
|
|
if(myrank == 0)
|
|
{
|
|
if(ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
else
|
|
cout<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#elif 1
|
|
// evolve ThetaAB and g00
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
if(RK_count==0)
|
|
{
|
|
if(f_NullEvol_Thetag00(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g00->sgfn],cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of ThetaAB in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if(f_NullEvol_Thetag00(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g22->sgfn],cg->fgfs[g23->sgfn],cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g00->sgfn],cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],Rmin))
|
|
{
|
|
cout<<"find NaN of ThetaAB in NullShell domain: sst = "<<sPp->data->sst<<", ("<<cg->bbox[0]<<":"<<cg->bbox[3]<<","
|
|
<<cg->bbox[1]<<":"<<cg->bbox[4]<<","<<cg->bbox[2]<<":"<<cg->bbox[5]<<")"<<endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
Synch(ThetaList,Symmetry,Thetawt,3);
|
|
if(RK_count==3) Dump_Data(ThetaList,0,PhysTime,dT);
|
|
Synch(g00List,Symmetry,g00wt,1);
|
|
if(RK_count==3) Dump_Data(g00List,0,PhysTime,dT);
|
|
//check error information
|
|
{int erh=ERROR;MPI_Allreduce(&erh,&ERROR,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); }
|
|
if(ERROR)
|
|
{
|
|
Dump_Data(pg0AList,0,PhysTime,dT);
|
|
if(myrank == 0)
|
|
{
|
|
if(ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
else
|
|
cout<<"find NaN in beta on NullShell Patches at t = "<<PhysTime<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#endif
|
|
#endif
|
|
}
|
|
#else
|
|
void NullShellPatch2::HyperSlice(double dT, double PhysTime, monitor *ErrorMonitor, int RK_count)
|
|
{
|
|
int ERROR = 0;
|
|
Null_Boundary(PhysTime);
|
|
|
|
MyList<ss_patch> *sPp;
|
|
|
|
// evolve g01
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
if (RK_count == 0)
|
|
{
|
|
if (f_NullEvol_g01(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g220->sgfn], cg->fgfs[g230->sgfn], cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g01->sgfn], Rmin))
|
|
{
|
|
cout << "find NaN of g01 in NullShell domain: sst = " << sPp->data->sst << ", (" << cg->bbox[0] << ":" << cg->bbox[3] << ","
|
|
<< cg->bbox[1] << ":" << cg->bbox[4] << "," << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (f_NullEvol_g01(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g22->sgfn], cg->fgfs[g23->sgfn], cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g01->sgfn], Rmin))
|
|
{
|
|
cout << "find NaN of g01 in NullShell domain: sst = " << sPp->data->sst << ", (" << cg->bbox[0] << ":" << cg->bbox[3] << ","
|
|
<< cg->bbox[1] << ":" << cg->bbox[4] << "," << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
Synch(g01List, Symmetry, g01wt, 1);
|
|
// if(RK_count==3) Dump_Data(g01List,0,PhysTime,dT);
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
if (RK_count == 0)
|
|
Dump_Data(StateList, 0, PhysTime, dT);
|
|
else
|
|
Dump_Data(SynchList_pre, 0, PhysTime, dT);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in beta on NullShell Patches at t = " << PhysTime << endl;
|
|
else
|
|
cout << "find NaN in beta on NullShell Patches at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
// evolve p02, p03, g02 and g03
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
if (RK_count == 0)
|
|
{
|
|
if (f_NullEvol_pg0A(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g220->sgfn], cg->fgfs[g230->sgfn], cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn], cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g02->sgfn], cg->fgfs[g03->sgfn], Rmin))
|
|
{
|
|
cout << "find NaN of pg0A in NullShell domain: sst = " << sPp->data->sst << ", (" << cg->bbox[0] << ":" << cg->bbox[3] << ","
|
|
<< cg->bbox[1] << ":" << cg->bbox[4] << "," << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (f_NullEvol_pg0A(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g22->sgfn], cg->fgfs[g23->sgfn], cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn], cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g02->sgfn], cg->fgfs[g03->sgfn], Rmin))
|
|
{
|
|
cout << "find NaN of pg0A in NullShell domain: sst = " << sPp->data->sst << ", (" << cg->bbox[0] << ":" << cg->bbox[3] << ","
|
|
<< cg->bbox[1] << ":" << cg->bbox[4] << "," << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
Synch(pg0AList, Symmetry, pg0Awt, 2);
|
|
// if(RK_count==3) Dump_Data(pg0AList,0,PhysTime,dT);
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Dump_Data(g01List, 0, PhysTime, dT);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in beta on NullShell Patches at t = " << PhysTime << endl;
|
|
else
|
|
cout << "find NaN in beta on NullShell Patches at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
// for gauge variable g00
|
|
{
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_get_g00_with_t(PhysTime, cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g00->sgfn], Rmin, sPp->data->sst);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
// if(RK_count==3) Dump_Data(g00List,0,PhysTime,dT);
|
|
}
|
|
// evolve ThetaAB
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
if (RK_count == 0)
|
|
{
|
|
if (f_NullEvol_Theta2(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g220->sgfn], cg->fgfs[g230->sgfn], cg->fgfs[g330->sgfn],
|
|
cg->fgfs[g00->sgfn], cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn], cg->fgfs[g03->sgfn],
|
|
cg->fgfs[p02->sgfn], cg->fgfs[p03->sgfn],
|
|
cg->fgfs[Theta22->sgfn], cg->fgfs[Theta23->sgfn], cg->fgfs[Theta33->sgfn], Rmin))
|
|
{
|
|
cout << "find NaN of ThetaAB in NullShell domain: sst = " << sPp->data->sst << ", (" << cg->bbox[0] << ":" << cg->bbox[3] << ","
|
|
<< cg->bbox[1] << ":" << cg->bbox[4] << "," << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (f_NullEvol_Theta2(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g22->sgfn], cg->fgfs[g23->sgfn], cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g00->sgfn], cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn], cg->fgfs[g03->sgfn],
|
|
cg->fgfs[p02->sgfn], cg->fgfs[p03->sgfn],
|
|
cg->fgfs[Theta22->sgfn], cg->fgfs[Theta23->sgfn], cg->fgfs[Theta33->sgfn], Rmin))
|
|
{
|
|
cout << "find NaN of ThetaAB in NullShell domain: sst = " << sPp->data->sst << ", (" << cg->bbox[0] << ":" << cg->bbox[3] << ","
|
|
<< cg->bbox[1] << ":" << cg->bbox[4] << "," << cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
}
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
Synch(ThetaList, Symmetry, Thetawt, 3);
|
|
// if(RK_count==3) Dump_Data(ThetaList,0,PhysTime,dT);
|
|
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Dump_Data(pg0AList, 0, PhysTime, dT);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor && ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in beta on NullShell Patches at t = " << PhysTime << endl;
|
|
else
|
|
cout << "find NaN in beta on NullShell Patches at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
#if 0
|
|
void NullShellPatch2::Null_Boundary(double PhysTime)
|
|
{
|
|
MyList<ss_patch> *sPp;
|
|
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
f_get_null_boundary2(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g22->sgfn],cg->fgfs[g23->sgfn],cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn],cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g02->sgfn],cg->fgfs[g03->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],Rmin);
|
|
// for Theta_AB
|
|
f_get_gauge_g00_K(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[g00->sgfn],Rmin);
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
// boundary for Theta_AB
|
|
//Synchronize K
|
|
Synch(g00List,Symmetry,g00wt,1);
|
|
sPp=PatL;
|
|
int IZEO=1;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
f_get_gauge_g00(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[Theta22->sgfn],cg->fgfs[Theta23->sgfn],cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[g00->sgfn],Rmin,IZEO);
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
Synch(ThetaList,Symmetry,Thetawt,3);
|
|
//Synch(ThetaList,Symmetry,g00wt,1);
|
|
// boundary condition is independent of angular direction, do not need synch
|
|
// Synch(pg0AList,Symmetry,pg0Awt,2,-1);
|
|
// Synch(g00List,Symmetry,g00wt,1,-1);
|
|
// Synch(ThetaList,Symmetry,Thetawt,3,-1);
|
|
}
|
|
#else
|
|
void NullShellPatch2::Null_Boundary(double PhysTime)
|
|
{
|
|
MyList<ss_patch> *sPp;
|
|
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_get_null_boundary3(PhysTime, cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g22->sgfn], cg->fgfs[g23->sgfn], cg->fgfs[g33->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn], cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g02->sgfn], cg->fgfs[g03->sgfn],
|
|
cg->fgfs[Theta22->sgfn], cg->fgfs[Theta23->sgfn], cg->fgfs[Theta33->sgfn], Rmin, sPp->data->sst);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
/*
|
|
// check Synch
|
|
Synch(g01List,Symmetry,g01wt,1);
|
|
Dump_Data(g01List,0,PhysTime,1);
|
|
Synch(pg0AList,Symmetry,pg0Awt,2);
|
|
Dump_Data(pg0AList,0,PhysTime,1);
|
|
Synch(StateList,Symmetry,Thetawt,3);
|
|
Dump_Data(StateList,0,PhysTime,1);
|
|
Synch(ThetaList,Symmetry,Thetawt,3);
|
|
Dump_Data(ThetaList,0,PhysTime,1);
|
|
if(myrank==0) MPI_Abort(MPI_COMM_WORLD,1);
|
|
*/
|
|
}
|
|
// 0: real L2 norm; 1: root mean squar
|
|
#define L2m 0
|
|
double NullShellPatch2::Error_Check(double PhysTime)
|
|
{
|
|
MyList<ss_patch> *sPp;
|
|
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_get_null_boundary3(PhysTime, cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[g221->sgfn], cg->fgfs[g231->sgfn], cg->fgfs[g331->sgfn],
|
|
cg->fgfs[g01->sgfn],
|
|
cg->fgfs[p02->sgfn], cg->fgfs[p03->sgfn],
|
|
cg->fgfs[g22_rhs->sgfn], cg->fgfs[g03->sgfn],
|
|
cg->fgfs[Theta22->sgfn], cg->fgfs[Theta23->sgfn], cg->fgfs[Theta33->sgfn], Rmin, sPp->data->sst);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
|
|
double tvf, dtvf = 0;
|
|
int tN, dtN = 0;
|
|
int BDW = ghost_width, OBDW = overghost;
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_array_subtract(cg->shape, cg->fgfs[g22_rhs->sgfn], cg->fgfs[g02->sgfn]);
|
|
#if (L2m == 0)
|
|
f_l2normhelper_sh(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
sPp->data->bbox[0], sPp->data->bbox[1], sPp->data->bbox[2], sPp->data->bbox[3], sPp->data->bbox[4], sPp->data->bbox[5],
|
|
cg->fgfs[g22_rhs->sgfn], tvf, BDW, OBDW, Symmetry);
|
|
#elif (L2m == 1)
|
|
f_l2normhelper_sh_rms(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
sPp->data->bbox[0], sPp->data->bbox[1], sPp->data->bbox[2], sPp->data->bbox[3], sPp->data->bbox[4], sPp->data->bbox[5],
|
|
cg->fgfs[g22_rhs->sgfn], tvf, BDW, OBDW, Symmetry, dtN);
|
|
dtN += dtN;
|
|
#endif
|
|
|
|
dtvf += tvf;
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
|
|
// Dump_Data(RHSList,0,PhysTime,1);
|
|
// Dump_Data(ThetaList,0,PhysTime,1);
|
|
// if(myrank==0) MPI_Abort(MPI_COMM_WORLD,1);
|
|
|
|
MPI_Allreduce(&dtvf, &tvf, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
|
#if (L2m == 0)
|
|
tvf = sqrt(tvf);
|
|
#elif (L2m == 1)
|
|
MPI_Allreduce(&dtN, &tN, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
|
tvf = sqrt(tvf / tN);
|
|
#endif
|
|
|
|
return tvf;
|
|
}
|
|
#undef L2m
|
|
#endif
|
|
|
|
void NullShellPatch2::Compute_News(double PhysTime)
|
|
{
|
|
MyList<ss_patch> *sPp;
|
|
|
|
// get omega and dtomega
|
|
// for RT
|
|
#if 0
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
f_get_omega_and_dtomega_pre(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[omega->sgfn],cg->fgfs[dtomega->sgfn],Rmin);
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
// Synch
|
|
{
|
|
MyList<var> * DG_List;
|
|
DG_List=new MyList<var>(omega);
|
|
Synch(DG_List,Symmetry,g00wt,1);
|
|
DG_List->clearList();
|
|
DG_List=new MyList<var>(dtomega);
|
|
Synch(DG_List,Symmetry,g00wt,1);
|
|
DG_List->clearList();
|
|
}
|
|
// get dtomega
|
|
sPp=PatL;
|
|
while(sPp)
|
|
{
|
|
MyList<Block> *BP=sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while(BP)
|
|
{
|
|
Block *cg=BP->data;
|
|
if(myrank == cg->rank)
|
|
{
|
|
f_get_dtomega(cg->shape,cg->X[0],cg->X[1],cg->X[2],
|
|
cg->fgfs[g220->sgfn],cg->fgfs[g230->sgfn],cg->fgfs[g330->sgfn],
|
|
cg->fgfs[omega->sgfn],cg->fgfs[dtomega->sgfn],Rmin);
|
|
}
|
|
if(BP==sPp->data->ble) break;
|
|
BP=BP->next;
|
|
}
|
|
sPp=sPp->next;
|
|
}
|
|
// Synch
|
|
{
|
|
MyList<var> * DG_List;
|
|
DG_List=new MyList<var>(dtomega);
|
|
Synch(DG_List,Symmetry,g00wt,1);
|
|
DG_List->clearList();
|
|
}
|
|
#else
|
|
// for linear wave
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_get_omega_and_dtomega_LN(PhysTime, cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[omega->sgfn], cg->fgfs[dtomega->sgfn], Rmin, sPp->data->sst);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
// Synch
|
|
{
|
|
MyList<var> *DG_List;
|
|
DG_List = new MyList<var>(omega);
|
|
Synch(DG_List, Symmetry, g00wt, 1);
|
|
DG_List->clearList();
|
|
DG_List = new MyList<var>(dtomega);
|
|
Synch(DG_List, Symmetry, g00wt, 1);
|
|
DG_List->clearList();
|
|
}
|
|
#endif
|
|
// calculate News
|
|
sPp = PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
int fngfs = sPp->data->fngfs;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
f_get_null_news2(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
cg->fgfs[omega->sgfn], cg->fgfs[dtomega->sgfn],
|
|
cg->fgfs[g00->sgfn], cg->fgfs[g01->sgfn],
|
|
cg->fgfs[g02->sgfn], cg->fgfs[g03->sgfn],
|
|
cg->fgfs[g220->sgfn], cg->fgfs[g230->sgfn], cg->fgfs[g330->sgfn],
|
|
cg->fgfs[Theta22->sgfn], cg->fgfs[Theta23->sgfn], cg->fgfs[Theta33->sgfn],
|
|
cg->fgfs[RNews->sgfn], cg->fgfs[INews->sgfn], Rmin, sPp->data->sst);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|