hopfcontinue.edp

//
// hopfcontinue.edp
// Chris Douglas
// christopher.douglas@duke.edu
//
// EXAMPLE USAGE:
//  Initialize with fold from file, solve on same mesh
// ff-mpirun -np 4 hopfcontinue.edp -param <PARAM1> -param2 <PARAM2> -fi <FILEIN> -fo <FILEOUT>
//
//  Initialize with fold from file, adapt mesh/solution
// ff-mpirun -np 4 hopfcontinue.edp -param <PARAM1> -param2 <PARAM2> -fi <FILEIN> -fo <FILEOUT> -mo <MESHOUT>
//
// NOTE: This file should not be changed unless you know what you're doing.
//
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);
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;
int monotone = getARGV("-mono", 0);
real eps = getARGV("-eps", 1.0e-7);
real eps2 = getARGV("-eps2", 1.0e-7);
int snesmaxit = getARGV("-snes_max_it", 10);
string sneslinesearchtype = getARGV("-snes_linesearch_type","basic");
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");
Th = readmeshN(workdir + meshin);
Thg = Th;
buildDmesh(Th);
restu = restrict(XMh, XMhg, n2o);
XMh<complex> defu(ub), defu(um), defu(uma), defu(yb), defu(um2), defu(um3);
if (count == 0){
  if( fileext == "hopf"){
    ub[].re = loadhopf(fileroot, meshin, um[], uma[], sym1, omega, alpha, beta);
  }
  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, gamma1, gamma2, gamma12, gamma13, 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, gamma1, gamma2, gamma12, gamma13, gamma22, gamma23);
    }
    else if(select == 2){
      ub[].re = loadhoho(fileroot, meshin, qNm, qNma, um[], uma[], sym, sym1, omegaN, omega, alphaN, alpha, betaN, beta, gamma1, gamma2, gamma12, gamma13, 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 lambda1 = getlambda(param);
real lambda2 = getlambda(param2);
real paramdiff1 = lambda1 - paramtarget;
real paramdiff2 = lambda2 - 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 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> H(J), 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), qap(Jaa.n);
int it, internalit, adaptflag;
real f, kappa, cosalpha, res, resp, delta, deltap, maxdelta, omega0;
complex alpha0, beta0;
// FUNCTIONS
  func complex[int] funcRa(complex[int]& qa) {
      ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true); // PETSc to FreeFEM
      if(mpirank == 0) {
        lambda1 = real(qa(qa.n-(3-zerofreq))); // Extract parameter value from state vector on proc 0
        omega = zerofreq ? 0.0 : real(qa(qa.n-2)); // Extract frequency value from state vector on proc 0
        lambda2 = real(qa(qa.n-1)); // Extract parameter value from state vector on proc 0
      }
      broadcast(processor(0), lambda1);
      broadcast(processor(0), omega);
      broadcast(processor(0), lambda2);
      updatelambda(param, lambda1);
      updatelambda(param2, lambda2);
      sym = 0;
      R = vR(0, XMh, tgv = -1);
      complex[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);
      complex ginv, ginvl = (qP'*qm);
      mpiAllReduce(ginvl, ginv, mpiCommWorld, mpiSUM);
      qm /= ginv; // rescale direct mode
      qma /= conj(ginv); // rescale adjoint mode
      if(mpirank == 0) {
        Ra.resize(J.n+(3-zerofreq)); // Append 0 to residual vector on proc 0
        Ra(Ra.n-(3-zerofreq)) = real(1.0/ginv);
        if (!zerofreq) Ra(Ra.n-2) = imag(1.0/ginv);
        Ra(Ra.n-1) = 0.0;
      }
      StepAdaptMonitors(Ra, qa, qap, yqP, yqP0);
      if(mpirank == 0) cout << "  " + text1 + ":\t||R|| = " << res << (it == 0 ? (",\th0 = " + h0) : (",\t||dx|| = " + delta + ",\tangle = " + (sign(cosalpha)*acos(abs(cosalpha))*180./pi))) << ",\t" + param + " = " << lambda1 << ",\t" + param2 + " = " << lambda2 << ",\tomega = " << omega << "." << endl;
      return Ra;
  }

  func int funcJa(complex[int]& qa) {
      ++it;
      internalit = 0;
      qap = qa;
      resp = res;
      deltap = delta;
      ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true); // PETSc to FreeFEM
      if(mpirank == 0) {
        lambda1 = real(qa(qa.n-(3-zerofreq))); // Extract parameter value from state vector on proc 0
        omega = zerofreq ? 0.0 : real(qa(qa.n-2)); // Extract frequency value from state vector on proc 0
        lambda2 = real(qa(qa.n-1)); // Extract parameter value from state vector on proc 0
      }
      broadcast(processor(0), lambda1);
      broadcast(processor(0), omega);
      broadcast(processor(0), lambda2);
      ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
      ChangeNumbering(J, uma[], qma, inverse = true);
      ik.im = sym1;
      iomega = 1i*omega;
      updatelambda(param, lambda1 + eps);
      updatelambda(param2, lambda2);
      sym = 0;
      complex[int] Jl1 = vR(0, XMh, tgv = -1);
      sym = sym1;
      complex[int] Hl1 = vJ(0, XMh, tgv = -10);
      updatelambda(param, lambda1);
      updatelambda(param2, lambda2 + eps2);
      sym = 0;
      complex[int] Jl2 = vR(0, XMh, tgv = -1);
      sym = sym1;
      complex[int] Hl2 = vJ(0, XMh, tgv = -10);
      updatelambda(param2, lambda2);
      Jl1 -= R;
      Jl2 -= R;
      Jl1 /= eps;
      Jl2 /= eps2;
      complex[int] temp1(J.n), temp2(J.n), bm = vJ(0, XMh, tgv = -10);
      Hl1 -= bm;
      Hl2 -= bm;
      Hl1 /= eps;
      Hl2 /= eps2;
      ChangeNumbering(J, Jl1, temp1); // FreeFEM to PETSc
      temp1.im = 0.0;
      matrix<complex> tempPms;
      if(zerofreq) tempPms = [[temp1]];
      else tempPms = [[temp1, 0]]; // dense array to sparse matrix
      ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
      sym = 0;
      H = vH(XMh, XMh, tgv = 0); // form the matrix (dL/dq*w)
      sym = sym1;
      MatMultHermitianTranspose(H, qma, temp1); // gqr,i
      if(!zerofreq) temp2.re = -temp1.im;
      temp1.im = 0.0;
      if(zerofreq) tempPms = [[temp1]];
      else tempPms = [[temp1, temp2]]; // dense array to sparse matrix
      ChangeOperator(gqPM, tempPms, parent = Ja); // send to Mat
      complex gl1 = J(uma[], Hl1);
      if(zerofreq) tempPms = [[real(gl1)]];
      else {
        bm = vM(0, XMh, tgv = -10);
        complex gw = J(uma[], bm);
        tempPms = [[real(gl1), -imag(gw)], [imag(gl1), real(gw)]];
      }
      ChangeOperator(glPM, tempPms, parent = Ja); // send to Mat
      ChangeNumbering(J, Jl2, temp1); // FreeFEM to PETSc
      yqP(0:J.n-1) = temp1;
      complex gl2 = J(uma[], Hl2);
      if(mpirank == 0) {
        yqP(yqP.n-(2-zerofreq)) = real(gl2);
        if(!zerofreq) yqP(yqP.n-1) = imag(gl2);
      }
      tempPms = [[yqP]]; // dense array to sparse matrix
      ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
      ik = 0.0;
      iomega = 0.0;
      sym = 0;
      J = vJ(XMh, XMh, tgv = -1);
      sym = sym1;
      KSPSolve(Ja, yqP, yqP);
      tempPms = [[yqP]]; // dense array to sparse matrix
      ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
      return 0;
  }
// 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 -fieldsplit_1_ksp_type preonly -fieldsplit_1_pc_type redundant -fieldsplit_1_redundant_pc_type lu", setup = 1);
set(Ja, sparams = "-prefix_push fieldsplit_0_ -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition full -fieldsplit_1_ksp_type preonly -fieldsplit_1_pc_type redundant -fieldsplit_1_redundant_pc_type lu -prefix_pop", prefix = "fieldsplit_0_", parent = Jaa, setup = 1);
set(J, IFMACRO(Jsetargs) Jsetargs, ENDIFMACRO sparams = "-prefix_push fieldsplit_0_fieldsplit_0_ " + KSPparams + " -prefix_pop", prefix = "fieldsplit_0_fieldsplit_0_", parent = Ja);
// PREDICTOR
complex[int] qa;
ChangeNumbering(J, ub[], qa);
ChangeNumbering(J, um[], qm);
ChangeNumbering(J, uma[], qma);
ChangeNumbering(J, uma[], qma, inverse = true);
if(mpirank == 0) {
  qa.resize(Jaa.n);
  qa(qa.n-(3-zerofreq)) = lambda1;
  if(!zerofreq) qa(qa.n-2) = omega;
  qa(qa.n-1) = lambda2;
}
{
  sym = 0;
  R = vR(0, XMh, tgv = -1);
  sym = sym1;
  ik.im = sym1;
  iomega = 1i*omega;
  updatelambda(param, lambda1 + eps);
  updatelambda(param2, lambda2);
  sym = 0;
  complex[int] Jl1 = vR(0, XMh, tgv = -1);
  sym = sym1;
  complex[int] Hl1 = vJ(0, XMh, tgv = -10);
  updatelambda(param, lambda1);
  updatelambda(param2, lambda2 + eps2);
  sym = 0;
  complex[int] Jl2 = vR(0, XMh, tgv = -1);
  sym = sym1;
  complex[int] Hl2 = vJ(0, XMh, tgv = -10);
  updatelambda(param2, lambda2);
  Jl1 -= R;
  Jl2 -= R;
  Jl1 /= eps;
  Jl2 /= eps2;
  complex[int] temp1(J.n), temp2(J.n), bm = vJ(0, XMh, tgv = -10);
  Hl1 -= bm;
  Hl2 -= bm;
  Hl1 /= eps;
  Hl2 /= eps2;
  ChangeNumbering(J, Jl1, temp1); // FreeFEM to PETSc
  temp1.im = 0.0;
  matrix<complex> tempPms;
  if(zerofreq) tempPms = [[temp1]];
  else tempPms = [[temp1, 0]]; // dense array to sparse matrix
  ChangeOperator(JlPM, tempPms, parent = Ja); // send to Mat
  sym = 0;
  H = vH(XMh, XMh, tgv = 0); // form the matrix (dL/dq*w)
  sym = sym1;
  MatMultHermitianTranspose(H, qma, temp1); // gqr,i
  if(!zerofreq) temp2.re = -temp1.im;
  temp1.im = 0.0;
  if(zerofreq) tempPms = [[temp1]];
  else tempPms = [[temp1, temp2]]; // dense array to sparse matrix
  ChangeOperator(gqPM, tempPms, parent = Ja); // send to Mat
  complex gl1 = J(uma[], Hl1);
  if (zerofreq) tempPms = [[real(gl1)]];
  else {
    bm = vM(0, XMh, tgv = -10);
    complex gw = J(uma[], bm);
    tempPms = [[real(gl1), -imag(gw)], [imag(gl1), real(gw)]];
  }
  ChangeOperator(glPM, tempPms, parent = Ja); // send to Mat
  ChangeNumbering(J, Jl2, temp1); // FreeFEM to PETSc
  yqP(0:J.n-1) = temp1;
  complex gl2 = J(uma[], Hl2);
  if (mpirank == 0) {
    yqP(yqP.n-(2-zerofreq)) = real(gl2);
    if(!zerofreq) yqP(yqP.n-1) = imag(gl2);
  }
  tempPms = [[yqP]]; // dense array to sparse matrix
  ChangeOperator(JlPMa, tempPms, parent = Jaa); // send to Mat
  ik = 0.0;
  iomega = 0.0;
  sym = 0;
  J = vJ(XMh, XMh, tgv = -1);
  sym = sym1;
  KSPSolve(Ja, yqP, yqP);
  tempPms = [[yqP]]; // dense array to sparse matrix
  ChangeOperator(yqPMa, tempPms, parent = Jaa); // send to Mat
  ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
  bm = vM(0, XMh, tgv = -10);
  ChangeNumbering(J, bm, qP);
  ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
  bm = vM(0, XMh, tgv = -10);
  ChangeNumbering(J, bm, pP);
  yqP0 = yqP;
  omega0 = omega;
  alpha0 = alpha[paramnames[0]];
  beta0 = beta;
}
while (!stopflag){
  complex[int] qa0 = qa;
  real h, hl = real(yqP'*yqP);
  mpiAllReduce(hl, h, mpiCommWorld, mpiSUM);
  h = h0/sqrt(h + 1.0);
  qa(0:Ja.n-1) -= (h*yqP);
  if (mpirank == 0) {
    qa(qa.n-1) += h;
    lambda1 = real(qa(qa.n-(3-zerofreq)));
    omega = zerofreq ? 0.0 : real(qa(qa.n-2)); // Extract frequency value from state vector on proc 0
    lambda2 = real(qa(qa.n-1));
  }
  broadcast(processor(0), lambda1);
  broadcast(processor(0), omega);
  broadcast(processor(0), lambda2);
  updatelambda(param, lambda1);
  updatelambda(param2, lambda2);
  // CORRECTOR LOOP
  int ret;
  it = 0;
  internalit = 0;
  adaptflag = 0;
  SNESSolve(Jaa, funcJa, funcRa, qa, 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 ((lambda1 - paramtarget)*paramdiff1 <= 0 || (lambda2 - param2target)*paramdiff2 <= 0) stopflag = true;
    h0 /= f;
    if (cosalpha < 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) {
        lambda1 = real(qa(qa.n-(3-zerofreq)));
        omega = zerofreq ? 0.0 : real(qa(qa.n-2)); // Extract frequency value from state vector on proc 0
        lambda2 = real(qa(qa.n-1));
      }
      broadcast(processor(0), lambda1);
      broadcast(processor(0), omega);
      broadcast(processor(0), lambda2);
      updatelambda(param, lambda1);
      updatelambda(param2, lambda2);
      ChangeNumbering(J, um[], qm, inverse = true);
      ChangeNumbering(J, uma[], qma, inverse = true);
      ChangeNumbering(J, yb[], yqP(0:J.n-1), inverse = true);
      complex ylambda1, yomega;
      if(mpirank == 0) {
        ylambda1 = yqP(yqP.n-(2-zerofreq));
        if(!zerofreq) yomega = yqP(yqP.n-1);
      }
      XMhg defu(uG), defu(umrG), defu(umiG), defu(umarG), defu(umaiG), defu(yG), defu(tempu);
      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) = yb[].re; // populate local portion of global soln
      mpiAllReduce(tempu[], yG[], mpiCommWorld, mpiSUM);
      if(mpirank == 0) {  // Perform mesh adaptation (serially) on processor 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);
          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);
        ENDIFMACRO
        IFMACRO(dimension,3)
      cout << "NOTE: 3D mesh adaptation is still under development." << endl;
          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));
          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
      } // TODO: add adaptation to endogeneity and structural sensitivity
      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(yb) = 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);
      um[].re = umrG[](restu);
      um[].im = umiG[](restu);
      uma[].re = umarG[](restu);
      uma[].im = umaiG[](restu);
      Mat<complex> Adapt0(J), 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
      H = Adapt0;
      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 sparams = "-prefix_push fieldsplit_0_fieldsplit_0_ " + KSPparams + " -prefix_pop", prefix = "fieldsplit_0_fieldsplit_0_", parent = Ja);
      set(Ja, sparams = "-prefix_push fieldsplit_0_ -ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition full -fieldsplit_1_ksp_type preonly -fieldsplit_1_pc_type redundant -fieldsplit_1_redundant_pc_type lu -prefix_pop", prefix = "fieldsplit_0_", parent = Jaa);
      qa.resize(J.n);
      ChangeNumbering(J, ub[], qa);
      if(mpirank == 0) {
        qa.resize(J.n+(3-zerofreq));
        qa(qa.n-(3-zerofreq)) = lambda1;
        if(!zerofreq) qa(qa.n-2) = omega;
        qa(qa.n-1) = lambda2;
      }
      qm.resize(J.n);
      ChangeNumbering(J, um[], qm);
      qma.resize(J.n);
      ChangeNumbering(J, uma[], qma);
      qP.resize(J.n);
      ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
      ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
      complex[int] Mq = vM(0, XMh, tgv = -10);
      ChangeNumbering(J, Mq, qP);
      pP.resize(J.n);
      ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
      Mq = vM(0, XMh, tgv = -10);
      ChangeNumbering(J, Mq, pP);
      R.resize(ub[].n);
      yqP.resize(J.n);
      ChangeNumbering(J, yb[], yqP);
      if (mpirank==0) {
        yqP.resize(Ja.n);
        yqP(yqP.n-(2-zerofreq)) = ylambda1;
        if(!zerofreq) yqP(yqP.n-1) = yomega;
      }
      yqP0.resize(Ja.n);
      yqP0 = yqP;
      qa0.resize(qa.n);
      qap.resize(qa.n);
      it = 0;
      internalit = 0;
      adaptflag = 1;
      SNESSolve(Jaa, funcJa, funcRa, qa, 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) {
      lambda1 = real(qa(qa.n-(3-zerofreq)));
      omega = zerofreq ? 0.0 : real(qa(qa.n-2)); // Extract frequency value from state vector on proc 0
      lambda2 = real(qa(qa.n-1));
    }
    broadcast(processor(0), lambda1);
    broadcast(processor(0), omega);
    broadcast(processor(0), lambda2);
    updatelambda(param, lambda1);
    updatelambda(param2, lambda2);
    ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
    complex[int] Mq = vM(0, XMh, tgv = -10);
    ChangeNumbering(J, um[], qm, inverse = true);
    ChangeNumbering(J, uma[], qma, inverse = true);
    real Mnorm = sqrt(real(J(um[], Mq)));
    um[] /= Mnorm; // so that <um[],M*um[]> = 1
    uma[] *= (Mnorm/J(Mq, uma[])); // so that <uma[],M*um[]> = 1
    ChangeNumbering(J, um[], qm);
    ChangeNumbering(J, uma[], qma);
    if (normalform){
      complex[int] qAAs(J.n), qAA(J.n), tempP(J.n), temp(um[].n), temp1(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 = -1);
      if(paramnames[0] != ""){
        temp1 = vR(0, XMh, tgv = -1);
        for (int k = 0; k < paramnames.n; ++k){
          lambda1 = getlambda(paramnames[k]);
          updatelambda(paramnames[k], lambda1 + eps);
          temp = vR(0, XMh, tgv = -1);
          updatelambda(paramnames[k], lambda1);
          temp -= temp1;
          temp /= -eps;
          ChangeNumbering(J, temp, tempP); // FreeFEM to PETSc
          KSPSolve(J, tempP, qAAs);
          qDa(k, :) = qAAs;
        }
      }
      //  B: base modification due to quadratic nonlinear interaction
      ik.im = sym1;
      ik2.im = -sym1;
      iomega = 1i*omega;
      iomega2 = -1i*omega;
      ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
      um2[] = conj(um[]);
      temp = vH(0, XMh, tgv = -10);
      temp.re *= -1.0; // -2.0/2.0
      temp.im = 0.0;
      ChangeNumbering(J, temp, tempP); // FreeFEM to PETSc
      KSPSolve(J, tempP, qAAs);
      //  C: harmonic generation due to quadratic nonlinear interaction
      ik.im = sym1;
      ik2.im = sym1;
      iomega = 1i*omega;
      iomega2 = 1i*omega;
      um2[] = um[];
      sym = 2*sym1;
      temp = vH(0, XMh, tgv = -10);
      temp *= -0.5; // -1.0/2.0
      ChangeNumbering(J, temp, tempP); // FreeFEM to PETSc
      ik.im = 2*sym1;
      iomega = 2i*omega;
      J = vJ(XMh, XMh, tgv = -1);
      KSPSolve(J, tempP, qAA);
      // 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] != ""){
        temp1 = vJ(0, XMh, tgv = -10);
        for (int k = 0; k < paramnames.n; ++k){
          lambda1 = getlambda(paramnames[k]);
          updatelambda(paramnames[k], lambda1 + eps);
          temp = vJ(0, XMh, tgv = -10);
          updatelambda(paramnames[k], lambda1);
          temp -= temp1;
          temp /= eps;
          ChangeNumbering(J, um2[], qDa(k, :), inverse = true, exchange = true); // FreeFEM to PETSc
          temp += vH(0, XMh, tgv = -10); // 2.0/2.0
          alpha[paramnames[k]] = -J(uma[], temp);
        }
      }
      //  B: fundamental modification due to cubic self-interaction of fundamental
      ik.im = sym1;
      ik2.im = sym1;
      ik3.im = -sym1;
      iomega = 1i*omega;
      iomega2 = 1i*omega;
      iomega3 = -1i*omega;
      um2[] = um[];
      um3[] = conj(um[]);
      temp = vT(0, XMh, tgv = -10);
      temp *= 0.5; //3.0/6.0
      //  C: fundamental modification due to quadratic interaction of fundamental with 2nd order modification B
      ik.im = sym1;
      ik2 = 0.0;
      iomega = 1i*omega;
      iomega2 = 0.0;
      ChangeNumbering(J, um2[], qAAs, inverse = true, exchange = true); // FreeFEM to PETSc
      temp += vH(0, XMh, tgv = -10);
      //  D: fundamental modification due to quadratic interaction of fundamental with 2nd order modification C
      ik.im = -sym1;
      ik2.im = 2*sym1;
      iomega = -1i*omega;
      iomega2 = 2i*omega;
      um[] = conj(um[]);
      ChangeNumbering(J, um2[], qAA, inverse = true, exchange = true); // FreeFEM to PETSc
      temp += vH(0, XMh, tgv = -10);
      beta = -J(uma[], temp);
      ik2 = 0.0;
      ik3 = 0.0;
      iomega2 = 0.0;
      iomega3 = 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][], qAAs, inverse = true); // FreeFEM to PETSc
        savemode(fileout + "_" + count + "_wnl_AAs", "", fileout + ".hopf", meshout, vec, val, sym, true);
        ChangeNumbering(J, vec[0][], qAA, 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) {
      if(mpirank == 0) {
        if(real(alpha0)*real(alpha[paramnames[0]]) < 0) cout << "\tFold-Hopf bifurcation detected." << endl;
        else cout << "\tBautin bifurcation detected." << endl;
      }
      forcesave = true;
    }
    if (omega*omega0 < 0){
      if(mpirank == 0) cout << "\tBogdanov-Takens bifurcation (or zero-frequency point) detected." << endl;
      forcesave = true;
    }
    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);
    forcesave = false;
    ChangeNumbering(J, ub[], qa(0:J.n-1), inverse = true, exchange = true);
    ChangeNumbering(J, um[], qm, inverse = true, exchange = true);
    Mq = vM(0, XMh, tgv = -10);
    ChangeNumbering(J, Mq, qP);
    ChangeNumbering(J, um[], qma, inverse = true, exchange = true);
    Mq = vM(0, XMh, tgv = -10);
    ChangeNumbering(J, Mq, pP);
    yqP0 = yqP;
    omega0 = omega;
    alpha0 = alpha[paramnames[0]];
    beta0 = beta;
  }
  else {
    if(mpirank == 0){
      if(res*(monotone!=0) >= resp) cout << "\tResidual norm failed to decrease. Reattempting with smaller step." << endl;
      if(kappa >= kappamax) cout << "\tContraction rate exceeds " << kappamax << ". Reattempting with smaller step." << endl;
      if(it >= snesmaxit) cout << "\tFailed to converge within limit of " + snesmaxit + " iterations. Reattempting with smaller step." << endl;
      if(maxdelta >= deltamax) cout << "\tStep size exceeds " << deltamax << "." << endl;
      if(acos(abs(cosalpha)) >= anglemax) cout << "\tAngle exceeds " << (anglemax*180./pi) << " degrees." << endl;
    }
    h0 /= fmax;
    qa = qa0;
  }
}