Author: Chris Douglas (@cmdoug) christopher.douglas@duke.edu
ff-mpirun -np 4 basecontinue.md -fi <FILEIN> -param <PARAM>ff-mpirun -np 4 basecontinue.md -fi <FILEIN> -param <PARAM> -mo <MESHOUT>NOTE: This file should not be changed unless you know what you're doing.
SEE ALSO: basecompute.md, basedeflate.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", "");
string filein = getARGV("-fi", "");
string branchfilein = getARGV("-fi2", "");
string fileout = getARGV("-fo", filein);
int contorder = getARGV("-contorder", 1);
int count = getARGV("-count", 0);
int savecount = getARGV("-scount", 1);
int maxcount = getARGV("-maxcount", 100);
real h0 = getARGV("-h0", 1.0);
string param = getARGV("-param", "");
string adaptto = getARGV("-adaptto", "b");
real fmax = getARGV("-fmax", 2.0);
real kappamax = getARGV("-kmax", 0.5);
real deltamax = getARGV("-dmax", 4.0);
real anglemax = getARGV("-amax", 30.)*pi/180.0;
real monotone = getARGV("-mono", 1.0);
real eps = getARGV("-eps", 1.0e-7);
bool stricttangent = bool(getARGV("-stricttangent", 1));
int snesmaxit = getARGV("-snes_max_it", 10);
string sneslinesearchtype = getARGV("-snes_linesearch_type", "none");
int refactor = getARGV("-refact", snesmaxit);
real paramtarget = getARGV("-paramtarget", 1.0);
real TGV = getARGV("-tgv", -1);
bool stopflag = false;
bool forcesave = false;
// Load mesh, make FE basis
string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension
string branchfileroot, branchfileext = parsefilename(branchfilein, branchfileroot); //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(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);
if(count == 0) {
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;
complex beta1, gamma12, gamma13;
real[string] alpha2;
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 (filein == "") {
setparams(paramnames,params);
defu(ub) = InitialConditions;
}
savebase(filein, (savecount > 0 ? fileout : ""), meshin, false, false);
}
else {
ub[] = loadbase(fileroot + "_" + count, meshin);
}
if (branchfilein != ""){
if(branchfileext == "base") {
um[] = loadbase(branchfileroot, meshin);
}
else if(branchfileext == "mode") {
complex eigenvalue;
complex[int] qm = loadmode(fileroot, meshin, sym, eigenvalue);
um[] = qm.re;
}
}
real paramval = getparam(param);
real paramdiff = paramval - paramtarget;
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"
bool adaptflag, adapt = false;
if(meshout != "") adapt = true; // if output meshfile is given, adapt mesh
meshout = meshin; // if no adaptation
// Build bordered block matrix from only Mat components
Mat JlPM(J.n, mpirank == 0 ? 1 : 0), yqPM(J.n, mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
Mat Ja = [[J, JlPM], [yqPM', -1.0]]; // make dummy Jacobian
real[int] R(ub[].n), yqP(J.n), yqP0(J.n);
int ret, it = 0;
real f, kappa, cosalpha, res, delta, maxdelta;
// FUNCTIONS
func PetscScalar[int] funcRa(PetscScalar[int]& qa) {
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true); // PETSc to FreeFEM
if(mpirank == 0) paramval = qa(Ja.n-1); // Extract parameter value from state vector on proc 0
broadcast(processor(0), paramval);
updateparam(param, paramval);
R = vR(0, XMh, tgv = TGV);
PetscScalar[int] Ra;
ChangeNumbering(J, R, Ra); // FreeFEM to PETSc
Ra.resize(Ja.n); // Append 0 to residual vector on proc 0
if(mpirank == 0) Ra(Ja.n-1) = 0.0;
return Ra;
}
func int funcJa(PetscScalar[int]& qa) {
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true); // PETSc to FreeFEM
if(mpirank == 0) paramval = qa(Ja.n-1); // Extract parameter value from state vector on proc 0
broadcast(processor(0), paramval);
updateparam(param, paramval);
if (it <= refactor) J = vJ(XMh, XMh, tgv = TGV);
if (contorder > 0){
updateparam(param, paramval + eps);
um2[] = vR(0, XMh, tgv = TGV);
updateparam(param, paramval);
um2[] -= R;
um2[] /= eps;
ChangeNumbering(J, um2[], yqP); // FreeFEM to PETSc
matrix<PetscScalar> tempPms = [[yqP]]; // dense array to sparse matrix
ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
KSPSolve(J, yqP, yqP); // compute tangent vector in PETSc numbering
tempPms = [[yqP]]; // dense array to sparse matrix
ChangeOperator(yqPM, tempPms, parent = Ja); // send to Mat
}
return 0;
}
ConvergenceCheck(param + " = " + paramval);
// set up Mat parameters
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(Ja, sparams = "-ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition full"
+ " -prefix_push fieldsplit_1_ -ksp_type preonly -pc_type redundant -redundant_pc_type lu -prefix_pop"
+ " -prefix_push fieldsplit_0_ " + KSPparams + " -prefix_pop", setup = 1);
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO prefix = "fieldsplit_0_", parent = Ja);
// PREDICTOR
real[int] qa(Ja.n), qa0;
ChangeNumbering(J, ub[], qa0);
qa0.resize(Ja.n);
if(mpirank == 0) qa0(Ja.n-1) = paramval;
if (contorder > 0) {
R = vR(0, XMh, tgv = TGV);
funcJa(qa0);
if (branchfilein != "") {
ChangeNumbering(J, um[], yqP);
if (contorder > 1) J = vJ(XMh, XMh, tgv = TGV);
}
}
else {
matrix<PetscScalar> tempPms = [[yqP]];
ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
ChangeOperator(yqPM, tempPms, parent = Ja); // send to Mat
}
yqP0 = yqP;
while (!stopflag){
qa = qa0;
if (contorder == 0 && mpirank == 0) qa(Ja.n-1) += h0;
else if (contorder > 0) {
real h, hl = (yqP0'*yqP0);
mpiAllReduce(hl, h, mpiCommWorld, mpiSUM);
h = h0/sqrt(h + 1.0);
qa(0:J.n-1) -= h*yqP0;
if (mpirank == 0) qa(Ja.n-1) += h;
if (contorder > 1) {
real[int,int] yqnP(J.n, contorder-1);
yqP = -h*yqP0;
real[int] NLq(ub[].n), temp(ub[].n), hn(contorder-1);
real cbrteps = cbrt(eps), sqrteps = sign(eps)*sqrt(abs(eps));
for (int i = 0; i < contorder-1; i++){
yqnP(:, i) = yqP;
hn(i) = h;
NLq = 0.0;
for (int j = 0; j < i+1; j++){
ChangeNumbering(J, um[], yqnP(:, j), inverse = true, exchange = true);
ChangeNumbering(J, um2[], yqnP(:, i-j), inverse = true, exchange = true);
um3[] = vH(0, XMh, tgv = -10);
NLq += um3[];
updateparam(param, paramval + 2.0*sqrteps);
um2[] = vR(0, XMh, tgv = TGV);
updateparam(param, paramval + sqrteps);
um3[] = vR(0, XMh, tgv = TGV);
um2[] += (R - 2.0*um3[]);
NLq += (hn(j)*hn(i-j)/eps)*um2[];
updateparam(param, paramval + eps);
um3[] = vJ(0, XMh, tgv = -10);
updateparam(param, paramval);
um2[] = vJ(0, XMh, tgv = -10);
um3[] -= um2[];
NLq += (2.0*hn(i-j)/eps)*um3[];
for (int k = 0; k < i-j; k++){
ChangeNumbering(J, um2[], yqnP(:, k), inverse = true, exchange = true);
IFMACRO(cubic)
ChangeNumbering(J, um3[], yqnP(:, i-j-k-1), inverse = true, exchange = true);
temp = vT(0, XMh, tgv = -10);
NLq += temp/3.0;
ENDIFMACRO
updateparam(param, paramval + eps);
um3[] = vH(0, XMh, tgv = -10);
updateparam(param, paramval);
temp = vH(0, XMh, tgv = -10);
um3[] -= temp;
NLq += (hn(i-j-k-1)/eps)*um3[];
um3[] = vJ(0, XMh, tgv = -10);
updateparam(param, paramval + 2.0*sqrteps);
um2[] = vJ(0, XMh, tgv = -10);
um2[] += um3[];
updateparam(param, paramval + sqrteps);
um3[] = vJ(0, XMh, tgv = -10);
um2[] -= 2.0*um3[];
NLq += (hn(k)*hn(i-j-k-1)/eps)*um2[];
updateparam(param, paramval + 3.0*cbrteps);
um2[] = vR(0, XMh, tgv = TGV);
updateparam(param, paramval + 2.0*cbrteps);
um3[] = vR(0, XMh, tgv = TGV);
um2[] -= 3.0*um3[];
updateparam(param, paramval + cbrteps);
um3[] = vR(0, XMh, tgv = TGV);
um2[] += (3.0*um3[] - R);
NLq += (hn(j)*hn(k)*hn(i-j-k-1)/(6.0*eps))*um2[];
updateparam(param, paramval);
}
}
NLq *= -0.5;
ChangeNumbering(J, NLq, yqP);
KSPSolve(J, yqP, yqP);
hl = (hn(0)/h0)^2.0*(yqP'*yqP0);
mpiAllReduce(hl, h, mpiCommWorld, mpiSUM);
yqP -= h*yqP0;
qa(0:J.n-1) += yqP;
if (mpirank == 0) qa(Ja.n-1) += h;
}
}
}
// CORRECTOR LOOP
adaptflag = false;
SNESSolve(Ja, funcJa, funcRa, qa, convergence = funcConvergence, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -snes_converged_reason -options_left no -snes_max_it " + snesmaxit); // solve nonlinear problem with SNES
if (ret > 0) {
++count;
if (maxcount > 0) stopflag = (count >= maxcount);
else if ((paramval - paramtarget)*paramdiff <= 0) stopflag = true;
h0 /= f;
if (cosalpha < 0 && contorder > 0) {
h0 *= -1.0;
if(mpirank == 0) cout << "\tFold bifurcation detected. Orientation reversed." << endl;
forcesave = true;
}
if (adapt && (count % savecount == 0)){
meshout = meshroot + "_" + count + "." + meshext;
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true);
if(mpirank == 0) paramval = qa(Ja.n-1);
broadcast(processor(0), paramval);
updateparam(param, paramval);
ChangeNumbering(J, um[], yqP, inverse = true);
XMhg defu(uG), defu(yG), defu(tempu);
tempu[](restu) = ub[]; // populate local portion of global soln
mpiAllReduce(tempu[], uG[], mpiCommWorld, mpiSUM);
tempu[](restu) = um[]; // populate local portion of global tangent vector
mpiAllReduce(tempu[], yG[], mpiCommWorld, mpiSUM);
if(mpirank == 0) {
IFMACRO(dimension,2)
if(adaptto == "b") Thg = adaptmesh(Thg, adaptu(uG), adaptmeshoptions);
else if(adaptto == "by") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), 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((bool(anisomax > 1) ? 6 : 1)*Thg.nv);
if(adaptto == "b") 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));
else if(adaptto == "by") met = mshmet(Thg, adaptu(uG), adaptu(yG), 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);
defu(uG) = defu(uG);
defu(yG) = defu(yG);
Th = Thg;
Mat Adapt;
createMatu(Th, Adapt, Pk);
J = Adapt;
defu(ub) = initu(0.0);
defu(um) = initu(0.0);
defu(um2) = initu(0.0);
defu(um3) = initu(0.0);
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);
um[] = yG[](restu);
Mat Adapt1(J.n, mpirank == 0 ? 1 : 0), Adapt2(J.n, mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
JlPM = Adapt1;
yqPM = Adapt2;
Ja = [[J, JlPM], [yqPM', -1.0]]; // make dummy Jacobian
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO prefix = "fieldsplit_0_", parent = Ja);
ChangeNumbering(J, ub[], qa);
qa.resize(Ja.n);
if(mpirank == 0) qa(Ja.n-1) = paramval;
R.resize(ub[].n);
ChangeNumbering(J, um[], yqP);
yqP0.resize(J.n);
yqP0 = yqP;
qa0.resize(Ja.n);
if (contorder == 0) {
matrix<PetscScalar> tempPms = [[yqP]];
ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
ChangeOperator(yqPM, tempPms, parent = Ja); // send to Mat
}
adaptflag = true;
SNESSolve(Ja, funcJa, funcRa, qa, convergence = funcConvergence, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -snes_converged_reason -options_left no"); // solve nonlinear problem with SNES
assert(ret > 0);
if(mpirank==0) { // Save adapted mesh
cout << " Saving adapted mesh '" + meshout + "' in '" + workdir + "'." << endl;
savemesh(Thg, workdir + meshout);
}
}
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true);
if(mpirank == 0) paramval = qa(Ja.n-1);
broadcast(processor(0), paramval);
updateparam(param, paramval);
savebase(fileout + "_" + count + (forcesave ? "specialpt" : ""), (savecount > 0 ? fileout : ""), meshout, ((count % savecount == 0) || forcesave || stopflag), true);
if (stopflag) break;
forcesave = false;
it = 0;
if (stricttangent && contorder > 0) funcJa(qa);
else ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
yqP0 = yqP;
qa0 = qa;
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
}
else h0 /= fmax;
}