# Unsteady Navier-Stokes Large Mesh **URL:** \ **Category:** General Discussion\ **Created:** [June 7, 2021, 4:53pm UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018 "2021-06-07T16:53:09Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![Robert](https://avatars.discourse-cdn.com/v4/letter/r/e95f7d/32.png) [@Robert](https://community.freefem.org/u/Robert)\ **Post date:** [June 7, 2021, 4:53pm UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/1 "2021-06-07T16:53:09Z")
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: ```auto 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 ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [June 7, 2021, 5:45pm UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/2 "2021-06-07T17:45:19Z")
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`. ---
**Author:** ![Robert](https://avatars.discourse-cdn.com/v4/letter/r/e95f7d/32.png) [@Robert](https://community.freefem.org/u/Robert)\ **Post date:** [June 11, 2021, 2:42pm UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/3 "2021-06-11T14:42:02Z")
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](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/af956d15d36fc5e6cf120fd5cdbde0bde22a8575.png) ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [June 11, 2021, 2:46pm UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/4 "2021-06-11T14:46:33Z")
I don’t know what I’m looking at, so I can’t help you. ---
**Author:** ![Robert](https://avatars.discourse-cdn.com/v4/letter/r/e95f7d/32.png) [@Robert](https://community.freefem.org/u/Robert)\ **Post date:** [June 15, 2021, 7:06pm UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/5 "2021-06-15T19:06:24Z")
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](https://github.com/prj-/moulin2019al) and mostly [Ex.2](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/navier-stokes-2d-PETSc.edp) 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 > +(ug_dx(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 > +(ug_dx(du)+vg_dy(du)+wg_dz(du)) \* hx > +(du_dx(vg)+dv_dy(vg)+dw_dz(vg)) \* hy > +(ug_dx(dv)+vg_dy(dv)+wg_dz(dv)) \* hy > +(du_dx(wg)+dv_dy(wg)+dw_dz(wg)) \* hy > +(ug_dx(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 ---
**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\ **Post date:** [June 16, 2021, 5:44am UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/6 "2021-06-16T05:44:04Z")
Could you please show me the modified `navier-stokes-2d-PETSc.edp` example? ---
**Author:** ![Qianrui](https://avatars.discourse-cdn.com/v4/letter/q/9fc348/32.png) [@Qianrui](https://community.freefem.org/u/Qianrui)\ **Post date:** [November 28, 2021, 2:19am UTC](https://community.freefem.org/t/unsteady-navier-stokes-large-mesh/1018/7 "2021-11-28T02:19:39Z")
Hello, I also met the same problem for solving klein-Gordon, is there a solution?