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
+(ugdx(vg)+vgdy(vg)+wgdz(vg)) * hy
+(ugdx(wg)+vgdy(wg)+wgdz(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) (
(dudx(ug)+dvdy(ug)+dwdz(ug)) * hx
+(ugdx(du)+vgdy(du)+wgdz(du)) * hx
+(dudx(vg)+dvdy(vg)+dwdz(vg)) * hy
+(ugdx(dv)+vgdy(dv)+wgdz(dv)) * hy
+(dudx(wg)+dvdy(wg)+dwdz(wg)) * hy
+(ugdx(dw)+vgdy(dw)+wgdz(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
