Author: Chris Douglas (@cmdoug) christopher.douglas@duke.edu
ff-mpirun -np 4 tdnscompute.md -fi <FILEIN> -fo <FILEOUT>ff-mpirun -np 4 tdnscompute.md -fi <FILEIN> -fo <FILEOUT> -mo <MESHOUT>NOTE: This file should not be changed unless you know what you're doing.
SEE ALSO: basecompute.md, porbcompute.md, modecompute.md, rslvcompute.md, respcompute.md, tdlscompute.md
load "iovtk"
load "PETSc"
include "settings.idp"
include "macros_bifbox.idp"
// arguments
string meshin = getARGV("-mi", ""); // input meshfile with extension
string meshout = getARGV("-mo", "");
string filein = getARGV("-fi", "");
string basefilein = getARGV("-bfi", "");
string fileout = getARGV("-fo", filein);
bool adj = getARGV("-adj", 0);
real amplitude = abs(getARGV("-amp", 0.0));
real phase = getARGV("-phase", 0.0);
int select = getARGV("-select", 1);
int count = getARGV("-count", 0);
int savecount = getARGV("-scount", 1);
int maxcount = getARGV("-maxcount", 100);
string tstype = getARGV("-ts_type", "beuler");
int tsmaxsnesfailures = getARGV("-ts_max_snes_failures", -1);
real tstimestep = getARGV("-ts_time_step", 0.01);
string tsadapttype = getARGV("-ts_adapt_type", "none");
string sneslinesearchtype = getARGV("-snes_linesearch_type", "none");
real TGV = getARGV("-tgv", -1);
// Load mesh, make FE basis
string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension
parsefilename(fileout, fileout); // trim extension from output file, if given
string basefileroot, basefileext = parsefilename(basefilein, basefileroot);
if((fileext == "mode" || fileext == "resp" || fileext == "rslv" || fileext == "tdls" || fileext == "floq") && basefilein == ""){
basefilein = readbasename(workdir + filein);
basefileext = parsefilename(basefilein, basefileroot);
}
if(meshin == "") meshin = readmeshname(workdir + filein); // get mesh file
string meshroot, meshext = parsefilename(meshin, meshroot);
parsefilename(meshout, meshroot); // trim extension from output mesh, if given
if(count > 0) {
fileroot = fileroot(0:fileroot.rfind("_" + count)-1); // get file root
meshroot = meshroot(0:meshroot.rfind("_" + count)-1); // get file root
}
Th = readmeshN(workdir + meshin);
Thg = Th;
DmeshCreate(Th);
restu = restrict(XMh, XMhg, n2o);
XMh defu(ub), defu(um), defu(um2), defu(um3);
real time = 0.0;
if (count == 0){
if(fileext == "base") {
ub[] = loadbase(fileroot, meshin);
}
else if(fileext == "mode") {
complex eigenvalue;
complex[int] qm(um[].n);
qm = loadmode(fileroot, meshin, sym, eigenvalue);
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(fileext == "resp") {
real omega;
complex[int] qm(um[].n);
qm = loadresp(fileroot, meshin, sym, omega);
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(fileext == "rslv") {
real omega, gain;
complex[int] fm, qm(um[].n);
qm = loadrslv(fileroot, meshin, fm, sym, omega, gain);
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(fileext == "tdls") {
real time;
complex[int] qm(um[].n);
qm = loadtdls(fileroot, meshin, sym, time);
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(fileext == "fold") {
real[string] alpha;
real beta;
real[int] qm(um[].n), qma(um[].n);
ub[] = loadfold(fileroot, meshin, qm, qma, alpha, beta);
if(adj) qm = qma;
um[] = qm;
}
else if(fileext == "cusp") {
real[string] alpha, alphaR;
real beta;
real[int] qm(um[].n), qma(um[].n);
ub[] = loadfold(fileroot, meshin, qm, qma, alpha, alphaR, beta);
if(adj) qm = qma;
um[] = qm;
}
else if(fileext == "hopf") {
real omega;
complex[string] alpha;
complex beta;
complex[int] qm(um[].n), qma(um[].n);
ub[] = loadhopf(fileroot, meshin, qm, qma, sym, omega, alpha, beta);
if(adj) qm = qma;
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(fileext == "bota") {
real[string] alpha1, alpha2;
real beta1, beta2, beta3, beta4;
real[int] qm(um[].n), qma(um[].n);
ub[] = loadbota(fileroot, meshin, qm, qma, alpha1, alpha2, beta1, beta2, beta3, beta4);
if(adj) qm = qma;
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(fileext == "foho") {
real omega;
complex[string] alpha1;
real[string] alpha2;
complex beta1, gamma12, gamma13;
real beta22, beta23, gamma22, gamma23;
complex[int] q1m(um[].n), q1ma(um[].n);
real[int] q2m(um[].n), q2ma(um[].n);
ub[] = loadfoho(fileroot, meshin, q1m, q1ma, q2m, q2ma, sym, omega, alpha1, alpha2, beta1, beta22, beta23, gamma12, gamma13, gamma22, gamma23);
if (select == 1){
if(adj) q1m = q1ma;
q1m *= exp(1i*phase);
um[] = 2.0*q1m.re;
}
else if(select == 2){
if(adj) q2m = q2ma;
um[] = q2m;
}
}
else if(fileext == "hoho") {
real[int] sym1(sym.n), sym2(sym.n);
real omega1, omega2;
complex[string] alpha1, alpha2;
complex beta1, beta2, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23;
complex[int] q1m(um[].n), q1ma(um[].n), q2m(um[].n), q2ma(um[].n);
ub[] = loadhoho(fileroot, meshin, q1m, q1ma, q2m, q2ma, sym1, sym2, omega1, omega2, alpha1, alpha2, beta1, beta2, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23);
if(select == 1 && adj) q1m = q1ma;
else if (select == 2 && !adj) q1m = q2m;
else if (select == 2 && adj) q1m = q2ma;
q1m *= exp(1i*phase);
um[] = 2.0*q1m.re;
}
else if(fileext == "tdns") {
real time;
ub[] = loadtdns(fileroot, meshin, time);
}
else if(fileext == "porb") {
int Nh;
real omega;
complex[int] qm(um[].n);
complex[int, int] qh(um[].n, Nh);
ub[] = loadporb(fileroot, meshin, qh, sym, omega, Nh);
for (int nh = 0; nh < Nh; nh++) {
qm = exp(1i*(nh+1)*phase)*qh(:,nh);
um[] += (2.0*qm.re);
}
}
else if(fileext == "floq") {
int Nh=1;
complex[int, int] qh(um[].n, 2);
complex eigenvalue;
real[int] symtemp(sym.n);
complex[int] qm(um[].n);
real omega;
qm = loadfloq(fileroot, meshin, qh, sym, eigenvalue, symtemp, omega, Nh);
qm *= exp(1i*phase);
um[] = 2.0*qm.re;
}
else if(filein == "" && basefilein == "") {
setparams(paramnames, params);
defu(ub) = InitialConditions;
}
if (basefileext == "base") {
ub[] = loadbase(basefileroot, meshin);
}
else if(basefileext == "fold") {
real[string] alpha;
real beta;
real[int] qm, qma;
ub[] = loadfold(basefileroot, meshin, qm, qma, alpha, beta);
}
else if(basefileext == "cusp") {
real[string] alpha, alphaR;
real beta;
real[int] qm, qma;
ub[] = loadcusp(basefileroot, meshin, qm, qma, alpha, alphaR, beta);
}
else if(basefileext == "hopf") {
real omega;
complex[string] alpha;
complex beta;
complex[int] qm, qma;
ub[] = loadhopf(basefileroot, meshin, qm, qma, sym, omega, alpha, beta);
}
else if(basefileext == "bota") {
real[string] alpha1, alpha2;
real beta1, beta2, beta3, beta4;
real[int] qm, qma;
ub[] = loadbota(basefileroot, meshin, qm, qma, alpha1, alpha2, beta1, beta2, beta3, beta4);
}
else if(basefileext == "foho") {
real omega;
complex[string] alpha1;
complex beta1, gamma12, gamma13;
real[string] alpha2;
real beta22, beta23, gamma22, gamma23;
complex[int] q1m, q1ma;
real[int] q2m, q2ma;
ub[] = loadfoho(basefileroot, meshin, q1m, q1ma, q2m, q2ma, sym, omega, alpha1, alpha2, beta1, beta22, beta23, gamma12, gamma13, gamma22, gamma23);
}
else if(basefileext == "hoho") {
real[int] sym2(sym.n);
real omega1, omega2;
complex[string] alpha1, alpha2;
complex beta1, beta2, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23;
complex[int] q1m, q1ma, q2m, q2ma;
ub[] = loadhoho(basefileroot, meshin, q1m, q1ma, q2m, q2ma, sym, sym2, omega1, omega2, alpha1, alpha2, beta1, beta2, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23);
}
else if(basefileext == "tdns") {
real time;
ub[] = loadtdns(basefileroot, meshin, time);
}
else if(basefileext == "porb") {
int Nh=1;
real omega;
complex[int, int] qh(um[].n, Nh);
ub[] = loadporb(basefileroot, meshin, qh, sym, omega, Nh);
}
um[] *= amplitude;
ub[] += um[];
}
else {
ub[] = loadtdns(fileroot + "_" + count, meshin, time);
}
// Create distributed Mat
Mat J;
createMatu(Th, J, Pk);
sym = 0;
real[int] ik(sym.n), ik2(sym.n), ik3(sym.n);
real iomega = 0.0, iomega2 = 0.0, iomega3 = 0.0;
include "eqns.idp" // load equations
bool adapt = false;
if(meshout == "") meshout = meshin; // if no adaptation
else adapt = true; // if output meshfile is given, adapt mesh
int adaptflag;
// initialize
real[int] q;
ChangeNumbering(J, ub[], q);
ChangeNumbering(J, ub[], q, inverse = true);
if (count == 0) savetdns(fileout + "_0", (savecount > 0 ? fileout : ""), meshin, filein, time, true, true);
real[int] FPETSc(J.n), qp = q;
// Function to build residual "F(t,u,\dot{u})" operator in PETSc numbering
func real[int] funcF(real t, real[int]& qPETSc, real[int]& qdotPETSc) {
ChangeNumbering(J, ub[], qPETSc, inverse = true, exchange = true);
ChangeNumbering(J, um[], qdotPETSc, inverse = true, exchange = true);
um2[] = vR(0, XMh, tgv = TGV);
um3[] = vM(0, XMh, tgv = -10);
um2[] += um3[];
ChangeNumbering(J, um2[], FPETSc);
return FPETSc;
}
// Function to build Jacobian operator in PETSc numbering
func int funcJ(real t, real[int]& qPETSc, real[int]& qdotPETSc, real a) {
ChangeNumbering(J, ub[], qPETSc, inverse = true, exchange = true);
iomega = a;
J = vJ(XMh, XMh, tgv = TGV);
return 0;
}
// Function to monitor solution progress
func int funcMon(int s, real t, real[int]& in) {
if(s > 0){
tstimestep = t - time;
count++;
if(mpirank == 0) cout << " " << ((adapt && adaptflag) ? "A" : "") << count + ":\tdt = " + tstimestep + ",\ttime = " + t << endl;
if( adapt ? (adaptflag ? false : (count % savecount != 0)) : true){
qp = in;
ChangeNumbering(J, ub[], in, inverse = true);
savetdns(fileout + "_" + count, ((savecount > 0) ? fileout : ""), meshout, filein, t, (count % savecount == 0), true)
}
}
time = t;
return 0;
}
// set up Mat parameters
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO sparams = KSPparams);
while (count < maxcount){
int ret;
adaptflag = 0;
// solve nonlinear problem with TS
TSSolve(J, funcJ, funcF, q, monitor = funcMon, reason = ret, sparams = " -ts_init_time " + time
+ " -ts_type " + tstype
+ " -ts_max_snes_failures " + tsmaxsnesfailures
+ " -ts_time_step " + tstimestep
+ " -ts_max_steps " + (adapt ? min(maxcount-count,savecount-(count % (savecount))) : (maxcount-count))
+ " -ts_adapt_type " + tsadapttype
+ " -snes_linesearch_type " + sneslinesearchtype
+ " -options_left no "
);
if (ret > 0 && adapt && (count % savecount == 0)) { // if output meshfile is given and solver converged, adapt mesh
ChangeNumbering(J, ub[], q, inverse = true);
ChangeNumbering(J, um[], qp, inverse = true);
meshout = meshroot + "_" + count + "." + meshext;
XMhg defu(uG), defu(umG), defu(tempu); // create private global FE functions
tempu[](restu) = ub[]; // populate local portion of global soln
mpiAllReduce(tempu[], uG[], mpiCommWorld, mpiSUM); //aggregate local solns into global soln
tempu[](restu) = um[]; // populate local portion of prev. global soln
mpiAllReduce(tempu[], umG[], mpiCommWorld, mpiSUM); //aggregate local solns into global soln
if(mpirank == 0) { // Perform mesh adaptation (serially) on processor 0
IFMACRO(dimension,2)
Thg = adaptmesh(Thg, adaptu(uG), adaptmeshoptions);
ENDIFMACRO
IFMACRO(dimension,3)
//NOTE: 3D mesh adaptation is still under development.
load "mshmet"
load "mmg"
real anisomax = getARGV("-anisomax",1.0);
real[int] met = mshmet(Thg, adaptu(uG), normalization = getARGV("-normalization",1), aniso = bool(anisomax > 1.0), hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), err = getARGV("-err", 1.0e-2));
if(anisomax > 1.0) {
load "aniso"
boundaniso(6, met, anisomax);
}
Thg = mmg3d(Thg, metric = met, hmin = getARGV("-hmin", 1.0e-6), hmax = getARGV("-hmax", 1.0e+2), hgrad = -1, verbose = verbosity-(verbosity==0));
ENDIFMACRO
}
broadcast(processor(0), Thg); // broadcast global mesh to all processors
defu(uG) = defu(uG); //interpolate global solution from old mesh to new mesh
defu(umG) = defu(umG); //interpolate global solution from old mesh to new mesh
Th = Thg; //Reinitialize local mesh with global mesh
Mat Adapt;
createMatu(Th, Adapt, Pk); // Partition new mesh and update the PETSc numbering
J = Adapt;
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO sparams = KSPparams);
defu(ub) = initu(0.0); // set local values to zero
defu(um) = initu(0.0); // set local values to zero
defu(um2) = initu(0.0); // set local values to zero
defu(um3) = initu(0.0); // set local values to zero
restu.resize(ub[].n); // Change size of restriction operator
restu = restrict(XMh, XMhg, n2o); // Compute new restriction from global mesh to local mesh
ub[] = uG[](restu); //restrict global solution to each local mesh
um[] = umG[](restu); //restrict global solution to each local mesh
Mat Adapt0(J);
q.resize(J.n); // Resize state vector
qp.resize(J.n); // Resize state vector
FPETSc.resize(J.n);
ChangeNumbering(J, um[], q);
adaptflag = 1;
count--;
time -= tstimestep;
TSSolve(J, funcJ, funcF, q, monitor = funcMon, reason = ret, sparams = " -ts_init_time " + time
+ " -ts_type " + tstype
+ " -ts_max_snes_failures " + tsmaxsnesfailures
+ " -ts_time_step " + tstimestep
+ " -snes_linesearch_type " + sneslinesearchtype
+ " -ts_max_steps 1 -ts_adapt_type none -options_left no "
);
assert(ret > 0);
if(mpirank == 0) { // Save adapted mesh
cout << " Saving adapted mesh '" + meshout + "' in '" + workdir + "'." << endl;
savemesh(Thg, workdir + meshout);
}
ChangeNumbering(J, ub[], q, inverse = true);
savetdns(fileout + "_" + count, (savecount > 0 ? fileout : ""), meshout, filein, time, (count % savecount == 0), true);
qp = q;
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
}
}