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

load "PETSc"                        // PETSc plugin
include "macro_ddm.idp"             // additional DDM functions

func Pk = P1;                       // finite element space
func Pkdc = P1dc;                   // finite element space

mesh Th = square(getARGV("-global", 40), getARGV("-global", 40));
Mat A;
macro ThRefinementFactor()getARGV("-split", 1)//
MatCreate(Th, A, Pk);

fespace Vh(Th, Pk);
fespace Ph(Th, Pkdc);
Vh b = 1;
real a = 0.001;
func real f(real u) { return u*a+u-log(1+u); }
func real df(real u) { return a+u/(1+u); }
func real ddf(real u) { return 1/((1+u)*(1+u)); }

Ph alpha;
Ph dfalpha;
Ph ddfalpha;
Vh u = 0;

varf vJ(uh, vh) = int2d(Th)(dfalpha*(dx(u)*dx(vh) + dy(u)*dy(vh)))
    + on(1, 2, 3, 4, uh = 0);
varf vC(uh, vh) = int2d(Th)(b*vh)
    + on(1, 2, 3, 4, uh = 0);
varf vhJ(uh, vh) = int2d(Th)(dfalpha*(dx(uh)*dx(vh) + dy(uh)*dy(vh)) + ddfalpha*(dx(u)*dx(vh) + dy(u)*dy(vh)) * (dx(u)*dx(uh) + dy(u)*dy(uh)))
    + on(1, 2, 3, 4, uh = 0);

func real[int] funcRes(real[int]& inPETSc) {
    ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
    alpha = dx(u)*dx(u) + dy(u)*dy(u);
    dfalpha = df(alpha);
    ddfalpha = 2 * ddf(alpha);
    real[int] v = vJ(0, Vh, tgv = -1);
    real[int] outPETSc;
    ChangeNumbering(A, v, outPETSc);
    return outPETSc;
}
func int funcJ(real[int]& inPETSc) {
    ChangeNumbering(A, u[], inPETSc, inverse = true, exchange = true);
    A = vhJ(Vh, Vh, tgv = -1);
    return 0;
}

Vh xu;
xu = sqrt((x - 0.5)^2 + (y - 0.5)^2) < 0.4 ? 0.2 : 1.0;
macro def(u)u//
plotMPI(Th, xu, Pk, def, real, cmm = "Upper bound");
real[int] xPETSc, bPETSc, xuPETSc;
u[] = vC(0, Vh, tgv = -1);
ChangeNumbering(A, u[], bPETSc);
ChangeNumbering(A, xu[], xuPETSc);
xPETSc.resize(bPETSc.n);
xPETSc = 0.1;
SNESSolve(A, funcJ, funcRes, bPETSc, xPETSc, xu = xuPETSc, sparams = "-snes_monitor -ksp_converged_reason -snes_view -snes_vi_monitor -snes_type vinewtonrsls -snes_rtol 1.0e-6 -pc_type lu");
ChangeNumbering(A, u[], xPETSc, inverse = true, exchange = true);
plotMPI(Th, u, Pk, def, real, cmm = "Global solution");
