// // foldcompute.edp // Chris Douglas // cdouglas33@gatech.edu // load "PETSc" macro dimension()2//EOM include "macro_ddm.idp" int[int] n2o, rest; macro ThN2O()n2o// EOM mesh Th = square(100,10,[10*x,y-0.5]); mesh Thg = Th; macro def(u)[u,u#y,u#p]//EOM macro init(i)[i, i, i]//EOM func Pk = [P2,P2,P1];// Define finite element space buildDmesh(Th); fespace XMh(Th, Pk); fespace XMhg(Thg, Pk); rest = restrict(XMh, XMhg, n2o); XMh def(dub), def(um), def(ub); Mat J, H, A; createMat(Th, J, Pk); createMat(Th, H, Pk); real nu = 0.1; macro div(u) ( dx(u) + dy(u#y) )//EOM macro ugradu(v, U, u) ( v*(U*dx(u) + U#y*dy(u)) + v#y*(U*dx(u#y) + U#y*dy(u#y)))//EOM macro visc(v, u) ( dx(v)*dx(u) + dy(v)*dy(u) + dx(v#y)*dx(u#y) + dy(v#y)*dy(u#y) )// EOM // Residual varf vR(def(dub), def(v)) = int2d(Th)( ugradu(v, ub, ub) - ubp*div(v) + nu*visc(v, ub) - vp*div(ub) ) + on(1, dub = ub - (2.0 - 8.0*y^2), duby = uby) + on(2,4, dub = ub, duby = uby); // Jacobian varf vJ(def(dub), def(v)) = int2d(Th)( ugradu(v, dub, ub) + ugradu(v, ub, dub) - dubp*div(v) + nu*visc(v, dub) - vp*div(dub) ) + on(1, dub = ub - (2.0 - 8.0*y^2), duby = uby) + on(2,4, dub = ub, duby = uby); varf vJq(def(dum), def(v)) = int2d(Th)( ugradu(v, um, ub) + ugradu(v, ub, um) - ump*div(v) + nu*visc(v, um) - vp*div(um) ) + on(1, 2, 4, dum = 0, dumy = 0); // Hessian varf vH(def(dub), def(v)) = int2d(Th)( ugradu(v, um, dub) + ugradu(v, dub, um) ) + on(1, 2, 4, dub = 0, duby = 0); varf vHq(def(dum), def(v)) = int2d(Th)( ugradu(v, um, dub) + ugradu(v, dub, um) ) + on(1, 2, 4, dum = 0, dumy = 0); // Parameter gradients varf vJl(def(dub), def(v)) = int2d(Th)( visc(v, ub) ) + on(1, 2, 4, dub = 0, duby = 0); varf vHl(def(dub), def(v)) = int2d(Th)( visc(v, um) ) + on(1, 2, 4, dub = 0, duby = 0); set(J, sparams = " -ksp_type preonly -pc_type lu "); real[int] q; changeNumbering([J, J], [ub[], um[]], q); if(mpirank == 0) q.resize(q.n+1); A = [[J,0,0],[0,J,0],[0,0,1]]; func real[int] funcResidual(real[int]& inPETSc) { changeNumbering([J, J], [ub[], um[]], inPETSc(0:inPETSc.n - (mpirank == 0 ? 2 : 1)), inverse = true, exchange = true); real[int] out1 = vR(0, XMh, tgv = -1); real[int] out2 = vJq(0, XMh, tgv = -1); real out3 = J(um[],um[]) - 1.0; real[int] outPETSc; changeNumbering([J, J], [out1, out2], outPETSc); if(mpirank == 0) { outPETSc.resize(outPETSc.n+1); outPETSc(outPETSc.n-1) = out3; } return outPETSc; } func int funcJacobian(real[int]& inPETSc) { changeNumbering([J, J], [ub[],um[]], inPETSc(0:inPETSc.n - (mpirank == 0 ? 2 : 1)), inverse = true, exchange = true); J = vJ(XMh, XMh, tgv = -1); H = vH(XMh, XMh, tgv = -10); real[int] Jl = vJl(0, XMh, tgv = -10); real[int] Hl = vHl(0, XMh, tgv = -10); real[int] qT; changeNumbering(J, um[], qT); A = [[J,0,Jl],[H,J,Hl],[0,qT',0]]; return 0; } set(A, sparams = "-pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_detect_saddle_point -ksp_view -ksp_monitor"); string SNESparams = " -snes_monitor -snes_converged_reason "; SNESSolve(A, funcJacobian, funcResidual, q, sparams = SNESparams); // solve nonlinear problem with SNES