// run with MPI: ff-mpirun -np 4 script.edp // NBPROC 4 load "PETSc" macro dimension()2// include "macro_ddm.idp" ////////////////////////////// macro deff(i)i//EOM ////////////////////////////////////// func Pk = [P2,P2,P1]; ////////////////////////////////////////// real nu=0.1, T=1., dt = 0.1; int m, M= T/dt; mesh ThNo; // Here ThNo is the domain decomposition without overlapping used in J mesh Th = square(getARGV("-global", 50), getARGV("-global", 50)); savemesh(Th, "meshTh.msh"); plot(Th,wait=1,cmm="mesh Th before creatMat"); /////////////////////////////////////////////////////// fespace Vh(Th,Pk); Mat H; { macro def(i)[i,i#BB,i#PP]// EOM macro init(i)[i,i,i]// EOM createMat(Th, H, Pk); } plot(Th,wait=1,cmm="mesh Th after creatMat"); fespace XXh(Th,P2), Qh(Th,P1); { fespace Ph(Th, P0); Ph part; plot(part, fill=1, value=1,wait=1,cmm="part befor creatPartition"); createPartition(Th, part[], P0) plot(part, fill=1, value=1,wait=1,cmm="part after creatPartition"); ThNo = trunc(Th, abs(part - 1.0) < 1e-1); plot(ThNo,wait=1,cmm="mesh ThNo"); } Vh [w1,w2,wp]=[0,0,0]; // used as the initial value for the first time step XXh uold=w1,vold=w2; func real J(real[int]& X) { Vh [uAA,uBB,uPP]; changeNumbering(H, uAA[], X, inverse = true, exchange = true); real glob; real loc = int2d(ThNo)( 0.5*(uAA^2 + uBB^2)/dt + 0.5*nu*(dx(uAA)*dx(uAA) + dy(uAA)*dy(uAA) + dx(uBB)*dx(uBB) + dy(uBB)*dy(uBB)) - 0.5*(uPP^2)*1.e-1 - (uold*uAA+vold*uBB)/dt); mpiAllReduce(loc, glob, mpiCommWorld, mpiSUM); return glob; } func real[int] DJ(real[int]& X) { Vh [uAA,uBB,uPP]; changeNumbering(H, uAA[], X, inverse = true, exchange = true); varf aa ([uk,vk,pk],[uu,vv,pp]) = int2d(Th)( (uAA*uu+uBB*vv)/dt + nu*(dx(uAA)*dx(uu) + dy(uAA)*dy(uu) + dx(uBB)*dx(vv) + dy(uBB)*dy(vv)) - uPP*pp*1.e-1 ) - int2d(Th)((uold*uu+vold*vv)/dt); real[int] G= aa(0,Vh); real[int] outPETSc; changeNumbering(H, G, outPETSc); return outPETSc; } func int HJ(real[int]& X) { Vh [uAA,uBB,uPP]; changeNumbering(H, uAA[], X, inverse = true, exchange = true); varf AK ([uu,vv,pp],[du,dv,dp]) = int2d(Th)( (du*uu+dv*vv)/dt + nu*(dx(du)*dx(uu) + dy(du)*dy(uu) + dx(dv)*dx(vv) + dy(dv)*dy(vv)) - dp*pp*1.e-1 ); matrix Loc = AK(Vh, Vh); H = Loc; return 0; } varf vGamma1(u,uu) = on(1,2,4,u=1); varf vGamma3(u,uu) = on(3,u=1); real[int] onGamma1 = vGamma1(0,XXh); real[int] onGamma3 = vGamma3(0,XXh); real[int] lbPETSc, ubPETSc; { XXh lbu=-1.e19; XXh lbv=-1.e19; Qh lbp=-1.e19; XXh ubu=1.e19; XXh ubv=1.e19; Qh ubp=1.e19; lbu[] = onGamma1 ? 0.:lbu[]; // for u lbv[] = onGamma1 ? 0.:lbv[]; // for v ubu[] = onGamma1 ? 0.:ubu[]; // for u ubv[] = onGamma1 ? 0.:ubv[]; // for v lbu[] = onGamma3 ? 1.:lbu[]; // for u lbv[] = onGamma3 ? 2.7:lbv[]; // for v ubu[] = onGamma3 ? 1.:ubu[]; // for u ubv[] = onGamma3 ? 2.7:ubv[]; // for v XXh iniu=0.0; XXh iniv=0.0; Qh inip=0.0; iniu[] = onGamma1 ? 0.:iniu[]; // for u iniv[] = onGamma1 ? 0.:iniv[]; // for v iniu[] = onGamma3 ? 1.:iniu[]; // for u iniv[] = onGamma3 ? 2.7:iniv[]; // for v Vh [lb1,lb2,lb3]=[lbu,lbv,lbp]; Vh [ub1,ub2,ub3]=[ubu,ubv,ubp]; [w1,w2,wp]=[iniu,iniv,inip]; changeNumbering(H, lb1[], lbPETSc); changeNumbering(H, ub1[], ubPETSc); } /////////////////////////////////////////// real[int] uPETSc; XXh u,v; Qh p; for(m=0;m