Hello,
I am trying to write Freefem++ code for the unsteady Navier Stokes equations in 3D using Newton’s method. I found a similar example in the FreeFem hpddm examples. I changed it into 3d version.
The problem is when I set the MaxLayer = 5, everything goes fine and the result looks reasonable. But when I change the MaxLayer bigger than 20, the solver didn’t converge. Something wrong with the parameter in SNESSolve?
The code is as follows:
///////////////////////////////////////////////////////////////////////////////////
load “PETSc”
macro dimension()3//
include “macro_ddm.idp”
verbosity=0;
real dt = 0.001;
real T = 200;
int saveEach = 0.005/dt;
int[int] orderOut = [1, 1, 1, 1, 1, 1];
string outputFileName = “output/sub-3d.vtu”;
/*
int Wall = 1;
int Inlet = 3;
int Outlet = 2;
int InnerWall = 4;
*/
real uin = 0.3;
real nu = 0.05;
macro velu(u) [u#x, u#y,u#z]//
macro div(u) (dx(u#x) + dy(u#y) + dz(u#z))//
macro grad(u) [dx(u), dy(u), dz(u)]//
macro Grad(u) [dx(u#x),dy(u#x),dz(u#x),dx(u#y),dy(u#y),dz(u#y),dx(u#z),dy(u#z),dz(u#z)]// EOM
macro UgradV(u,v) [ [u#x,u#y,u#z]‘[dx(v#x),dy(v#x),dz(v#x)] ,
[u#x,u#y,u#z]'[dx(v#y),dy(v#y),dz(v#y)] ,
[u#x,u#y,u#z]’*[dx(v#z),dy(v#z),dz(v#z)] ]//
mesh3 Mesh;
{
real r=0.0475;
real pitch=0.063;
real gap=2*(0.063-0.0475);
int nn=10;
if (mpirank == 0)
{
border C01(t=0, pi/2){x=-pitch+rcos(t); y=-pitch+rsin(t); label=1;}
border C02(t=0, gap){x=-pitch; y=-pitch+r+t; label=2;}
border C03(t=0, pi/2){x=-pitch+rsin(t); y=pitch-rcos(t); label=3;}
border C04(t=0, gap){x=-pitch+r+t; y=pitch; label=4;}
border C05(t=0, pi/2){x=pitch-rcos(t); y=pitch-rsin(t); label=5;}
border C06(t=0, gap){x=pitch; y=pitch-r-t; label=6;}
border C07(t=0, pi/2){x=pitch-rsin(t); y=-pitch+rcos(t); label=7;}
border C08(t=0, gap){x=pitch-r-t; y=-pitch; label=8;}
mesh Th = buildmesh(C01(-20)+C02(-10)+C03(-20)+C04(-10)+C05(-20)+C06(-10)+C07(-20)+C08(-10)
);
int[int] rup=[0,10], rdown=[0,9], rmid=[1,1,2,2,3,3,4,4];
int[int] r1 = [0, 11];
func zmin = 0.;
func zmax = 0.5;
int MaxLayer = 5;
Mesh = buildlayers(Th, MaxLayer, zbound=[zmin, zmax], region=r1,labelmid=rmid,labelup = rup,labeldown = rdown);
}
broadcast(processor(0), Mesh);
}
if (mpirank == 0)
cout << "Number of Elements: " + Mesh.nt << endl;
func PkVector = [P2, P2, P2, P1];
buildDmesh(Mesh);
Mat J;
{
macro def(i)[i, i#B, i#C, i#D]//
macro init(i)[i, i, i, i]//
createMat(Mesh, J, PkVector)
}
fespace SpaceVector(Mesh, PkVector);
SpaceVector [ucx, ucy, ucz, pc] = [0,0,uin, 0];
SpaceVector [uckx, ucky, uckz, pck]= [0,0,uin, 0];
SpaceVector [ucoldx, ucoldy, ucoldz, pcold]= [0,0,uin, 0];
if (mpirank == 0)
cout << "Finite Element DOF (in each partition): " + SpaceVector.ndof << endl;
varf vRes([ux, uy, uz, p], [vx, vy, vz, q]) = int3d(Mesh)(
(1/dt*velu(uck)’ * velu(v))
- (1/dt*velu(ucold)’ * velu(v)) +
nu * Grad(uc)'*Grad(v)- UgradV(uc, uc)’ * velu(v)
- pc * div(v) - div(uc) * q
//+ 0.000000001pcq
)
- on(1, 3,5,7, ux = ucx-0, uy = ucy-0, uz = ucz-0)
- on(2,6, ux = ucx-0)
+on(4,8,uy=ucy-0) - on(9, ux = ucx-0, uy = ucy-0, uz = ucz-uin);
varf vJ([ux, uy, uz, p], [vx, vy, vz, q]) = int3d(Mesh)(
(1/dtvelu(u)’ * velu(v)) +
(UgradV(uc, u) + UgradV(u, uc))’ * velu(v) +
nu * Grad(u)’ Grad(v)
- p * div(v) - div(u) * q
//+0.000000001pq
)
+ on(1, 3,5,7, ux = ucx-0, uy = ucy-0, uz = ucz-0)
+ on(2,6, ux = ucx-0)
+on(4,8,uy=ucy-0)
+ on(9, ux = ucx-0, uy = ucy-0, uz = ucz-uin);
set(J, sparams = “-pc_type lu”);
func real[int] funcRes(real[int]& inPETSc) {
ChangeNumbering(J, ucx, inPETSc, inverse = true, exchange = true);
[uckx, ucky, uckz, pck]=[ucx, ucy, ucz, pc];
real[int] out(SpaceVector.ndof);
out = vRes(0, SpaceVector, tgv = -1);
ChangeNumbering(J, out, inPETSc);
return inPETSc;
}
func int funcJ(real[int]& inPETSc) {
ChangeNumbering(J, ucx, inPETSc, inverse = true, exchange = true);
J = vJ(SpaceVector, SpaceVector, tgv = -1);
return 0;
}
for(int i = 0; i < T/dt; i++)
{
if (mpirank == 0)
cout << "time step: " << i << endl;
[uckx, ucky, uckz, pck] = [ucx, ucy, ucz, pc];
real[int] xPETSc;
ChangeNumbering(J, ucx, xPETSc);
SNESSolve(J, funcJ, funcRes, xPETSc, sparams = “-snes_monitor”);
ChangeNumbering(J, ucx, xPETSc, inverse = true, exchange = true);
[ucoldx, ucoldy, ucoldz, pcold]=[ucx, ucy, ucz, pc];
if (i % saveEach == 0)
savevtk(outputFileName, Mesh, ucx, ucy, ucz,[ucx,ucy,ucz], pc, dataname=“u v w velocity p”, order=orderOut,
append = i ? true : false);
}