Author: Chris Douglas (@cmdoug) christopher.douglas@duke.edu
ff-mpirun -np 4 hopfcontinue.md -param <PARAM1> -param2 <PARAM2> -fi <FILEIN> -fo <FILEOUT>ff-mpirun -np 4 hopfcontinue.md -param <PARAM1> -param2 <PARAM2> -fi <FILEIN> -fo <FILEOUT> -mo <MESHOUT>NOTE: This file should not be changed unless you know what you're doing.
SEE ALSO: modecompute.md, hopfcompute.md, fohocompute.md, ./botacompute.md, hohocompute.md, porbcontinue.md
load "iovtk"
load "PETSc-complex"
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 fileout = getARGV("-fo", filein);
int contorder = getARGV("-contorder", 1);
bool normalform = getARGV("-nf", 1);
bool wnlsave = getARGV("-wnl", 0);
int select = getARGV("-select", 1);
bool zerofreq = getARGV("-zero", 0);
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 param2 = getARGV("-param2", "");
string adaptto = getARGV("-adaptto", "b");
real fmax = getARGV("-fmax", 2.0);
real kappamax = getARGV("-kmax", 1.0);
real deltamax = getARGV("-dmax", 4.0);
real anglemax = getARGV("-amax", 30.)*pi/180.0;
real monotone = getARGV("-mono", 0.0);
real eps = getARGV("-eps", 1.0e-7);
real eps2 = getARGV("-eps2", 1.0e-7);
bool stricttangent = bool(getARGV("-stricttangent", 1));
int snesmaxit = getARGV("-snes_max_it", 10);
string sneslinesearchtype = getARGV("-snes_linesearch_type", "none");
real TGV = getARGV("-tgv", -1);
real[int] sym1(sym.n);
real omega;
complex[string] alpha;
complex beta;
real paramtarget = getARGV("-paramtarget",1.0);
real param2target = getARGV("-param2target",1.0);
bool stopflag = false;
bool forcesave = false;
// 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
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
}
assert(fileext == "hopf" || fileext == "foho" || fileext == "hoho" || fileext == "bota" || fileext == "baut");
Th = readmeshN(workdir + meshin);
Thg = Th;
DmeshCreate(Th);
restu = restrict(XMh, XMhg, n2o);
XMh<complex> defu(ub), defu(um), defu(uma), defu(um2), defu(um3);
if (count == 0){
if( fileext == "hopf"){
ub[].re = loadhopf(fileroot, meshin, um[], uma[], sym1, omega, alpha, beta);
}
else if(fileext == "bota") {
real[string] alpha1, alpha2;
real beta1, beta2, beta3, beta4;
ub[].re = loadbota(fileroot, meshin, um[].re, uma[].re, alpha1, alpha2, beta1, beta2, beta3, beta4);
}
else if (fileext == "foho") {
real[string] alpha2;
real beta22, beta23, gamma22, gamma23;
complex gamma12, gamma13;
real[int] q2m, q2ma;
ub[].re = loadfoho(fileroot, meshin, um[], uma[], q2m, q2ma, sym1, omega, alpha, alpha2, beta, beta22, beta23, gamma12, gamma13, gamma22, gamma23);
}
else if(fileext == "hoho") {
real omegaN;
complex[string] alphaN;
complex betaN, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23;
complex[int] qNm, qNma;
if(select == 1){
ub[].re = loadhoho(fileroot, meshin, um[], uma[], qNm, qNma, sym1, sym, omega, omegaN, alpha, alphaN, beta, betaN, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23);
}
else if(select == 2){
ub[].re = loadhoho(fileroot, meshin, qNm, qNma, um[], uma[], sym, sym1, omegaN, omega, alphaN, alpha, betaN, beta, gamma11, gamma12, gamma13, gamma21, gamma22, gamma23);
}
}
savehopf(filein, (savecount > 0 ? fileout : ""), meshin, sym1, omega, alpha, beta, false, false);
}
else {
ub[].re = loadhopf(fileroot + "_" + count, meshin, um[], uma[], sym1, omega, alpha, beta);
}
real[int] paramvals(3-zerofreq);
paramvals(0) = getparam(param);
paramvals(1) = zerofreq ? 0.0 : omega;
paramvals(2-zerofreq) = getparam(param2);
real paramdiff1 = paramvals(0) - paramtarget;
real paramdiff2 = paramvals(2-zerofreq) - param2target;
// Create distributed Mat
Mat<complex> J;
createMatu(Th, J, Pk);
complex[int] ik(sym.n), ik2(sym.n), ik3(sym.n);
complex iomega, 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<complex> JlPM(J.n, mpirank == 0 ? (2-zerofreq) : 0), gqPM(J.n, mpirank == 0 ? (2-zerofreq) : 0), glPM(mpirank == 0 ? (2-zerofreq) : 0, mpirank == 0 ? (2-zerofreq) : 0); // Initialize Mat objects for bordered matrix
Mat<complex> JlPMa(J.n + (mpirank == 0 ? (2-zerofreq) : 0), mpirank == 0 ? 1 : 0), yqPMa(J.n + (mpirank == 0 ? (2-zerofreq) : 0), mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
Mat<complex> Ja = [[J, JlPM], [gqPM', glPM]], Jaa = [[Ja, JlPMa], [yqPMa', -1.0]]; // make dummy Jacobian
complex[int] R(um[].n), qm(J.n), qma(J.n), pP(J.n), qP(J.n), yqP(Ja.n), yqP0(Ja.n);
int ret, it = 0;
real f, kappa, cosalpha, res, delta, maxdelta, omega0;
complex alpha0, beta0;
// 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) paramvals = qa(J.n:Jaa.n-1).re;
broadcast(processor(0), paramvals);
updateparam(param, paramvals(0));
omega = zerofreq ? 0.0 : paramvals(1);
updateparam(param2, paramvals(2-zerofreq));
sym = 0;
R = vR(0, XMh, tgv = TGV);
PetscScalar[int] Ra;
ChangeNumbering(J, R, Ra); // FreeFEM to PETSc
sym = sym1;
iomega = 1i*omega;
ik.im = sym1;
J = vJ(XMh, XMh, tgv = -2);
KSPSolve(J, pP, qm);
KSPSolveHermitianTranspose(J, qP, qma);
PetscScalar ginv, ginvl = (qP'*qm);
mpiAllReduce(ginvl, ginv, mpiCommWorld, mpiSUM);
qm /= ginv; // rescale direct mode
qma /= conj(ginv); // rescale adjoint mode
Ra.resize(Jaa.n); // Append 0 to residual vector on proc 0
if(mpirank == 0) {
Ra(J.n) = real(1.0/ginv);
if (!zerofreq) Ra(Jaa.n-2) = imag(1.0/ginv);
Ra(Jaa.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) paramvals = qa(J.n:Jaa.n-1).re;
broadcast(processor(0), paramvals);
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
ChangeNumbering(J, uma[], qma, inverse = true);
updateparam(param, paramvals(0) + eps);
omega = zerofreq ? 0.0 : paramvals(1);
updateparam(param2, paramvals(2-zerofreq));
sym = 0;
um2[] = vR(0, XMh, tgv = TGV);
um2[] -= R;
um2[] /= eps;
ChangeNumbering(J, um2[], qm); // FreeFEM to PETSc
matrix<PetscScalar> tempPms;
if(zerofreq) tempPms = [[qm]];
else tempPms = [[qm, 0]]; // dense array to sparse matrix
ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
sym = sym1;
ik.im = sym1;
iomega = 1i*omega;
um2[] = vJ(0, XMh, tgv = -10);
updateparam(param, paramvals(0));
um3[] = vJ(0, XMh, tgv = -10);
um2[] -= um3[];
PetscScalar gl = J(uma[], um2[])/eps;
if(zerofreq) tempPms = [[real(gl)]];
else {
um2[] = vM(0, XMh, tgv = -10);
PetscScalar gw = J(uma[], um2[]);
tempPms = [[real(gl), -imag(gw)], [imag(gl), real(gw)]];
}
ChangeOperator(glPM, tempPms, parent = Ja); // send to Mat
sym = 0;
J = vH(XMh, XMh, tgv = 0); // form the matrix (dL/dq*w)
MatMultHermitianTranspose(J, qma, qm); // gqr,i
if(!zerofreq) qma.re = -qm.im;
qm.im = 0.0;
qma.im = 0.0;
if(zerofreq) tempPms = [[qm]];
else tempPms = [[qm, qma]]; // dense array to sparse matrix
ChangeOperator(gqPM, tempPms, parent = Ja); // send to Mat
ik = 0.0;
iomega = 0.0;
J = vJ(XMh, XMh, tgv = TGV);
if (contorder > 0) {
updateparam(param2, paramvals(2-zerofreq) + eps2);
um2[] = vR(0, XMh, tgv = TGV);
um2[] -= R;
ChangeNumbering(J, um2[], qm); // FreeFEM to PETSc
yqP(0:J.n-1) = qm/eps2;
sym = sym1;
ik.im = sym1;
iomega = 1i*omega;
um2[] = vJ(0, XMh, tgv = -10);
updateparam(param2, paramvals(2-zerofreq));
um2[] -= um3[];
gl = J(uma[], um2[])/eps2;
if(mpirank == 0) {
yqP(J.n) = real(gl);
if(!zerofreq) yqP(Ja.n-1) = imag(gl);
}
tempPms = [[yqP]]; // dense array to sparse matrix
ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
KSPSolve(Ja, yqP, yqP);
tempPms = [[yqP]]; // dense array to sparse matrix
ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
}
return 0;
}
ConvergenceCheck(param + " = " + paramvals(0) + ",\t" + param2 + " = " + paramvals(2-zerofreq) + ",\tomega = " + (zerofreq ? 0.0 : paramvals(1)));
// set up Mat parameters
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(Jaa, 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_ -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition full -prefix_pop"
+ " -prefix_push fieldsplit_0_fieldsplit_1_ -ksp_type preonly -pc_type redundant -redundant_pc_type lu -prefix_pop"
+ " -prefix_push fieldsplit_0_fieldsplit_0_ " + KSPparams + " -prefix_pop", setup = 1);
set(Ja, prefix = "fieldsplit_0_", setup = 1, parent = Jaa);
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO prefix = "fieldsplit_0_fieldsplit_0_", parent = Ja);
// PREDICTOR
complex[int] qa(Jaa.n), qa0;
ChangeNumbering(J, ub[], qa0);
ChangeNumbering(J, um[], qm);
ChangeNumbering(J, uma[], qma);
qa0.resize(Jaa.n);
if(mpirank == 0) qa0(J.n:Jaa.n-1).re = paramvals;
sym = sym1;
ik.im = sym1;
J = vM(XMh, XMh, tgv = 0);
MatMult(J, qm, qP);
MatMultHermitianTranspose(J, qma, pP);
if (contorder > 0) {
sym = 0;
R = vR(0, XMh, tgv = TGV);
funcJa(qa0);
}
else {
matrix<PetscScalar> tempPms = [[yqP]];
ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
}
yqP0 = yqP;
omega0 = omega;
alpha0 = alpha[param];
beta0 = beta;
while (!stopflag){
qa = qa0;
if (contorder == 0 && mpirank == 0) qa(Jaa.n-1) += h0;
else if (contorder > 0) {
real h, hl = real(yqP0'*yqP0);
mpiAllReduce(hl, h, mpiCommWorld, mpiSUM);
h = h0/sqrt(h + 1.0);
qa(0:Ja.n-1) -= (h*yqP0);
if (mpirank == 0) qa(Jaa.n-1) += h;
}
// CORRECTOR LOOP
adaptflag = false;
SNESSolve(Jaa, funcJa, funcRa, qa, convergence = funcConvergence, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -options_left no -snes_converged_reason -snes_max_it " + snesmaxit); // solve nonlinear problem with SNES
if (ret > 0) {
++count;
if (maxcount > 0) stopflag = (count >= maxcount);
else if ((paramvals(0) - paramtarget)*paramdiff1 <= 0 || (paramvals(2-zerofreq) - param2target)*paramdiff2 <= 0) stopflag = true;
h0 /= f;
if (cosalpha < 0 && contorder > 0) {
h0 *= -1.0;
if(mpirank == 0) cout << "\tOrientation 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) paramvals = qa(J.n:Jaa.n-1).re;
broadcast(processor(0), paramvals);
updateparam(param, paramvals(0));
omega = zerofreq ? 0.0 : paramvals(1);
updateparam(param2, paramvals(2-zerofreq));
ChangeNumbering(J, um[], qm, inverse = true);
ChangeNumbering(J, uma[], qma, inverse = true);
ChangeNumbering(J, um2[], yqP(0:J.n-1), inverse = true);
complex[int] yparamvals(2-zerofreq);
if(mpirank == 0) yparamvals = yqP(J.n:Ja.n-1);
XMhg defu(uG), defu(umrG), defu(umiG), defu(umarG), defu(umaiG), defu(yG), defu(tempu), defu(uoG);
tempu[](restu) = ub[].re; // populate local portion of global soln
mpiAllReduce(tempu[], uG[], mpiCommWorld, mpiSUM);
tempu[](restu) = um[].re; // populate local portion of global soln
mpiAllReduce(tempu[], umrG[], mpiCommWorld, mpiSUM);
tempu[](restu) = um[].im; // populate local portion of global soln
mpiAllReduce(tempu[], umiG[], mpiCommWorld, mpiSUM);
tempu[](restu) = uma[].re; // populate local portion of global soln
mpiAllReduce(tempu[], umarG[], mpiCommWorld, mpiSUM);
tempu[](restu) = uma[].im; // populate local portion of global soln
mpiAllReduce(tempu[], umaiG[], mpiCommWorld, mpiSUM);
tempu[](restu) = um2[].re; // populate local portion of global soln
mpiAllReduce(tempu[], yG[], mpiCommWorld, mpiSUM);
if(mpirank == 0) { // Perform mesh adaptation (serially) on processor 0
if(adaptto == "bo") {
defu(uoG) = initu(defu(umarG)'*defu(umarG));
defu(tempu) = initu(defu(umaiG)'*defu(umaiG));
tempu[] += uoG[];
tempu[] = sqrt(tempu[]);
uoG[] = (umrG[].*umrG[]);
uoG[] += (umiG[].*umiG[]);
uoG[] = sqrt(uoG[]);
uoG[] .*= tempu[];
}
IFMACRO(dimension,2)
if(adaptto == "b") Thg = adaptmesh(Thg, adaptu(uG), adaptmeshoptions);
else if(adaptto == "by") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptmeshoptions);
else if(adaptto == "bd") Thg = adaptmesh(Thg, adaptu(uG), adaptu(umrG), adaptu(umiG), adaptmeshoptions);
else if(adaptto == "ba") Thg = adaptmesh(Thg, adaptu(uG), adaptu(umarG), adaptu(umaiG), adaptmeshoptions);
else if(adaptto == "byd") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptu(umrG), adaptu(umiG), adaptmeshoptions);
else if(adaptto == "bya") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptu(umarG), adaptu(umaiG), adaptmeshoptions);
else if(adaptto == "bda") Thg = adaptmesh(Thg, adaptu(uG), adaptu(umrG), adaptu(umiG), adaptu(umarG), adaptu(umaiG), adaptmeshoptions);
else if(adaptto == "byda") Thg = adaptmesh(Thg, adaptu(uG), adaptu(yG), adaptu(umrG), adaptu(umiG), adaptu(umarG), adaptu(umaiG), adaptmeshoptions);
else if(adaptto == "bo") Thg = adaptmesh(Thg, adaptu(uG), adaptu(uoG), 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));
else if(adaptto == "bd") met = mshmet(Thg, adaptu(uG), adaptu(umrG), adaptu(umiG), 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 == "ba") met = mshmet(Thg, adaptu(uG), adaptu(umarG), adaptu(umaiG), 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 == "byd") met = mshmet(Thg, adaptu(uG), adaptu(yG), adaptu(umrG), adaptu(umiG), 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 == "bya") met = mshmet(Thg, adaptu(uG), adaptu(yG), adaptu(umarG), adaptu(umaiG), 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 == "bda") met = mshmet(Thg, adaptu(uG), adaptu(umrG), adaptu(umiG), adaptu(umarG), adaptu(umaiG), 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 == "byda") met = mshmet(Thg, adaptu(uG), adaptu(yG), adaptu(umrG), adaptu(umiG), adaptu(umarG), adaptu(umaiG), 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 == "bo") met = mshmet(Thg, adaptu(uG), adaptu(uoG), 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(yG) = defu(yG); //interpolate global solution from old mesh to new mesh
defu(umrG) = defu(umrG);
defu(umiG) = defu(umiG);
defu(umarG) = defu(umarG);
defu(umaiG) = defu(umaiG);
Th = Thg;
Mat<complex> Adapt;
createMatu(Th, Adapt, Pk);
J = Adapt;
defu(ub) = initu(0.0);
defu(um) = initu(0.0);
defu(uma) = 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[].re = uG[](restu);
um2[].re = yG[](restu);
um[].re = umrG[](restu);
um[].im = umiG[](restu);
uma[].re = umarG[](restu);
uma[].im = umaiG[](restu);
Mat<complex> Adapt1(J.n, mpirank == 0 ? (2-zerofreq) : 0), Adapt2(J.n, mpirank == 0 ? (2-zerofreq) : 0); // Initialize Mat objects for bordered matrix
Mat<complex> Adapt3(J.n + (mpirank == 0 ? (2-zerofreq) : 0), mpirank == 0 ? 1 : 0), Adapt4(J.n + (mpirank == 0 ? (2-zerofreq) : 0), mpirank == 0 ? 1 : 0); // Initialize Mat objects for bordered matrix
JlPM = Adapt1;
gqPM = Adapt2;
JlPMa = Adapt3;
yqPMa = Adapt4;
Ja = [[J, JlPM], [gqPM', glPM]]; // make dummy Jacobian
Jaa = [[Ja, JlPMa], [yqPMa', -1.0]]; // make dummy Jacobian
IFMACRO(Jprecon) Jprecon(0); ENDIFMACRO
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO prefix = "fieldsplit_0_fieldsplit_0_", parent = Ja);
set(Ja, prefix = "fieldsplit_0_", parent = Jaa);
ChangeNumbering(J, ub[], qa);
qa.resize(Jaa.n);
if(mpirank == 0) qa(J.n:Jaa.n-1).re = paramvals;
ChangeNumbering(J, um[], qm);
ChangeNumbering(J, uma[], qma);
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um3[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um3[], qP);
ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
um3[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um3[], pP);
R.resize(ub[].n);
ChangeNumbering(J, um2[], yqP);
yqP.resize(Ja.n);
if (mpirank==0) yqP(J.n:Ja.n-1) = yparamvals;
yqP0.resize(Ja.n);
yqP0 = yqP;
qa0.resize(Jaa.n);
if (contorder == 0) {
matrix<PetscScalar> tempPms = [[yqP]];
ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
}
adaptflag = true;
SNESSolve(Jaa, funcJa, funcRa, qa, convergence = funcConvergence, reason = ret,
sparams = "-snes_linesearch_type " + sneslinesearchtype + " -options_left no -snes_converged_reason"); // 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, exchange = true);
if(mpirank == 0) paramvals = qa(J.n:Jaa.n-1).re;
broadcast(processor(0), paramvals);
updateparam(param, paramvals(0));
omega = zerofreq ? 0.0 : paramvals(1);
updateparam(param2, paramvals(2-zerofreq));
J = vM(XMh, XMh, tgv = 0);
ChangeNumbering(J, um[], qm);
MatMult(J, qm, qP);
complex phaseref, phaserefl = qP.sum;
mpiAllReduce(phaserefl, phaseref, mpiCommWorld, mpiSUM);
qm /= phaseref;
qP /= phaseref;
real Mnorm, local = real(qm'*qP);
mpiAllReduce(local, Mnorm, mpiCommWorld, mpiSUM);
local = sqrt(Mnorm);
qP /= local;
qm /= local;
phaserefl = (qP'*qma);
mpiAllReduce(phaserefl, phaseref, mpiCommWorld, mpiSUM);
qma /= phaseref;
ChangeNumbering(J, uma[], qma, inverse = true);
if (normalform){
complex[int] temp(um[].n);
complex[int,int] qDa(paramnames.n, J.n);
// 2nd-order
// A: base modification due to parameter changes
ik = 0.0;
ik2 = 0.0;
iomega = 0.0;
iomega2 = 0.0;
sym = 0;
J = vJ(XMh, XMh, tgv = TGV);
if(paramnames[0] != ""){
for (int k = 0; k < paramnames.n; ++k){
real paramval = getparam(paramnames[k]);
updateparam(paramnames[k], paramval + eps);
um[] = vR(0, XMh, tgv = TGV);
updateparam(paramnames[k], paramval);
um[] -= R;
um[] /= -eps;
ChangeNumbering(J, um[], qP); // FreeFEM to PETSc
KSPSolve(J, qP, qP);
qDa(k, :) = qP;
}
}
// B: base modification due to quadratic nonlinear interaction
ik.im = sym1;
ik2.im = -sym1;
iomega = 1i*omega;
iomega2 = -iomega;
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
um2[] = conj(um[]);
um3[] = vH(0, XMh, tgv = -10);
um3[].re *= -1.0; // -2.0/2.0
um3[].im = 0.0;
ChangeNumbering(J, um3[], qP); // FreeFEM to PETSc
KSPSolve(J, qP, qP);
// C: harmonic generation due to quadratic nonlinear interaction
ik2.im = sym1;
iomega2 = iomega;
um2[] = -0.5*um[];
sym = 2*sym1;
um3[] = vH(0, XMh, tgv = -10);
ChangeNumbering(J, um3[], pP); // FreeFEM to PETSc
ik.im = 2*sym1;
iomega = 2i*omega;
J = vJ(XMh, XMh, tgv = TGV);
KSPSolve(J, pP, pP);
// 3rd-order
// A: fundamental modification due to parameter change and quadratic interaction of fundamental with 2nd order modification A.
sym = sym1;
ik.im = sym1;
ik2 = 0.0;
iomega = 1i*omega;
iomega2 = 0.0;
if(paramnames[0] != ""){
temp = vJ(0, XMh, tgv = -10);
for (int k = 0; k < paramnames.n; ++k){
real paramval = getparam(paramnames[k]);
ChangeNumbering(J, um2[], qDa(k, :), inverse = true, exchange = true); // FreeFEM to PETSc
um3[] = vH(0, XMh, tgv = -10); // 2.0/2.0
updateparam(paramnames[k], paramval + eps);
um2[] = vJ(0, XMh, tgv = -10);
updateparam(paramnames[k], paramval);
um2[] -= temp;
um3[] += um2[]/eps;
alpha[paramnames[k]] = J(uma[], um3[]);
}
}
IFMACRO(cubic)
// B: fundamental modification due to cubic self-interaction of fundamental
ik2.im = sym1;
ik3.im = -sym1;
iomega2 = iomega;
iomega3 = -iomega;
um2[] = 0.5*um[];
um3[] = conj(um[]);
temp = vT(0, XMh, tgv = -10);
ENDIFMACRO
// C: fundamental modification due to quadratic interaction of fundamental with 2nd order modification B
ik2 = 0.0;
iomega2 = 0.0;
ChangeNumbering(J, um2[], qP, inverse = true, exchange = true); // FreeFEM to PETSc
IFMACRO(cubic)
um3[] = vH(0, XMh, tgv = -10);
temp += um3[];
ENDIFMACRO
IFMACRO(!cubic)
temp = vH(0, XMh, tgv = -10);
ENDIFMACRO
// D: fundamental modification due to quadratic interaction of fundamental with 2nd order modification C
ik.im = -sym1;
ik2.im = 2*sym1;
iomega = -iomega;
iomega2 = 2i*omega;
um[] = conj(um[]);
ChangeNumbering(J, um2[], pP, inverse = true, exchange = true); // FreeFEM to PETSc
um3[] = vH(0, XMh, tgv = -10);
temp += um3[];
beta = J(uma[], temp);
ik2 = 0.0;
iomega2 = 0.0;
if(wnlsave){
complex[int] val(1);
XMh<complex>[int] defu(vec)(1);
sym = 0;
val = 0.0;
if(paramnames[0] != ""){
for (int k = 0; k < paramnames.n; ++k){
ChangeNumbering(J, vec[0][], qDa(k, :), inverse = true); // FreeFEM to PETSc
savemode(fileout + "_" + count + "_wnl_param" + k, "", fileout + "_" + count + ".hopf", meshout, vec, val, sym, true);
}
}
ChangeNumbering(J, vec[0][], qP, inverse = true); // FreeFEM to PETSc
savemode(fileout + "_" + count + "_wnl_AAs", "", fileout + ".hopf", meshout, vec, val, sym, true);
ChangeNumbering(J, vec[0][], pP, inverse = true); // FreeFEM to PETSc
val = 2i*omega;
sym = 2*sym1;
savemode(fileout + "_" + count + "_wnl_AA", "", fileout + ".hopf", meshout, vec, val, sym, true);
}
}
else {
if(paramnames[0] != ""){
for (int k = 0; k < paramnames.n; ++k){
alpha[paramnames[k]] = 0.0;
}
}
beta = 0.0;
}
if (real(beta)*real(beta0) < 0) {
forcesave = true;
if(mpirank == 0) {
if(real(alpha0)*real(alpha[param]) < 0) {
if ( omega*omega0 < 0 || abs(omega) < 1.0e-6)
cout << "\tBogdanov-Takens bifurcation (or zero-frequency point) detected." << endl;
else cout << "\tFold-Hopf bifurcation detected." << endl;
}
else cout << "\tBautin bifurcation detected." << endl;
}
}
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true);
ChangeNumbering(J, um[], qm, inverse = true);
savehopf(fileout + "_" + count + (forcesave ? "specialpt" : ""), fileout, meshout, sym1, omega, alpha, beta, ((count % savecount == 0) || forcesave || stopflag), true);
if (stopflag) break;
forcesave = false;
ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
sym = sym1;
ik.im = sym1;
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um2[], qP);
ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
um2[] = vM(0, XMh, tgv = 0);
ChangeNumbering(J, um2[], pP);
it = 0;
if (stricttangent && contorder > 0) funcJa(qa);
yqP0 = yqP;
qa0 = qa;
omega0 = omega;
alpha0 = alpha[param];
beta0 = beta;
}
else h0 /= fmax;
}