Unsteady Navier-Stokes Large Mesh

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)(
 			+ 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?


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++?


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);
XXXXFESBackup [ub, vb, wb, pb];
XXXXFESBackup [ured, vred, wred, pred];
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
dx(wg)+vgdy(wg)+wgdz(wg)) * hz
+(1/Re)*(dx(ug)*dx(hx)+dy(ug)dy(hx)+dz(ug) * dz(hx))
(dx(vg)*dx(hy)+dy(vg)dy(hy)+dz(vg) * dz(hy))
(dx(wg)*dx(hz)+dy(wg)dy(hz)+dz(wg) * dz(hz))
- pg
+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
dx(du)+vgdy(du)+wgdz(du)) * hx
+(dudx(vg)+dvdy(vg)+dwdz(vg)) * hy
dx(dv)+vgdy(dv)+wgdz(dv)) * hy
+(dudx(wg)+dvdy(wg)+dwdz(wg)) * hy
dx(dw)+vgdy(dw)+wgdz(dw)) * hy
+(1/Re)*(dx(du)*dx(hx)+dy(du)dy(hx)+dz(du) * dz(hx))
(dx(dv)*dx(hy)+dy(dv)dy(hy)+dz(dv) * dz(hy))
(dx(dw)dx(hz)+dy(dw)dy(hz)+dz(dw) * 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.

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