Unsteady Navier-Stokes Large Mesh

Hej,
I am doing a mesh sensitivity analysis of an unsteady three-dimensional laminar flow. I am running the code in parallel with the preconditioner “-pc_type lu”. Now I reached a mesh with ~200k nodes (~1 million elements) and it seems that I’m having an “error in optimization”:

nan != nan diff nan => Sorry error in Optimization (q) add: int2d(Th,optimize=0)(…)
remark if you add (… , optimize=2) then you remove this check (be carefull);
current line = 129 mpirank 0 / 28
Exec error : In Optimized version

Attached are the navier-stokes equations that I’m using:

varf LNS ([u, v, w, p],[hx, hy, hz, q])			
	= int3d(Th)(
 			   (1./dt)*(u*hx+v*hy+w*hz))
 	 +int3d(Th)(
 			   +(1./Re)*(dx(u)*dx(hx)+dy(u)*dy(hx)+dz(u)*dz(hx))		
 			   +(1./Re)*(dx(v)*dx(hy)+dy(v)*dy(hy)+dz(v)*dz(hy))		
 			   +(1./Re)*(dx(w)*dx(hz)+dy(w)*dy(hz)+dz(w)*dz(hz))		
 			   -p*(dx(hx)+dy(hy)+dz(hz))								
 			   -(dx(u)+dy(v)+dz(w))*q				  					   
 			   -1e-10*p*q												
 			   )
 	 +int3d(Th)(
 			   +(1./dt)*convect([up,vp,wp],-dt,up)*hx
               +(1./dt)*convect([up,vp,wp],-dt,vp)*hy
 		       +(1./dt)*convect([up,vp,wp],-dt,wp)*hz		      
			   )
 			
 			+ on(Wall, u=0., v=0., w=0.)
 			+ on(Inlet, u=0., v=0., w=0.5);

I used the informations from the theory, tutorials and other posts on this forum but I don’t manage to find an optimal combination for the preconditioners. I used several time steps dt varying from 10 to 0.01. If I add “optimize=0” to the integrals and change the pcs to: “-ksp_type fgmres -pc_type lu -pc_factor_mat_solver_type mumps” the simulation starts but I get null values. Is there something I could improve here, or in this case with a large mesh I should use fieldsplit conditioners?

Regards,
Robert

This probably comes from the convect operator you have in your varf, which is not appropriate to use in parallel unless you have a very small dt.

Thank you, Pierre. I lowered the time step as much as I could and I repetead the simulations. I attached a screenshot of the results. Do you know what could cause this post-processing error? The strange part is that the values from the color-scale are the expected ones but even if I plot the results over a line it doesn’t see any values. Another question I have would be about the Courant number. What are the tolerated limits for this value in FreeFEM++?

Regards,
Robert
image

I don’t know what I’m looking at, so I can’t help you.

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

Could you please show me the modified navier-stokes-2d-PETSc.edp example?