minimal-surface-Tao-2d-PETSc.edp

//  run with MPI:  ff-mpirun -np 4 script.edp
// NBPROC 4

load "PETSc"
macro dimension()2//
include "macro_ddm.idp"

macro def(i)i// EOM
macro init(i)i// EOM
func Pk = P2;
macro grad(u)[dx(u), dy(u)]// EOM

mesh ThNo, Th = square(getARGV("-global", 20), getARGV("-global", 20), [1 + x, y]);
fespace Vh(Th, Pk);
Mat H;
MatCreate(Th, H, Pk);
{
    fespace Ph(Th, P0);
    Ph part;
    PartitionCreate(Th, part[], P0);
    ThNo = trunc(Th, abs(part - 1.0) < 1e-1);
}
func g = cos(pi*x)*cos(pi*y);
func real J(real[int]& X) {
    Vh u;
    ChangeNumbering(H, u[], X, inverse = true, exchange = true);
    real glob, loc = int2d(ThNo)(sqrt(1 + grad(u)'*grad(u)));
    mpiAllReduce(loc, glob, mpiCommWorld, mpiSUM);
    return glob;
}

func real[int] DJ(real[int]& X) {
    Vh u;
    ChangeNumbering(H, u[], X, inverse = true, exchange = true);
    varf vG(w, v) = int2d(Th)((grad(u)'*grad(v)) / sqrt(1 + grad(u)'*grad(u)));
    real[int] out = vG(0, Vh);
    real[int] outPETSc;
    ChangeNumbering(H, out, outPETSc);
    return outPETSc;
}

func int HJ(real[int]& X) {
    Vh u;
    ChangeNumbering(H, u[], X, inverse = true, exchange = true);
    varf vH(v, w) = int2d(Th)((grad(w)'*grad(v)) / sqrt(1 + grad(u)'*grad(u))
            - (grad(w)'*grad(u)) * (grad(v)'*grad(u)) *(1 + grad(u)'*grad(u))^-1.5);
    matrix Loc = vH(Vh, Vh);
    H = Loc;
    return 0;
}
varf onGamma(u, v) = on(1, 2, 3, 4, u = 1);
Vh onG;
onG[] = onGamma(0, Vh, tgv = 1);
real[int] lbPETSc, ubPETSc;
{
    Vh lb = onG != 0 ? g : -1e19;
    Vh ub = onG != 0 ? g :  1e19;
    lb = max(lb, 3-100*(((x-1.5)^2 + (y-0.5)^2))^2);
    ChangeNumbering(H, lb[], lbPETSc);
    ChangeNumbering(H, ub[], ubPETSc);
}

{
    Vh u = onG != 0 ? g : 0;
    real[int] uPETSc;
    ChangeNumbering(H, u[], uPETSc);
    TaoSolve(H, J, DJ, uPETSc, xl = lbPETSc, xu = ubPETSc, sparams = "-tao_monitor -tao_view -tao_type bqnls -tao_max_it 40 -tao_gatol 1e-4");
    ChangeNumbering(H, u[], uPETSc, inverse = true, exchange = true);
    real Ju = J(uPETSc);
    plotMPI(Th, u, Pk, def, real, cmm = "Global solution, J(u) = " + Ju);
}
{
    Vh u = onG != 0 ? g : 0;
    real[int] uPETSc;
    ChangeNumbering(H, u[], uPETSc);
    TaoSolve(H, J, DJ, uPETSc, xl = lbPETSc, xu = ubPETSc, sparams = "-tao_monitor -tao_view -tao_type bnls -tao_max_it 40 -tao_gatol 1e-4", HessianRoutine = HJ);
    ChangeNumbering(H, u[], uPETSc, inverse = true, exchange = true);
    real Ju = J(uPETSc);
    plotMPI(Th, u, Pk, def, real, cmm = "Global solution, J(u) = " + Ju);
}