update parallel steady codes on reference domain

Author: DongweiYe

Example1-Stenosis/ReferNS.edp renamed to Example1-Stenosis/ReferNS_steady_parallel.edp

@@ -1,4 +1,9 @@
-for(int sample = 0; sample<100; sample++){
+load "PETSc"
+macro dimension()2//
+include "macro_ddm.idp"
+
+
+for(int sample = 0; sample<500; sample++){
     cout<<"Sample:"+sample<<endl;
 
     /////// Import weights for each CP
@@ -78,36 +83,54 @@ for(int sample = 0; sample<100; sample++){
 
 
     mesh Mesh = readmesh("data/reference/reference.msh");
+    mesh MeshBackup = readmesh("data/reference/reference.msh");
 
     int[int] orderOut = [1, 1, 1, 1];
 
-    real umax = 2.51327412286999996*1e3 / (4 * pi); //mm/s
-
-    // real uin = -umax*(y-1)^2+umax;
+    real umax = 2.51327412286999996*1e3 / (8 * 4 * pi); //mm/s
     real niu = 3.5*1e-3; // g mm^-1 s^-1
     real rho = 1.06*1e-3; // g mm^-3
     real nu = niu / rho;
-    
-    
-    cout << "Number of Elements: " + Mesh.nt << endl;
+
+
+
+    // fespace Uh(Th, P2);
+    // fespace Vh(Th, P1);
+    // Uh u, v, uh, vh, uold, vold, du, dv;
+    // Vh p, ph, dp, pold;
+
+    if (mpirank == 0)
+      cout << "Number of Elements: " + Mesh.nt << endl;
 
     func PkVector = [P2, P2, P1];
 
+    int[int] myN2O;
+    macro MeshN2O() myN2O//
+
+    buildDmesh(Mesh);
+    Mat J;
+
+    {
+      macro def(i)[i, i#B, i#C]//
+      macro init(i)[i, i, i]//
+      createMat(Mesh, J, PkVector)
+    }
+
     fespace SpaceVector(Mesh, PkVector);
     fespace Uh(Mesh,P2);
 
-
-    SpaceVector [ucx, ucy, pc];
-    SpaceVector [ux, uy, p];
-    SpaceVector [vx, vy, q];
+    SpaceVector [ucx, ucy, pc] = [-umax*(y-1)^2+umax, 0, 0];
+ 
     Uh FX = funcT(weights,cps,int(cps.n/2))[0];
     Uh FY = funcT(weights,cps,int(cps.n/2))[1];
 
-    cout << "Finite Element DOF (in each partition): " + SpaceVector.ndof << endl;
 
-    problem vRes([ux, uy, p], [vx, vy, q]) = int2d(Mesh)(
+    if (mpirank == 0)
+      cout << "Finite Element DOF (in each partition): " + SpaceVector.ndof << endl;
 
-        + JF(FX,FY)*nu*((dx(ucx)*Finv11(FX,FY)+dy(ucx)*Finv21(FX,FY))*(dx(vx)*Finv11(FX,FY)+dy(vx)*Finv21(FX,FY))
+
+    varf vRes([ux, uy, p], [vx, vy, q]) = int2d(Mesh)(
+        JF(FX,FY)*nu*((dx(ucx)*Finv11(FX,FY)+dy(ucx)*Finv21(FX,FY))*(dx(vx)*Finv11(FX,FY)+dy(vx)*Finv21(FX,FY))
                         +(dx(ucx)*Finv12(FX,FY)+dy(ucx)*Finv22(FX,FY))*(dx(vx)*Finv12(FX,FY)+dy(vx)*Finv22(FX,FY)) 
                         +(dx(ucy)*Finv11(FX,FY)+dy(ucy)*Finv21(FX,FY))*(dx(vy)*Finv11(FX,FY)+dy(vy)*Finv21(FX,FY))
                         +(dx(ucy)*Finv12(FX,FY)+dy(ucy)*Finv22(FX,FY))*(dx(vy)*Finv12(FX,FY)+dy(vy)*Finv22(FX,FY)))
@@ -118,12 +141,11 @@ for(int sample = 0; sample<100; sample++){
         - JF(FX,FY)/rho*pc*((dx(vx)*Finv11(FX,FY)+dy(vx)*Finv21(FX,FY))+(dx(vy)*Finv12(FX,FY)+dy(vy)*Finv22(FX,FY))) 
         - JF(FX,FY)*q*((dx(ucx)*Finv11(FX,FY)+dy(ucx)*Finv21(FX,FY))+(dx(ucy)*Finv12(FX,FY)+dy(ucy)*Finv22(FX,FY))))
 
-        +on(1,3,ux = 0, uy = 0)
-        +on(4, ux = -umax*(y-1)^2+umax, uy = 0);
-
-    problem vJ([ux, uy, p], [vx, vy, q]) = int2d(Mesh)(
+        +on(1,3,ux = ucx-0, uy = ucy-0)
+        +on(4, ux = 0, uy = 0);
 
-        + JF(FX,FY)*vx*(ucx*(dx(ux)*Finv11(FX,FY)+dy(ux)*Finv21(FX,FY))
+    varf vJ([ux, uy, p], [vx, vy, q]) = int2d(Mesh)(
+        JF(FX,FY)*vx*(ucx*(dx(ux)*Finv11(FX,FY)+dy(ux)*Finv21(FX,FY))
                             +ux*(dx(ucx)*Finv11(FX,FY)+dy(ucx)*Finv21(FX,FY))
                             +ucy*(dx(ux)*Finv12(FX,FY)+dy(ux)*Finv22(FX,FY))
                             +uy*(dx(ucx)*Finv12(FX,FY)+dy(ucx)*Finv22(FX,FY))) 
@@ -141,33 +163,80 @@ for(int sample = 0; sample<100; sample++){
         - JF(FX,FY)*q*((dx(ux)*Finv11(FX,FY)+dy(ux)*Finv21(FX,FY))+(dx(uy)*Finv12(FX,FY)+dy(uy)*Finv22(FX,FY))))
 
 
-        + on(1,3, ux = 0, uy = 0)
-        + on(4, ux = -umax*(y-1)^2+umax , uy = 0);
+        + on(1,3, ux = ucx-0, uy = ucy-0)
+        + on(4, ux = 0 , uy = 0);
+
+    set(J, sparams = "-pc_type lu");
+
 
 
-    real arrns = 1e-9;
+    func real[int] funcRes(real[int]& inPETSc) {
+        ChangeNumbering(J, ucx[], inPETSc, inverse = true, exchange = true);
+        real[int] out(SpaceVector.ndof);
+        out = vRes(0, SpaceVector, tgv = -1);
+        ChangeNumbering(J, out, inPETSc);
+        return inPETSc;
+    }
+
+    func int funcJ(real[int]& inPETSc) {
+        ChangeNumbering(J, ucx[], inPETSc, inverse = true, exchange = true);
+        J = vJ(SpaceVector, SpaceVector, tgv = -1);
+        return 0;
+    }
+
+
+    // Save to global solution
+    fespace SpaceVectorGlobal(MeshBackup, PkVector);
+    int[int] rst = restrict(SpaceVector, SpaceVectorGlobal, myN2O);
+    SpaceVectorGlobal [globux, globuy, globp], [sumux, sumuy, sump]; 
+    SpaceVector [uxTemp, uyTemp, pTemp];
+    fespace SpaceP1Global(MeshBackup, P1);
+    SpaceP1Global uxi, uyi, pi;
+
+
+
+    for(int i = 0; i<3;i++){
+        ucx[] = ucx[]*2;
+        real[int] xPETSc;
+        ChangeNumbering(J, ucx[], xPETSc);
+        SNESSolve(J, funcJ, funcRes, xPETSc, sparams = "-snes_monitor ");
+        ChangeNumbering(J, ucx[], xPETSc, inverse = true, exchange = true);
+    }
+
+    globux[] = 0;
+    globuy[] = 0;
+    globp[] = 0;
+
+    [uxTemp, uyTemp, pTemp] = [ucx, ucy, pc];
+
+    uxTemp[] .*= J.D;
+    uyTemp[] .*= J.D;
+    pTemp[] .*= J.D;
+
+    for[i, v : rst] globux[][v] = uxTemp[][i];
+    for[i, v : rst] globuy[][v] = uyTemp[][i];
+    for[i, v : rst] globp[][v] = pTemp[][i];
+
+    mpiAllReduce(globux[], sumux[], mpiCommWorld, mpiSUM);
+    mpiAllReduce(globuy[], sumuy[], mpiCommWorld, mpiSUM);
+    mpiAllReduce(globp[], sump[], mpiCommWorld, mpiSUM);
+
+    uxi = sumux;
+    uyi = sumuy;
+    pi = sump;
+
+    // if(mpirank == 0){
+    //     savevtk("result/output_"+sample+".vtu", MeshBackup, uxi, uyi, pi, dataname="u v p", order=orderOut);
+    // }
 
-    macro ns() {
-        int n;
-        real err=0;
-        vRes;
-
-        for(n=0; n< 15; n++) {
-            vJ;
-
-            ucx[] -= ux[];
-            ucy[] -= uy[];
-            pc[] -= p[];
-            real Lu1=ucx[].linfty, Lu2=ucy[].linfty, Lp=pc[].linfty;
-            err = ux[].linfty/Lu1 + uy[].linfty/Lu2 + p[].linfty/Lp;
-            cout << n << " err = " << err << " " << arrns << " rey = " << 1./nu << endl;
-            if(err < arrns) break; 
-            }
-    } //EOF
-    ns;
-
-    // Vectorial [ux,uy,p] savev the data together
-    // ofstream uxfile("data/snapshots/sample_"+sample+"/uvp.txt");
-    // uxfile << ux[] << endl;
+    if(mpirank == 0){
+        ofstream uxfile("data/snapshots/sample_"+sample+"/u.txt");
+        uxfile << uxi[] << endl;
+        ofstream uyfile("data/snapshots/sample_"+sample+"/v.txt");
+        uyfile << uyi[] << endl;
+        ofstream pfile("data/snapshots/sample_"+sample+"/p.txt");
+        pfile << pi[] << endl;
+    }
+
 
 }