Hello, Pierre! In the previous post I performed an unsteady three-dimensional laminar flow in parallel. I lowered the time step very much to see if I still get the optimizing error. Sorry for the bad crop of the printscreen. I plotted the velocity along the geometry (a simple pipe). As you can see I chose the color scale for the velocity values in white-black but the geometry has the color red. There is no red on the color-scale. Then I plotted the velocity values along a line which is located inside the pipe, from the inlet to the outlet, but the graph shows no values. The strange part is that in the color-scale I have some values which I expect but on the geometry (even if I cut a slice or something) the velocity is shown with “red” there is no flow developing etc. so I believe there is an error in the visualization which I don’t know how to fix.

Now, I’m trying to run in parallel the same flow but in a steady regime to avoid `convect`

. I used these examples Ex.1 and mostly Ex.2 but I am missing something as I get a flow equal to 0. I attach below the code I wrote.

load “PETSc”

macro dimension()3

macro def(j)[j, j#B, j#C, j#D]

macro init(j)[j, j, j, j]

include “macro_ddm.idp”

load “iovtk”

load “msh3”

load “medit”

mesh3 Th = readmesh3 …

func FES=[P2,P2,P2,P1];

fespace XXXXFES(Th,FES);

mesh3 ThBackup=Th;

fespace XXXXFESBackup(ThBackup, FES);

XXXXFES[u,v,w,p];

XXXXFES[hx,hy,hz,q];

XXXXFES[up,vp,wp,pp];

XXXXFESBackup [ub, vb, wb, pb];

XXXXFESBackup [ured, vred, wred, pred];

XXXXFES[du,dv,dw,dp];

XXXXFES[ug,vg,wg,pg];

int[int] n2o;

macro ThN2O()n2o

Mat NS;

createMat(Th, NS, FES);

real InletVelW = 1;

int iter = 2; //Numer of iterations

real Re = 500.; //Reynolds number

[ug, vg, wg, pg]=[0., 0., 1., 0.];

varf vRes([du,dv,dw,dp],[hx,hy,hz,q])

=int3d(Th) (

(ug * dx(ug)+vg * dy(ug)+wg * dz(ug)) * hx

+(ug*dx(vg)+vg*dy(vg)+wg*dz(vg)) * hy*

+(ugdx(wg)+vg*dy(wg)+wg*dz(wg)) * hz

+(1/Re)*(dx(ug)*dx(hx)+dy(ug)*dy(hx)+dz(ug) * dz(hx)) *

+(1/Re)(dx(vg)*dx(hy)+dy(vg)*dy(hy)+dz(vg) * dz(hy)) *

+(1/Re)(dx(wg)*dx(hz)+dy(wg)*dy(hz)+dz(wg) * dz(hz))*

- pg(dx(hx)+dy(hy)+dz(hz))

-(dx(ug)+dy(vg)+dz(wg))*q

)

+on(Wall, du=ug-0, dv=vg-0, dw=wg-0)

+on(Inlet, du=ug-0, dv=vg-0, dw=wg-1.);

varf vJ([du,dv,dw,dp],[hx,hy,hz,q])

=int3d(Th) (

(du*dx(ug)+dv*dy(ug)+dw*dz(ug)) * hx *

+(ugdx(du)+vg*dy(du)+wg*dz(du)) * hx

+(du*dx(vg)+dv*dy(vg)+dw*dz(vg)) * hy *

+(ugdx(dv)+vg*dy(dv)+wg*dz(dv)) * hy

+(du*dx(wg)+dv*dy(wg)+dw*dz(wg)) * hy *

+(ugdx(dw)+vg*dy(dw)+wg*dz(dw)) * hy

+(1/Re)*(dx(du)*dx(hx)+dy(du)*dy(hx)+dz(du) * dz(hx)) *

+(1/Re)(dx(dv)*dx(hy)+dy(dv)*dy(hy)+dz(dv) * dz(hy))*

+(1/Re)(dx(dw)*dx(hz)+dy(dw)**dy(hz)+dz(dw) * dz(hz)) *

-dp(dx(hx)+dy(hy)+dz(hz))

-(dx(du)+dy(dv)+dz(dw)) * q

)

+on(Wall, du=ug-0, dv=vg-0, dw=wg-0)

+on(Inlet, du=ug-0, dv=vg-0, dw=wg-1.);

set(NS, sparams = “-pc_type lu”);

func real[int] funcRes(real[int]&inPETSc){

ChangeNumbering(NS, ug[], inPETSc, inverse = true, exchange = true);

real[int] out(XXXXFES.ndof);

out = vRes(0, XXXXFES, tgv = -1);

ChangeNumbering(NS, out, inPETSc);

return inPETSc;

}

func int funcJ(real[int]& inPETSc){

ChangeNumbering(NS, ug[], inPETSc, inverse = true, exchange = true);

NS = vJ(XXXXFES, XXXXFES, tgv = -1);

return 0;

}

fespace Qh(Th, P1);

for (int i=0; i<=iter; i++)

{

real[int] xPETSc;

ChangeNumbering(NS, ug[], xPETSc);

SNESSolve(NS, funcJ, funcRes, xPETSc, sparams="-snes_monitor");

ChangeNumbering(NS, ug[], xPETSc, inverse = true, exchange = true);

ug[] .= NS.D;

int[int] rest = restrict(XXXXFES, XXXXFESBackup, n2o);

for[i, v : rest] ub[][v] = ug[][i];

mpiReduce(ub[], ured[], processor(0, mpiCommWorld), mpiSUM);

//Plot(ured, …)

}

Thank you advance.

Regards,

Robert