Author: Chris Douglas (@cmdoug) christopher.douglas@duke.edu
ff-mpirun -np 4 basedeflate.md -mi <FILEIN> -fo <FILEOUT>ff-mpirun -np 4 basedeflate.md -fi <FILEIN> -fo <FILEOUT>ff-mpirun -np 4 basedeflate.md -mi <MESHIN> -fi <FILEIN> -fo <FILEOUT>ff-mpirun -np 4 basedeflate.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, basecontinue.md, modecompute.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", ""); // output mesh without extension
string filein = getARGV("-fi", ""); // input file with extension
string deflatefilein = getARGV("-fi2", ""); // input file with extension
string fileout = getARGV("-fo", ""); // output file without extension
string sneslinesearchtype = getARGV("-snes_linesearch_type", "basic");
real TGV = getARGV("-tgv", -1.);
real defp = getARGV("-defp", 1.0); // deflation root order
real defa = getARGV("-defa", 1.0); // deflation zero bias
int ndeflate = getARGV("-ndeflate", 3);
real noiseamp = getARGV("-noise", 0.0);
string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension
string deflatefileroot, deflatefileext = parsefilename(deflatefilein, deflatefileroot); //extract file name and extension
parsefilename(fileout, fileout); // trim extension from output file, if given
if(fileext == "mode" || fileext == "resp" || fileext == "rslv" || fileext == "tdls" || fileext == "floq"){
filein = readbasename(workdir + filein);
fileext = parsefilename(filein, fileroot);
}
if(filein != "" && meshin == "") meshin = readmeshname(workdir + filein); // get mesh file
string meshroot, meshext = parsefilename(meshin, meshroot);
parsefilename(meshout, meshout); // trim extension from output mesh, if given
// Load mesh
Th = readmeshN(workdir + meshin);
Thg = Th;
// Partition mesh across processors
DmeshCreate(Th);
restu = restrict(XMh, XMhg, n2o);
// Make finite element basis
XMh defu(ub), defu(um), defu(um2), defu(um3);
if(fileext == "base") {
ub[] = loadbase(fileroot, meshin);
}
else if(fileext == "fold") {
real[string] alpha;
real beta;
real[int] qm, qma;
ub[] = loadfold(fileroot, meshin, qm, qma, alpha, beta);
}
else if(fileext == "hopf") {
real omega;
complex[string] alpha;
complex beta;
complex[int] qm, qma;
ub[] = loadhopf(fileroot, meshin, qm, qma, sym, omega, alpha, beta);
}
else if(fileext == "foho") {
real omega;
complex[string] alpha1;
real[string] alpha2;
complex beta1, gamma12, gamma13;
real beta22, beta23, gamma22, gamma23;
complex[int] q1m, q1ma;
real[int] q2m, q2ma;
ub[] = loadfoho(fileroot, meshin, q1m, q1ma, q2m, q2ma, sym, omega, alpha1, alpha2, beta1, beta22, beta23, gamma12, gamma13, gamma22, gamma23);
}
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, q1ma, q2m, q2ma;
ub[] = loadhoho(fileroot, meshin, q1m, q1ma, q2m, q2ma, sym1, sym2, omega1, omega2, alpha1, alpha2, beta1, beta2, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23);
}
else if(fileext == "tdns") {
real time;
ub[] = loadtdns(fileroot, meshin, time);
}
else if(fileext == "porb") {
int Nh=1;
real omega;
complex[int, int] qh(ub[].n, Nh);
ub[] = loadporb(fileroot, meshin, qh, sym, omega, Nh);
}
else {
setparams(paramnames,params);
defu(ub) = InitialConditions;
}
if(deflatefileext == "base") {
um[] = loadbase(deflatefileroot, meshin);
}
// Create distributed Mat
Mat J;
createMatu(Th, J, Pk);
// MESH ADAPTATION
bool adapt = false;
if(meshout == "") meshout = meshin; // if no adaptation
else { // if output meshfile is given, adapt mesh
adapt = true;
meshout = meshout + "." + meshext;
real[int] q;
ChangeNumbering(J, ub[], q);
ChangeNumbering(J, ub[], q, inverse = true);
ChangeNumbering(J, um[], q);
ChangeNumbering(J, um[], q, inverse = true);
XMhg defu(ubG), defu(umG), defu(tempu); // create private global FE functions
tempu[](restu) = ub[]; // populate local portion of global soln
mpiAllReduce(tempu[], ubG[], mpiCommWorld, mpiSUM); //aggregate local solns into global soln
tempu[](restu) = um[]; // populate local portion of 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(ubG), 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(ubG), 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(ubG) = defu(ubG); //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;
defu(ub) = initu(0.0); // set local values to zero
defu(um) = 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[] = ubG[](restu); //restrict global solution to each local mesh
um[] = umG[](restu); //restrict global solution to each local mesh
}
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
// Initialize
real[int] q, q0(J.n);
real[int, int] qdeflate(J.n, (deflatefilein != ""));
if (deflatefilein != ""){
ChangeNumbering(J, um[], q);
qdeflate(:, 0) = q;
}
ChangeNumbering(J, ub[], q);
if (noiseamp > 0.0) for [i, ai : q] ai += noiseamp*randreal1();
q0 = q;
// Function to build residual operator in PETSc numbering
func real[int] funcR(real[int]& qPETSc) {
ChangeNumbering(J, ub[], qPETSc, inverse = true, exchange = true);
real[int] RPETSc, R = vR(0, XMh, tgv = TGV);
ChangeNumbering(J, R, RPETSc);
return RPETSc;
}
// Function to build Jacobian operator in PETSc numbering
func int funcJ(real[int]& qPETSc) {
q = qPETSc;
ChangeNumbering(J, ub[], qPETSc, inverse = true, exchange = true);
J = vJ(XMh, XMh, tgv = TGV);
return 0;
}
func real[int] funcshell(PetscScalar[int]& RPETSc) {
PetscScalar[int] dq(RPETSc.n);
KSPSolve(J, RPETSc, dq);
real factor = 1.0;
for (int deflatej = 0; deflatej < qdeflate.m; deflatej++) {
PetscScalar[int] qdiff = q - qdeflate(:, deflatej);
real innerprod, normsq, loc = (qdiff'*qdiff);
mpiAllReduce(loc, normsq, mpiCommWorld, mpiSUM); //aggregate local solns into global soln
loc = (qdiff'*dq);
mpiAllReduce(loc, innerprod, mpiCommWorld, mpiSUM); //aggregate local solns into global soln
factor *= (1.0 + defa*normsq^(0.5*defp))/(1.0 + defa*normsq^(0.5*defp) - defp/normsq*innerprod);
}
return dq *= factor;
}
// set up Mat parameters
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO sparams = KSPparams);
Mat deflatedJ = [[J]];
set(deflatedJ, precon = funcshell, sparams = "-pc_type shell");
int ret;
ndeflate += (deflatefilein != "");
for (int deflatei = (deflatefilein != ""); deflatei < ndeflate; deflatei++) {
// solve nonlinear problem with SNES
q = q0;
SNESSolve(deflatedJ, funcJ, funcR, q, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -snes_monitor -snes_converged_reason -options_left no");
if(ret > 0) { // Save solution if solver converged and output file is given
if(mpirank == 0 && adapt) { // Save adapted mesh
cout << " Saving adapted mesh '" + meshout + "' in '" + workdir + "'." << endl;
savemesh(Thg, workdir + meshout);
}
qdeflate.resize(J.n, qdeflate.m + 1);
qdeflate(:, deflatei) = q;
ChangeNumbering(J, ub[], q, inverse = true);
savebase(fileout + "-" + deflatei, "", meshout, true, true);
adapt = false;
}
else break;
}