509 lines
18 KiB
C
509 lines
18 KiB
C
|
|
#ifdef newc
|
|
#include <iostream>
|
|
#include <iomanip>
|
|
#include <fstream>
|
|
#include <cstdlib>
|
|
#include <cstdio>
|
|
#include <string>
|
|
#include <cmath>
|
|
#include <map>
|
|
using namespace std;
|
|
#else
|
|
#include <iostream.h>
|
|
#include <iomanip.h>
|
|
#include <fstream.h>
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include <string.h>
|
|
#include <math.h>
|
|
#include <map.h>
|
|
#endif
|
|
|
|
#include <mpi.h>
|
|
|
|
#include "misc.h"
|
|
#include "macrodef.h"
|
|
|
|
#ifndef ABEtype
|
|
#error "not define ABEtype"
|
|
#endif
|
|
|
|
#if (ABEtype == 0)
|
|
|
|
#ifdef USE_GPU
|
|
#include "bssn_gpu_class.h"
|
|
#else
|
|
#include "bssn_class.h"
|
|
#endif
|
|
|
|
#elif (ABEtype == 1)
|
|
#include "bssnEScalar_class.h"
|
|
|
|
#elif (ABEtype == 2)
|
|
#include "Z4c_class.h"
|
|
|
|
#elif (ABEtype == 3)
|
|
#include "bssnEM_class.h"
|
|
|
|
#else
|
|
#error "not recognized ABEtype"
|
|
#endif
|
|
|
|
namespace parameters
|
|
{
|
|
map<string, int> int_par;
|
|
map<string, double> dou_par;
|
|
map<string, string> str_par;
|
|
}
|
|
|
|
//=================================================================================================
|
|
//=================================================================================================
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
int myrank = 0, nprocs = 1;
|
|
MPI_Init(&argc, &argv);
|
|
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
|
|
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
|
|
|
|
double Begin_clock, End_clock;
|
|
if (myrank == 0)
|
|
{
|
|
Begin_clock = MPI_Wtime();
|
|
}
|
|
|
|
if (argc > 1)
|
|
{
|
|
string sttr(argv[1]);
|
|
parameters::str_par.insert(map<string, string>::value_type("inputpar", sttr));
|
|
}
|
|
else
|
|
{
|
|
string sttr("input.par");
|
|
parameters::str_par.insert(map<string, string>::value_type("inputpar", sttr));
|
|
}
|
|
|
|
int checkrun;
|
|
char checkfilename[50];
|
|
int ID_type;
|
|
int Steps;
|
|
double StartTime, TotalTime;
|
|
double AnasTime, DumpTime, d2DumpTime, CheckTime;
|
|
double Courant;
|
|
double numepss, numepsb, numepsh;
|
|
int Symmetry;
|
|
int a_lev, maxl, decn;
|
|
double maxrex, drex;
|
|
// read parameter from file
|
|
{
|
|
map<string, string>::iterator iter;
|
|
string out_dir;
|
|
const int LEN = 256;
|
|
char pline[LEN];
|
|
string str, sgrp, skey, sval;
|
|
int sind;
|
|
char pname[50];
|
|
iter = parameters::str_par.find("inputpar");
|
|
if (iter != parameters::str_par.end())
|
|
{
|
|
out_dir = iter->second;
|
|
sprintf(pname, "%s", out_dir.c_str());
|
|
}
|
|
else
|
|
{
|
|
cout << "Error inputpar" << endl;
|
|
exit(0);
|
|
}
|
|
ifstream inf(pname, ifstream::in);
|
|
if (!inf.good() && myrank == 0)
|
|
{
|
|
cout << "Can not open parameter file " << pname << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
|
|
for (int i = 1; inf.good(); i++)
|
|
{
|
|
inf.getline(pline, LEN);
|
|
str = pline;
|
|
|
|
int status = misc::parse_parts(str, sgrp, skey, sval, sind);
|
|
if (status == -1)
|
|
{
|
|
cout << "error reading parameter file " << pname << " in line " << i << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
else if (status == 0)
|
|
continue;
|
|
|
|
if (sgrp == "ABE")
|
|
{
|
|
if (skey == "checkrun")
|
|
checkrun = atoi(sval.c_str());
|
|
else if (skey == "checkfile")
|
|
strcpy(checkfilename, sval.c_str());
|
|
else if (skey == "ID Type")
|
|
ID_type = atoi(sval.c_str());
|
|
else if (skey == "Steps")
|
|
Steps = atoi(sval.c_str());
|
|
else if (skey == "StartTime")
|
|
StartTime = atof(sval.c_str());
|
|
else if (skey == "TotalTime")
|
|
TotalTime = atof(sval.c_str());
|
|
else if (skey == "DumpTime")
|
|
DumpTime = atof(sval.c_str());
|
|
else if (skey == "d2DumpTime")
|
|
d2DumpTime = atof(sval.c_str());
|
|
else if (skey == "CheckTime")
|
|
CheckTime = atof(sval.c_str());
|
|
else if (skey == "AnalysisTime")
|
|
AnasTime = atof(sval.c_str());
|
|
else if (skey == "Courant")
|
|
Courant = atof(sval.c_str());
|
|
else if (skey == "Symmetry")
|
|
Symmetry = atoi(sval.c_str());
|
|
else if (skey == "small dissipation")
|
|
numepss = atof(sval.c_str());
|
|
else if (skey == "big dissipation")
|
|
numepsb = atof(sval.c_str());
|
|
else if (skey == "shell dissipation")
|
|
numepsh = atof(sval.c_str());
|
|
else if (skey == "Analysis Level")
|
|
a_lev = atoi(sval.c_str());
|
|
else if (skey == "Max mode l")
|
|
maxl = atoi(sval.c_str());
|
|
else if (skey == "detector number")
|
|
decn = atoi(sval.c_str());
|
|
else if (skey == "farest detector position")
|
|
maxrex = atof(sval.c_str());
|
|
else if (skey == "detector distance")
|
|
drex = atof(sval.c_str());
|
|
else if (skey == "output dir")
|
|
out_dir = sval;
|
|
}
|
|
}
|
|
inf.close();
|
|
|
|
iter = parameters::str_par.find("output dir");
|
|
if (iter != parameters::str_par.end())
|
|
{
|
|
out_dir = iter->second;
|
|
}
|
|
else
|
|
{
|
|
parameters::str_par.insert(map<string, string>::value_type("output dir", out_dir));
|
|
}
|
|
}
|
|
|
|
if (myrank == 0)
|
|
{
|
|
string out_dir;
|
|
char filename[50];
|
|
map<string, string>::iterator iter;
|
|
iter = parameters::str_par.find("output dir");
|
|
if (iter != parameters::str_par.end())
|
|
{
|
|
out_dir = iter->second;
|
|
}
|
|
sprintf(filename, "%s/setting.par", out_dir.c_str());
|
|
ofstream setfile;
|
|
setfile.open(filename, ios::trunc);
|
|
|
|
if (!setfile.good())
|
|
{
|
|
char cmd[100];
|
|
// sprintf(cmd,"rm %s -f",out_dir.c_str());
|
|
// system(cmd);
|
|
sprintf(cmd, "mkdir %s", out_dir.c_str());
|
|
system(cmd);
|
|
|
|
setfile.open(filename, ios::trunc);
|
|
}
|
|
|
|
time_t tnow;
|
|
time(&tnow);
|
|
struct tm *loc_time;
|
|
loc_time = localtime(&tnow);
|
|
setfile << "# File created on " << asctime(loc_time);
|
|
setfile << "#" << endl;
|
|
// echo the micro definition in "microdef.fh"
|
|
setfile << "macro definition used in microdef.fh" << endl;
|
|
|
|
#if (tetradtype == 0)
|
|
setfile << "my own tetrad type for psi4 calculation" << endl;
|
|
#elif (tetradtype == 1)
|
|
setfile << "Lousto's tetrad type for psi4 calculation" << endl;
|
|
#elif (tetradtype == 2)
|
|
setfile << "Frans' tetrad type for psi4 calculation" << endl;
|
|
#else
|
|
setfile << "not recognized tetrad type" << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
#endif
|
|
|
|
#ifdef Cell
|
|
setfile << "Cell center numerical grid structure" << endl;
|
|
#endif
|
|
#ifdef Vertex
|
|
setfile << "Vertex center numerical grid structure" << endl;
|
|
#endif
|
|
|
|
setfile << " ghost zone = " << ghost_width << endl;
|
|
|
|
setfile << " buffer zone = " << buffer_width << endl;
|
|
|
|
#ifdef CPBC
|
|
setfile << "constraint preserving boundary condition is used" << endl;
|
|
setfile << " ghost zone for CPBC = " << CPBC_ghost_width << endl;
|
|
#endif
|
|
|
|
setfile << " Gauge type = " << GAUGE << endl;
|
|
|
|
#if (ABV == 0)
|
|
setfile << "using BSSN variable for constraint violation and psi4 calculation" << endl;
|
|
#elif (tetradtype == 1)
|
|
setfile << "using ADM variable for constraint violation and psi4 calculation" << endl;
|
|
#else
|
|
setfile << "not recognized ABV type" << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
#endif
|
|
|
|
// echo the micro definition in "microdef.h"
|
|
setfile << "macro definition used in microdef.h" << endl;
|
|
setfile << " Sommerfeld boundary type = " << SommerType << endl;
|
|
#ifdef GaussInt
|
|
setfile << "using Gauss integral in waveshell" << endl;
|
|
#else
|
|
setfile << "using usual integral in waveshell" << endl;
|
|
#endif
|
|
setfile << " ABE type = " << ABEtype << endl;
|
|
setfile << " ID type = " << ID_type << endl;
|
|
#ifdef With_AHF
|
|
setfile << "Apparent Horizon Finder is turned on" << endl;
|
|
#endif
|
|
setfile << " Psi4 calculation type = " << Psi4type << endl;
|
|
#ifdef Point_Psi4
|
|
setfile << "Using point Psi4 calculation method" << endl;
|
|
#endif
|
|
setfile << " RestrictProlong time type = " << RPS << endl;
|
|
setfile << " RestrictProlong scheme type = " << RPB << endl;
|
|
setfile << "Enforce algebra constraint type = " << AGM << endl;
|
|
setfile << "Analysis and PBH treat type = " << MAPBH << endl;
|
|
setfile << " mesh level parallel type = " << PSTR << endl;
|
|
setfile << " regrid type = " << REGLEV << endl;
|
|
|
|
setfile << " dim = " << dim << endl;
|
|
setfile << " buffer_width = " << buffer_width << endl;
|
|
setfile << " SC_width = " << SC_width << endl;
|
|
setfile << " CS_width = " << CS_width << endl;
|
|
|
|
setfile.close();
|
|
}
|
|
|
|
// echo parameters
|
|
if (myrank == 0)
|
|
{
|
|
cout << endl;
|
|
cout << " /////////////////////////////////////////////////////////////// " << endl;
|
|
cout << " AMSS-NCKU Begin !!! " << endl;
|
|
cout << " /////////////////////////////////////////////////////////////// " << endl;
|
|
cout << endl;
|
|
|
|
if (checkrun)
|
|
cout << " checked run" << endl;
|
|
else
|
|
cout << " new run" << endl;
|
|
|
|
cout << " simulation with cpu numbers = " << nprocs << endl;
|
|
cout << " simulation time = (" << StartTime << ", " << TotalTime << ")" << endl;
|
|
cout << " simulation steps for this run = " << Steps << endl;
|
|
cout << " Courant number = " << Courant << endl;
|
|
|
|
switch (ID_type)
|
|
{
|
|
case -3:
|
|
cout << " Initial Data Type: Analytical NBH (Cao's Formula)" << endl;
|
|
break;
|
|
case -2:
|
|
cout << " Initial Data Type: Analytical Kerr-Schild" << endl;
|
|
break;
|
|
case -1:
|
|
cout << " Initial Data Type: Analytical NBH (Lousto's Formula)" << endl;
|
|
break;
|
|
case 0:
|
|
cout << " Initial Data Type: Numerical Ansorg TwoPuncture" << endl;
|
|
break;
|
|
case 1:
|
|
cout << " Initial Data Type: Numerical Pablo" << endl;
|
|
break;
|
|
default:
|
|
cout << " OOOOps, not supported Initial Data setting!" << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
|
|
switch (Symmetry)
|
|
{
|
|
case 0:
|
|
cout << " Symmetry setting: No_Symmetry" << endl;
|
|
break;
|
|
case 1:
|
|
cout << " Symmetry setting: Equatorial" << endl;
|
|
break;
|
|
case 2:
|
|
cout << " Symmetry setting: Octant" << endl;
|
|
break;
|
|
default:
|
|
cout << " OOOOps, not supported Symmetry setting!" << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
|
|
cout << " Courant = " << Courant << endl;
|
|
cout << " artificial dissipation for shell patches = " << numepsh << endl;
|
|
cout << " artificial dissipation for fixed levels = " << numepsb << endl;
|
|
cout << " artificial dissipation for moving levels = " << numepss << endl;
|
|
cout << " Dumpt Time = " << DumpTime << endl;
|
|
cout << " Check Time = " << CheckTime << endl;
|
|
cout << " Analysis Time = " << AnasTime << endl;
|
|
cout << " Analysis level = " << a_lev << endl;
|
|
cout << " checkfile = " << checkfilename << endl;
|
|
|
|
switch (ghost_width)
|
|
{
|
|
case 2:
|
|
cout << " second order finite difference is used" << endl;
|
|
break;
|
|
case 3:
|
|
cout << " fourth order finite difference is used" << endl;
|
|
break;
|
|
case 4:
|
|
cout << " sixth order finite difference is used" << endl;
|
|
break;
|
|
case 5:
|
|
cout << " eighth order finite difference is used" << endl;
|
|
break;
|
|
default:
|
|
cout << " Why are you using ghost width = " << ghost_width << endl;
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
|
|
cout << "///////////////////////////////////////////////////////////////" << endl;
|
|
}
|
|
|
|
//===========================the computation body====================================================
|
|
|
|
bssn_class *ADM;
|
|
|
|
#if (ABEtype == 0)
|
|
ADM = new bssn_class(Courant, StartTime, TotalTime, DumpTime, d2DumpTime, CheckTime, AnasTime,
|
|
Symmetry, checkrun, checkfilename, numepss, numepsb, numepsh,
|
|
a_lev, maxl, decn, maxrex, drex);
|
|
#elif (ABEtype == 1)
|
|
ADM = new bssnEScalar_class(Courant, StartTime, TotalTime, DumpTime, d2DumpTime, CheckTime, AnasTime,
|
|
Symmetry, checkrun, checkfilename, numepss, numepsb, numepsh,
|
|
a_lev, maxl, decn, maxrex, drex);
|
|
#elif (ABEtype == 2)
|
|
ADM = new Z4c_class(Courant, StartTime, TotalTime, DumpTime, d2DumpTime, CheckTime, AnasTime,
|
|
Symmetry, checkrun, checkfilename, numepss, numepsb, numepsh,
|
|
a_lev, maxl, decn, maxrex, drex);
|
|
#elif (ABEtype == 3)
|
|
ADM = new bssnEM_class(Courant, StartTime, TotalTime, DumpTime, d2DumpTime, CheckTime, AnasTime,
|
|
Symmetry, checkrun, checkfilename, numepss, numepsb, numepsh,
|
|
a_lev, maxl, decn, maxrex, drex);
|
|
#endif
|
|
|
|
ADM->Initialize();
|
|
// ADM->testRestrict();
|
|
// ADM->testOutBd();
|
|
|
|
// set up initial data
|
|
|
|
// old code manually
|
|
/*
|
|
#if (ABEtype == 0)
|
|
// set up initial data with analytical formula
|
|
// ADM->Setup_Initial_Data();
|
|
ADM->Read_Ansorg();
|
|
#elif (ABEtype == 1)
|
|
// ADM->Read_Pablo();
|
|
ADM->Read_Ansorg();
|
|
#elif (ABEtype == 2)
|
|
ADM->Read_Ansorg();
|
|
// ADM->Setup_KerrSchild();
|
|
#elif (ABEtype == 3)
|
|
ADM->Setup_Initial_Data();
|
|
// ADM->Read_Ansorg();
|
|
#endif
|
|
*/
|
|
|
|
// new code Xiao Qu
|
|
switch (ID_type)
|
|
{
|
|
case (-3):
|
|
// set up initial data with Cao's analytical formula
|
|
ADM->Setup_Initial_Data_Cao();
|
|
break;
|
|
case (-2):
|
|
// set up initial data with KerrSchild analytical formula
|
|
ADM->Setup_KerrSchild();
|
|
break;
|
|
case (-1):
|
|
// set up initial data with Lousto's analytical formula
|
|
ADM->Setup_Initial_Data_Lousto();
|
|
break;
|
|
case (0):
|
|
// set up initial data with Ansorg TwoPuncture Solver
|
|
ADM->Read_Ansorg();
|
|
break;
|
|
case (1):
|
|
// set up initial data with Pablo's Olliptic Solver
|
|
ADM->Read_Pablo();
|
|
// ADM->Write_Pablo();
|
|
break;
|
|
default:
|
|
if (myrank == 0)
|
|
{
|
|
cout << "not recognized ABE::InitialDataType = " << ID_type << endl;
|
|
}
|
|
MPI_Abort(MPI_COMM_WORLD, 1);
|
|
}
|
|
|
|
End_clock = MPI_Wtime();
|
|
if (myrank == 0)
|
|
{
|
|
cout << endl;
|
|
cout << " Before Evolve, it takes " << MPI_Wtime() - Begin_clock << " seconds" << endl;
|
|
cout << endl;
|
|
}
|
|
|
|
ADM->Evolve(Steps);
|
|
|
|
if (myrank == 0)
|
|
{
|
|
cout << endl;
|
|
cout << " Total Evolve Time: " << MPI_Wtime() - End_clock << " seconds!" << endl;
|
|
cout << " Total Running Time: " << MPI_Wtime() - Begin_clock << " seconds!" << endl;
|
|
cout << endl;
|
|
}
|
|
|
|
delete ADM;
|
|
|
|
//=======================caculation done=============================================================
|
|
|
|
if (myrank == 0)
|
|
{
|
|
cout << endl;
|
|
cout << " =============================================================== " << endl;
|
|
cout << " Simulation is successfully done!! " << endl;
|
|
cout << " =============================================================== " << endl;
|
|
cout << endl;
|
|
cout << " This run used " << MPI_Wtime() - Begin_clock << " seconds! " << endl;
|
|
cout << endl;
|
|
}
|
|
|
|
MPI_Finalize();
|
|
|
|
exit(0);
|
|
}
|
|
|
|
//===================================================================================================
|
|
//===================================================================================================
|