1943 lines
68 KiB
C
1943 lines
68 KiB
C
// includes, system
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include <iostream>
|
|
#include <sys/time.h>
|
|
|
|
#ifdef RESULT_CHECK
|
|
#include <fstream>
|
|
#endif
|
|
|
|
// include BSSN class files
|
|
#include "macrodef.h"
|
|
#include "fmisc.h"
|
|
#include "bssn_gpu_class.h"
|
|
#include "bssn_rhs.h"
|
|
#include "enforce_algebra.h"
|
|
#include "rungekutta4_rout.h"
|
|
#include "sommerfeld_rout.h"
|
|
|
|
// include gpu files
|
|
#include "bssn_gpu.h"
|
|
|
|
#if (PSTR == 0)
|
|
#if 1
|
|
void bssn_class::Step_GPU(int lev, int YN)
|
|
{
|
|
setpbh(BH_num, Porg0, Mass, BH_num_input);
|
|
|
|
double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
|
|
|
|
// new code 2013-2-15, zjcao
|
|
#if (MAPBH == 1)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
for (int ith = 0; ith < 3; ith++)
|
|
Porg1[ithBH][ith] = Porg0[ithBH][ith] + Porg_rhs[ithBH][ith] * dT_lev;
|
|
if (Symmetry > 0)
|
|
Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
|
|
Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
|
|
}
|
|
if (!finite(Porg1[ithBH][0]) || !finite(Porg1[ithBH][1]) || !finite(Porg1[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "predictor step finds NaN for BH's position from ("
|
|
<< Porg0[ithBH][0] << "," << Porg0[ithBH][1] << "," << Porg0[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
|
|
// data analysis part
|
|
// Warning NOTE: the variables1 are used as temp storege room
|
|
if (lev == a_lev)
|
|
{
|
|
AnalysisStuff(lev, dT_lev);
|
|
}
|
|
#endif
|
|
|
|
#ifdef With_AHF
|
|
AH_Step_Find(lev, dT_lev);
|
|
#endif
|
|
bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
|
|
double ndeps = numepss;
|
|
if (lev < GH->movls)
|
|
ndeps = numepsb;
|
|
double TRK4 = PhysTime;
|
|
int iter_count = 0; // count RK4 substeps
|
|
int pre = 0, cor = 1;
|
|
int ERROR = 0;
|
|
|
|
MyList<ss_patch> *sPp;
|
|
// Predictor
|
|
MyList<Patch> *Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
#if (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
|
|
cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_FIRST_TIME))
|
|
{
|
|
cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," << cg->bbox[1] << ":" << cg->bbox[4] << ","
|
|
<< cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
#if (SommerType == 0)
|
|
#ifndef WithShell
|
|
if (lev == 0) // sommerfeld indeed
|
|
f_sommerfeld_routbam(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
cg->fgfs[varlrhs->data->sgfn],
|
|
cg->fgfs[varl0->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
#endif
|
|
#endif
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
#ifndef WithShell
|
|
if (lev > 0) // fix BD point
|
|
#endif
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
|
|
#if (SommerType == 1)
|
|
#warning "shell part still bam type"
|
|
if (lev == 0) // Shibata type sommerfeld
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], varl0->data->SoA,
|
|
Symmetry, pre);
|
|
#endif
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny);
|
|
}
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in state variables at t = " << PhysTime << ", lev = " << lev << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
#ifdef WithShell
|
|
// evolve Shell Patches
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->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 (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
|
|
cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs_ss(CALLED_BY_STEP, myrank, RHS_SS_PARA_CALLED_FIRST_TIME))
|
|
{
|
|
cout << "find NaN in Shell 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;
|
|
}
|
|
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
// sommerfeld indeed for outter boudary while fix BD for inner boundary
|
|
f_sommerfeld_routbam_ss(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[varlrhs->data->sgfn],
|
|
cg->fgfs[varl0->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
#if 0
|
|
// check rhs
|
|
{
|
|
SH->Dump_Data(RHSList,0,PhysTime,dT_lev);
|
|
if(myrank == 0)
|
|
{
|
|
cout<<"check rhs"<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
|
|
if (ERROR)
|
|
{
|
|
SH->Dump_Data(StateList, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in state variables on Shell Patches at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
|
|
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->Synch(SynchList_pre, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "Shell stuff synchronization used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
#if (MAPBH == 0)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][0], Porg[ithBH][0], Porg_rhs[ithBH][0], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][1], Porg[ithBH][1], Porg_rhs[ithBH][1], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][2], Porg[ithBH][2], Porg_rhs[ithBH][2], iter_count);
|
|
if (Symmetry > 0)
|
|
Porg[ithBH][2] = fabs(Porg[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg[ithBH][0] = fabs(Porg[ithBH][0]);
|
|
Porg[ithBH][1] = fabs(Porg[ithBH][1]);
|
|
}
|
|
if (!finite(Porg[ithBH][0]) || !finite(Porg[ithBH][1]) || !finite(Porg[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "predictor step finds NaN for BH's position from ("
|
|
<< Porg0[ithBH][0] << "," << Porg0[ithBH][1] << "," << Porg0[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
// data analysis part
|
|
// Warning NOTE: the variables1 are used as temp storege room
|
|
if (lev == a_lev)
|
|
{
|
|
AnalysisStuff(lev, dT_lev);
|
|
}
|
|
#endif
|
|
|
|
// 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)
|
|
TRK4 += dT_lev / 2;
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
#if (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#elif (AGM == 1)
|
|
if (iter_count == 3)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_THEN))
|
|
{
|
|
cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," << cg->bbox[1] << ":" << cg->bbox[4] << ","
|
|
<< cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
#if (SommerType == 0)
|
|
#ifndef WithShell
|
|
if (lev == 0) // sommerfeld indeed
|
|
f_sommerfeld_routbam(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
cg->fgfs[varl1->data->sgfn],
|
|
cg->fgfs[varl->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
#endif
|
|
#endif
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
#ifndef WithShell
|
|
if (lev > 0) // fix BD point
|
|
#endif
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
|
|
#if (SommerType == 1)
|
|
if (lev == 1) // shibata type sommerfeld
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varl1->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
#endif
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varl1 = varl1->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny);
|
|
}
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
|
|
if (ERROR)
|
|
{
|
|
Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in RK4 substep#" << iter_count << " variables at t = " << PhysTime << ", lev = " << lev << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
#ifdef WithShell
|
|
// evolve Shell Patches
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->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 (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#elif (AGM == 1)
|
|
if (iter_count == 3)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs_ss(CALLED_BY_STEP, myrank, RHS_SS_PARA_CALLED_THEN))
|
|
{
|
|
cout << "find NaN in Shell 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;
|
|
}
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
// sommerfeld indeed for outter boudary while fix BD for inner boundary
|
|
f_sommerfeld_routbam_ss(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[varl1->data->sgfn],
|
|
cg->fgfs[varl->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varl1 = varl1->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
SH->Dump_Data(SynchList_pre, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN on Shell Patches in RK4 substep#" << iter_count << " variables at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
|
|
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->Synch(SynchList_cor, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "Shell stuff synchronization used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
#if (MAPBH == 0)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg, Porg1, Sfx, Sfy, Sfz, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][0], Porg1[ithBH][0], Porg_rhs[ithBH][0], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][1], Porg1[ithBH][1], Porg_rhs[ithBH][1], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][2], Porg1[ithBH][2], Porg_rhs[ithBH][2], iter_count);
|
|
if (Symmetry > 0)
|
|
Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
|
|
Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
|
|
}
|
|
if (!finite(Porg1[ithBH][0]) || !finite(Porg1[ithBH][1]) || !finite(Porg1[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << iter_count << " corrector step finds NaN for BH's position from ("
|
|
<< Porg[ithBH][0] << "," << Porg[ithBH][1] << "," << Porg[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
|
|
// swap time level
|
|
if (iter_count < 3)
|
|
{
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(SynchList_pre, SynchList_cor, myrank);
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(SynchList_pre, SynchList_cor, myrank);
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
#if (MAPBH == 0)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
Porg[ithBH][0] = Porg1[ithBH][0];
|
|
Porg[ithBH][1] = Porg1[ithBH][1];
|
|
Porg[ithBH][2] = Porg1[ithBH][2];
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
}
|
|
#if (RPS == 0)
|
|
// mesh refinement boundary part
|
|
RestrictProlong(lev, YN, BB);
|
|
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->CS_Inter(SynchList_cor, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "CS_Inter used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
#endif
|
|
// note the data structure before update
|
|
// SynchList_cor 1 -----------
|
|
//
|
|
// StateList 0 -----------
|
|
//
|
|
// OldStateList old -----------
|
|
// update
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(StateList, SynchList_cor, myrank);
|
|
cg->swapList(OldStateList, SynchList_cor, myrank);
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
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 0
|
|
// check StateList
|
|
{
|
|
SH->Dump_Data(StateList,0,PhysTime,dT_lev);
|
|
if(myrank == 0)
|
|
{
|
|
cout<<"check StateList"<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
#endif
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
Porg0[ithBH][0] = Porg1[ithBH][0];
|
|
Porg0[ithBH][1] = Porg1[ithBH][1];
|
|
Porg0[ithBH][2] = Porg1[ithBH][2];
|
|
}
|
|
}
|
|
}
|
|
#else // #if 1
|
|
// ICN for bam comparison
|
|
void bssn_class::Step_GPU(int lev, int YN)
|
|
{
|
|
double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
|
|
#ifdef With_AHF
|
|
AH_Step_Find(lev, dT_lev);
|
|
#endif
|
|
bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
|
|
double ndeps = numepss;
|
|
if (lev < GH->movls)
|
|
ndeps = numepsb;
|
|
double TRK4 = PhysTime;
|
|
int iter_count = 0; // count RK4 substeps
|
|
int pre = 0, cor = 1;
|
|
int ERROR = 0;
|
|
|
|
MyList<ss_patch> *sPp;
|
|
// Predictor
|
|
MyList<Patch> *Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
#if (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
|
|
cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_FIRST_TIME))
|
|
{
|
|
cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," << cg->bbox[1] << ":" << cg->bbox[4] << ","
|
|
<< cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
#ifndef WithShell
|
|
if (lev == 0) // sommerfeld indeed
|
|
f_sommerfeld_routbam(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
cg->fgfs[varlrhs->data->sgfn],
|
|
cg->fgfs[varl0->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
#endif
|
|
f_icn_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
#ifndef WithShell
|
|
if (lev > 0) // fix BD point
|
|
#endif
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny);
|
|
}
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in state variables at t = " << PhysTime << ", lev = " << lev << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
#ifdef WithShell
|
|
// evolve Shell Patches
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->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 (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
|
|
cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs_ss(CALLED_BY_STEP, myrank, RHS_SS_PARA_CALLED_FIRST_TIME))
|
|
{
|
|
cout << "find NaN in Shell 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;
|
|
}
|
|
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
// sommerfeld indeed for outter boudary while fix BD for inner boundary
|
|
f_sommerfeld_routbam_ss(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[varlrhs->data->sgfn],
|
|
cg->fgfs[varl0->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
f_icn_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
#if 0
|
|
// check rhs
|
|
{
|
|
SH->Dump_Data(RHSList,0,PhysTime,dT_lev);
|
|
if(myrank == 0)
|
|
{
|
|
cout<<"check rhs"<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
SH->Dump_Data(StateList, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in state variables on Shell Patches at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
|
|
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->Synch(SynchList_pre, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "Shell stuff synchronization used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
f_icn_scalar(dT_lev, Porg0[ithBH][0], Porg[ithBH][0], Porg_rhs[ithBH][0], iter_count);
|
|
f_icn_scalar(dT_lev, Porg0[ithBH][1], Porg[ithBH][1], Porg_rhs[ithBH][1], iter_count);
|
|
f_icn_scalar(dT_lev, Porg0[ithBH][2], Porg[ithBH][2], Porg_rhs[ithBH][2], iter_count);
|
|
if (Symmetry > 0)
|
|
Porg[ithBH][2] = fabs(Porg[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg[ithBH][0] = fabs(Porg[ithBH][0]);
|
|
Porg[ithBH][1] = fabs(Porg[ithBH][1]);
|
|
}
|
|
if (!finite(Porg[ithBH][0]) || !finite(Porg[ithBH][1]) || !finite(Porg[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "predictor step finds NaN for BH's position from ("
|
|
<< Porg0[ithBH][0] << "," << Porg0[ithBH][1] << "," << Porg0[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
// data analysis part
|
|
// Warning NOTE: the variables1 are used as temp storege room
|
|
if (lev == a_lev)
|
|
{
|
|
AnalysisStuff(lev, dT_lev);
|
|
}
|
|
// corrector
|
|
for (iter_count = 1; iter_count < 3; iter_count++)
|
|
{
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
#if (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#elif (AGM == 1)
|
|
if (iter_count == 3)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_THEN))
|
|
{
|
|
cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," << cg->bbox[1] << ":" << cg->bbox[4] << ","
|
|
<< cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
#ifndef WithShell
|
|
if (lev == 0) // sommerfeld indeed
|
|
f_sommerfeld_routbam(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
cg->fgfs[varl1->data->sgfn],
|
|
cg->fgfs[varl->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
#endif
|
|
f_icn_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
#ifndef WithShell
|
|
if (lev > 0) // fix BD point
|
|
#endif
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varl1 = varl1->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny);
|
|
}
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in RK4 substep#" << iter_count << " variables at t = " << PhysTime << ", lev = " << lev << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
#ifdef WithShell
|
|
// evolve Shell Patches
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->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 (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#elif (AGM == 1)
|
|
if (iter_count == 3)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs_ss(CALLED_BY_STEP, myrank, RHS_SS_PARA_CALLED_THEN))
|
|
{
|
|
cout << "find NaN in Shell 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;
|
|
}
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
// sommerfeld indeed for outter boudary while fix BD for inner boundary
|
|
f_sommerfeld_routbam_ss(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[varl1->data->sgfn],
|
|
cg->fgfs[varl->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varl1 = varl1->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
SH->Dump_Data(SynchList_pre, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN on Shell Patches in RK4 substep#" << iter_count << " variables at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
#endif
|
|
|
|
Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
|
|
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->Synch(SynchList_cor, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "Shell stuff synchronization used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg, Porg1, Sfx, Sfy, Sfz, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
f_icn_scalar(dT_lev, Porg0[ithBH][0], Porg1[ithBH][0], Porg_rhs[ithBH][0], iter_count);
|
|
f_icn_scalar(dT_lev, Porg0[ithBH][1], Porg1[ithBH][1], Porg_rhs[ithBH][1], iter_count);
|
|
f_icn_scalar(dT_lev, Porg0[ithBH][2], Porg1[ithBH][2], Porg_rhs[ithBH][2], iter_count);
|
|
if (Symmetry > 0)
|
|
Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
|
|
Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
|
|
}
|
|
if (!finite(Porg1[ithBH][0]) || !finite(Porg1[ithBH][1]) || !finite(Porg1[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << iter_count << " corrector step finds NaN for BH's position from ("
|
|
<< Porg[ithBH][0] << "," << Porg[ithBH][1] << "," << Porg[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
// swap time level
|
|
if (iter_count < 3)
|
|
{
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(SynchList_pre, SynchList_cor, myrank);
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(SynchList_pre, SynchList_cor, myrank);
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
#endif
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
Porg[ithBH][0] = Porg1[ithBH][0];
|
|
Porg[ithBH][1] = Porg1[ithBH][1];
|
|
Porg[ithBH][2] = Porg1[ithBH][2];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
#if (RPS == 0)
|
|
// mesh refinement boundary part
|
|
RestrictProlong(lev, YN, BB);
|
|
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->CS_Inter(SynchList_cor, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "CS_Inter used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
|
|
#endif
|
|
// note the data structure before update
|
|
// SynchList_cor 1 -----------
|
|
//
|
|
// StateList 0 -----------
|
|
//
|
|
// OldStateList old -----------
|
|
// update
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(StateList, SynchList_cor, myrank);
|
|
cg->swapList(OldStateList, SynchList_cor, myrank);
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
#ifdef WithShell
|
|
if (lev == 0)
|
|
{
|
|
sPp = SH->PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
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 0
|
|
// check StateList
|
|
{
|
|
SH->Dump_Data(StateList,0,PhysTime,dT_lev);
|
|
if(myrank == 0)
|
|
{
|
|
cout<<"check StateList"<<endl;
|
|
MPI_Abort(MPI_COMM_WORLD,1);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
#endif
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
Porg0[ithBH][0] = Porg1[ithBH][0];
|
|
Porg0[ithBH][1] = Porg1[ithBH][1];
|
|
Porg0[ithBH][2] = Porg1[ithBH][2];
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
|
|
#elif (PSTR == 1)
|
|
void bssn_class::Step_GPU(int lev, int YN)
|
|
{
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"start Step");
|
|
|
|
setpbh(BH_num, Porg0, Mass, BH_num_input);
|
|
|
|
double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
|
|
|
|
// new code 2013-2-15, zjcao
|
|
#if (MAPBH == 1)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
for (int ith = 0; ith < 3; ith++)
|
|
Porg1[ithBH][ith] = Porg0[ithBH][ith] + Porg_rhs[ithBH][ith] * dT_lev;
|
|
if (Symmetry > 0)
|
|
Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
|
|
Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
|
|
}
|
|
if (!finite(Porg1[ithBH][0]) || !finite(Porg1[ithBH][1]) || !finite(Porg1[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "predictor step finds NaN for BH's position from ("
|
|
<< Porg0[ithBH][0] << "," << Porg0[ithBH][1] << "," << Porg0[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
#endif //(MAPBH == 1)
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Predictor");
|
|
|
|
#ifdef With_AHF
|
|
AH_Step_Find(lev, dT_lev);
|
|
#endif
|
|
bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
|
|
double ndeps = numepss;
|
|
if (lev < GH->movls)
|
|
ndeps = numepsb;
|
|
double TRK4 = PhysTime;
|
|
int iter_count = 0; // count RK4 substeps
|
|
int pre = 0, cor = 1;
|
|
int ERROR = 0;
|
|
|
|
MyList<ss_patch> *sPp;
|
|
// Predictor
|
|
MyList<Patch> *Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
#if (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
|
|
cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_FIRST_TIME))
|
|
{
|
|
cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," << cg->bbox[1] << ":" << cg->bbox[4] << ","
|
|
<< cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
#if (SommerType == 0)
|
|
#ifndef WithShell
|
|
if (lev == 0) // sommerfeld indeed
|
|
f_sommerfeld_routbam(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
cg->fgfs[varlrhs->data->sgfn],
|
|
cg->fgfs[varl0->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
#endif
|
|
#endif
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
#ifndef WithShell
|
|
if (lev > 0) // fix BD point
|
|
#endif
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
|
|
#if (SommerType == 1)
|
|
#warning "shell part still bam type"
|
|
if (lev == 0) // Shibata type sommerfeld
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], varl0->data->SoA,
|
|
Symmetry, pre);
|
|
#endif
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny);
|
|
}
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after Predictor rhs calculation");
|
|
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, GH->Commlev[lev]);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Parallel::Dump_Data(GH->PatL[lev], StateList, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in state variables at t = " << PhysTime << ", lev = " << lev << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Predictor sync");
|
|
|
|
Parallel::Sync(GH->PatL[lev], SynchList_pre, Symmetry);
|
|
|
|
#if (MAPBH == 0)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg0, Porg_rhs, Sfx0, Sfy0, Sfz0, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][0], Porg[ithBH][0], Porg_rhs[ithBH][0], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][1], Porg[ithBH][1], Porg_rhs[ithBH][1], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][2], Porg[ithBH][2], Porg_rhs[ithBH][2], iter_count);
|
|
if (Symmetry > 0)
|
|
Porg[ithBH][2] = fabs(Porg[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg[ithBH][0] = fabs(Porg[ithBH][0]);
|
|
Porg[ithBH][1] = fabs(Porg[ithBH][1]);
|
|
}
|
|
if (!finite(Porg[ithBH][0]) || !finite(Porg[ithBH][1]) || !finite(Porg[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "predictor step finds NaN for BH's position from ("
|
|
<< Porg0[ithBH][0] << "," << Porg0[ithBH][1] << "," << Porg0[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Corrector");
|
|
|
|
// corrector
|
|
for (iter_count = 1; iter_count < 4; iter_count++)
|
|
{
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"head of Corrector");
|
|
|
|
// for RK4: t0, t0+dt/2, t0+dt/2, t0+dt;
|
|
if (iter_count == 1 || iter_count == 3)
|
|
TRK4 += dT_lev / 2;
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
if (myrank == cg->rank)
|
|
{
|
|
#if (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#elif (AGM == 1)
|
|
if (iter_count == 3)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs(CALLED_BY_STEP, myrank, RHS_PARA_CALLED_THEN))
|
|
{
|
|
cout << "find NaN in domain: (" << cg->bbox[0] << ":" << cg->bbox[3] << "," << cg->bbox[1] << ":" << cg->bbox[4] << ","
|
|
<< cg->bbox[2] << ":" << cg->bbox[5] << ")" << endl;
|
|
ERROR = 1;
|
|
}
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
#if (SommerType == 0)
|
|
#ifndef WithShell
|
|
if (lev == 0) // sommerfeld indeed
|
|
f_sommerfeld_routbam(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
cg->fgfs[varl1->data->sgfn],
|
|
cg->fgfs[varl->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
#endif
|
|
#endif
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
#ifndef WithShell
|
|
if (lev > 0) // fix BD point
|
|
#endif
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
|
|
#if (SommerType == 1)
|
|
if (lev == 1) // shibata type sommerfeld
|
|
f_sommerfeld_rout(cg->shape, cg->X[0], cg->X[1], cg->X[2],
|
|
Pp->data->bbox[0], Pp->data->bbox[1], Pp->data->bbox[2], Pp->data->bbox[3], Pp->data->bbox[4], Pp->data->bbox[5],
|
|
dT_lev, cg->fgfs[phi0->sgfn],
|
|
cg->fgfs[Lap0->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varl1->data->sgfn], varl0->data->SoA,
|
|
Symmetry, cor);
|
|
#endif
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varl1 = varl1->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny);
|
|
}
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Corrector error check");
|
|
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, GH->Commlev[lev]);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
Parallel::Dump_Data(GH->PatL[lev], SynchList_pre, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in RK4 substep#" << iter_count << " variables at t = " << PhysTime << ", lev = " << lev << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Corrector sync");
|
|
|
|
Parallel::Sync(GH->PatL[lev], SynchList_cor, Symmetry);
|
|
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after Corrector sync");
|
|
|
|
#if (MAPBH == 0)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
compute_Porg_rhs(Porg, Porg1, Sfx, Sfy, Sfz, lev);
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][0], Porg1[ithBH][0], Porg_rhs[ithBH][0], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][1], Porg1[ithBH][1], Porg_rhs[ithBH][1], iter_count);
|
|
f_rungekutta4_scalar(dT_lev, Porg0[ithBH][2], Porg1[ithBH][2], Porg_rhs[ithBH][2], iter_count);
|
|
if (Symmetry > 0)
|
|
Porg1[ithBH][2] = fabs(Porg1[ithBH][2]);
|
|
if (Symmetry == 2)
|
|
{
|
|
Porg1[ithBH][0] = fabs(Porg1[ithBH][0]);
|
|
Porg1[ithBH][1] = fabs(Porg1[ithBH][1]);
|
|
}
|
|
if (!finite(Porg1[ithBH][0]) || !finite(Porg1[ithBH][1]) || !finite(Porg1[ithBH][2]))
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << iter_count << " corrector step finds NaN for BH's position from ("
|
|
<< Porg[ithBH][0] << "," << Porg[ithBH][1] << "," << Porg[ithBH][2] << ")" << endl;
|
|
|
|
MyList<var> *DG_List = new MyList<var>(Sfx0);
|
|
DG_List->insert(Sfx0);
|
|
DG_List->insert(Sfy0);
|
|
DG_List->insert(Sfz0);
|
|
Parallel::Dump_Data(GH->PatL[lev], DG_List, 0, PhysTime, dT_lev);
|
|
DG_List->clearList();
|
|
}
|
|
}
|
|
}
|
|
misc::tillherecheck(GH->Commlev[lev], GH->start_rank[lev], "after Corrector of black hole position");
|
|
#endif
|
|
|
|
// swap time level
|
|
if (iter_count < 3)
|
|
{
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(SynchList_pre, SynchList_cor, myrank);
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"after pre cor swap");
|
|
|
|
#if (MAPBH == 0)
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
Porg[ithBH][0] = Porg1[ithBH][0];
|
|
Porg[ithBH][1] = Porg1[ithBH][1];
|
|
Porg[ithBH][2] = Porg1[ithBH][2];
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"tail of corrector");
|
|
}
|
|
#if (RPS == 0)
|
|
// mesh refinement boundary part
|
|
misc::tillherecheck(GH->Commlev[lev], GH->start_rank[lev], "before RestrictProlong");
|
|
RestrictProlong(lev, YN, BB);
|
|
#endif
|
|
// note the data structure before update
|
|
// SynchList_cor 1 -----------
|
|
//
|
|
// StateList 0 -----------
|
|
//
|
|
// OldStateList old -----------
|
|
// update
|
|
Pp = GH->PatL[lev];
|
|
while (Pp)
|
|
{
|
|
MyList<Block> *BP = Pp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(StateList, SynchList_cor, myrank);
|
|
cg->swapList(OldStateList, SynchList_cor, myrank);
|
|
if (BP == Pp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
Pp = Pp->next;
|
|
}
|
|
// for black hole position
|
|
if (BH_num > 0 && lev == GH->levels - 1)
|
|
{
|
|
for (int ithBH = 0; ithBH < BH_num; ithBH++)
|
|
{
|
|
Porg0[ithBH][0] = Porg1[ithBH][0];
|
|
Porg0[ithBH][1] = Porg1[ithBH][1];
|
|
Porg0[ithBH][2] = Porg1[ithBH][2];
|
|
// if(myrank==GH->start_rank[lev]) cout<<Porg0[ithBH][0]<<","<<Porg0[ithBH][1]<<","<<Porg0[ithBH][2]<<endl;
|
|
}
|
|
}
|
|
|
|
// if(myrank==GH->start_rank[lev]) cout<<GH->mylev<<endl;
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"complet GH Step");
|
|
}
|
|
#endif // PSTR == ?
|
|
|
|
//--------------------------With Shell--------------------------
|
|
|
|
#ifdef WithShell
|
|
void bssn_class::SHStep()
|
|
{
|
|
int lev = 0;
|
|
// #if (PSTR == 1)
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"start Step");
|
|
// #endif
|
|
|
|
setpbh(BH_num, Porg0, Mass, BH_num_input);
|
|
|
|
double dT_lev = dT * pow(0.5, Mymax(lev, trfls));
|
|
|
|
// #if (PSTR == 1)
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Predictor");
|
|
// #endif
|
|
|
|
#ifdef With_AHF
|
|
AH_Step_Find(lev, dT_lev);
|
|
#endif
|
|
bool BB = fgt(PhysTime, StartTime, dT_lev / 2);
|
|
double ndeps = numepss;
|
|
if (lev < GH->movls)
|
|
ndeps = numepsb;
|
|
double TRK4 = PhysTime;
|
|
int iter_count = 0; // count RK4 substeps
|
|
int pre = 0, cor = 1;
|
|
int ERROR = 0;
|
|
|
|
MyList<ss_patch> *sPp;
|
|
// Predictor
|
|
sPp = SH->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 (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx0->sgfn], cg->fgfs[gxy0->sgfn], cg->fgfs[gxz0->sgfn], cg->fgfs[gyy0->sgfn], cg->fgfs[gyz0->sgfn], cg->fgfs[gzz0->sgfn],
|
|
cg->fgfs[Axx0->sgfn], cg->fgfs[Axy0->sgfn], cg->fgfs[Axz0->sgfn], cg->fgfs[Ayy0->sgfn], cg->fgfs[Ayz0->sgfn], cg->fgfs[Azz0->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs_ss(RHS_SS_PARA_CALLED_FIRST_TIME))
|
|
{
|
|
cout << "find NaN in Shell 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;
|
|
}
|
|
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
// sommerfeld indeed for outter boudary while fix BD for inner boundary
|
|
f_sommerfeld_routbam_ss(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[varlrhs->data->sgfn],
|
|
cg->fgfs[varl0->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi->sgfn], chitiny);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
|
|
#if (PSTR == 1)
|
|
// misc::tillherecheck(GH->Commlev[lev],GH->start_rank[lev],"before Predictor's error check");
|
|
#endif
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
|
|
if (ERROR)
|
|
{
|
|
SH->Dump_Data(StateList, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN in state variables on Shell Patches at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->Synch(SynchList_pre, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "Shell stuff synchronization used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
|
|
// 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)
|
|
TRK4 += dT_lev / 2;
|
|
|
|
{
|
|
sPp = SH->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 (AGM == 0)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#elif (AGM == 1)
|
|
if (iter_count == 3)
|
|
f_enforce_ga(cg->shape,
|
|
cg->fgfs[gxx->sgfn], cg->fgfs[gxy->sgfn], cg->fgfs[gxz->sgfn], cg->fgfs[gyy->sgfn], cg->fgfs[gyz->sgfn], cg->fgfs[gzz->sgfn],
|
|
cg->fgfs[Axx->sgfn], cg->fgfs[Axy->sgfn], cg->fgfs[Axz->sgfn], cg->fgfs[Ayy->sgfn], cg->fgfs[Ayz->sgfn], cg->fgfs[Azz->sgfn]);
|
|
#endif
|
|
|
|
if (gpu_rhs_ss(RHS_SS_PARA_CALLED_THEN))
|
|
{
|
|
cout << "find NaN in Shell 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;
|
|
}
|
|
// rk4 substep and boundary
|
|
{
|
|
MyList<var> *varl0 = StateList, *varl = SynchList_pre, *varl1 = SynchList_cor, *varlrhs = RHSList; // we do not check the correspondence here
|
|
while (varl0)
|
|
{
|
|
// sommerfeld indeed for outter boudary while fix BD for inner boundary
|
|
f_sommerfeld_routbam_ss(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[varl1->data->sgfn],
|
|
cg->fgfs[varl->data->sgfn], varl0->data->propspeed, varl0->data->SoA,
|
|
Symmetry);
|
|
|
|
f_rungekutta4_rout(cg->shape, dT_lev, cg->fgfs[varl0->data->sgfn], cg->fgfs[varl1->data->sgfn], cg->fgfs[varlrhs->data->sgfn],
|
|
iter_count);
|
|
|
|
varl0 = varl0->next;
|
|
varl = varl->next;
|
|
varl1 = varl1->next;
|
|
varlrhs = varlrhs->next;
|
|
}
|
|
}
|
|
f_lowerboundset(cg->shape, cg->fgfs[phi1->sgfn], chitiny);
|
|
}
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
// check error information
|
|
{
|
|
int erh = ERROR;
|
|
MPI_Allreduce(&erh, &ERROR, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
|
}
|
|
if (ERROR)
|
|
{
|
|
SH->Dump_Data(SynchList_pre, 0, PhysTime, dT_lev);
|
|
if (myrank == 0)
|
|
{
|
|
if (ErrorMonitor->outfile)
|
|
ErrorMonitor->outfile << "find NaN on Shell Patches in RK4 substep#" << iter_count << " variables at t = " << PhysTime << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
}
|
|
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->Synch(SynchList_cor, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "Shell stuff synchronization used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
|
|
sPp = SH->PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(SynchList_pre, SynchList_cor, myrank);
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
#if (RPS == 0)
|
|
{
|
|
clock_t prev_clock, curr_clock;
|
|
if (myrank == 0)
|
|
curr_clock = clock();
|
|
SH->CS_Inter(SynchList_cor, Symmetry);
|
|
if (myrank == 0)
|
|
{
|
|
prev_clock = curr_clock;
|
|
curr_clock = clock();
|
|
cout << "CS_Inter used " << (double)(curr_clock - prev_clock) / ((double)CLOCKS_PER_SEC) << " seconds!" << endl;
|
|
}
|
|
}
|
|
#endif
|
|
// note the data structure before update
|
|
// SynchList_cor 1 -----------
|
|
//
|
|
// StateList 0 -----------
|
|
//
|
|
// OldStateList old -----------
|
|
// update
|
|
sPp = SH->PatL;
|
|
while (sPp)
|
|
{
|
|
MyList<Block> *BP = sPp->data->blb;
|
|
while (BP)
|
|
{
|
|
Block *cg = BP->data;
|
|
cg->swapList(StateList, SynchList_cor, myrank);
|
|
cg->swapList(OldStateList, SynchList_cor, myrank);
|
|
if (BP == sPp->data->ble)
|
|
break;
|
|
BP = BP->next;
|
|
}
|
|
sPp = sPp->next;
|
|
}
|
|
}
|
|
d
|
|
#endif // withshell
|