Hi!
I’m new to FreeFem and parallel computing. I attended the tutorial sessions of the FreeFem Days in December, when PETSc was extensively used.
I am currently on a coupled Vlasov-Navier-Stokes system and I want numerical simulations. I began with solving Navier-Stokes. I have a running code for sequential computing but it is too slow. Therefore, I am trying to switch to parallel computing.
My code generates an error at the line where I am building the PETSc matrix. To sum up, I have Pk=[P2,P2,P2,P1] and a 3d mesh called cylindre. The line “buildMat(cylindre,getARGV(”-split",1),A,Pk,mpiCommWorld)" generates : “Error line number 241, in file macro: partition in /usr/local/ff++/mpich3/lib/ff++/4.4-3/idp/macro_ddm.idp”
I found this syntax in the stokes-fieldsplit-3d-PETSc.edp example (which runs smoothly on my computer) and I do not understand what I did differently.
Thank you all for your help!
David
Here is the full (minimal) code.
macro dimension()3//
include "macro_ddm.idp"
load "msh3"
load "tetgen"
load "PETSc"
macro def(i)[i, i#B, i#C,i#D]// // vector field definition
//Mesh
real r00 = 9.00e-1;
real L00 = 9.40;
int m=2;
int nlayer = L00*m;
border c(t=0,2*pi){x=r00*cos(t);y=r00*sin(t);};
mesh cercle = buildmesh(c(2*pi*r00*m));
int[int] rup=[0,1], rdown=[0,2], rmid=[1,3];
mesh3 cylindre = buildlayers(cercle,nlayer,
coef = 1.,
zbound = [0,L00],
labelmid = rmid,
reffaceup = rup,
reffacelow = rdown);
buildDmesh(cylindre)
// Physical data
real u0 = 0.5e2;
real rhoair = 1.18e-3;
real nu = 1e-5*1e4;
real eps = 1e-8;
real dt = 0.01;
//Finite element spaces
func Pk=[P2,P2,P2,P1];
fespace Ph(cylindre,P0);
fespace Uh(cylindre,P1);
fespace Vh(cylindre,[P2,P2,P2,P1]);
//Navier-Stokes
func u1entree = 0;
func u2entree = 0;
func u3entree = u0*(1-(x^2+y^2)/(r00^2)); // Poiseuille at entry of the cylinder
Vh [u1,u2,u3,p];
Vh [v1,v2,v3,ptest];
// Initial data
Vh [up1,up2,up3,pini];
up1=0;
up2=0;
up3 = u0*(1-(x^2+y^2)/(r00^2));
varf bilinNS([u1,u2,u3,p],[v1,v2,v3,ptest]) = int3d(cylindre)(
rhoair*(u1*v1 + u2*v2 + u3*v3)/dt
+nu*(dx(u1)*dx(v1)+dy(u1)*dy(v1)+dz(u1)*dz(v1)
+dx(u2)*dx(v2)+dy(u2)*dy(v2)+dz(u2)*dz(v2)
+dx(u3)*dx(v3)+dy(u3)*dy(v3)+dz(u3)*dz(v3)
)
- p*(dx(v1)+dy(v2)+dz(v3))
- ptest*(dx(u1)+dy(u2)+dz(u3))
+ p*ptest*eps
)+
int3d(cylindre)(rhoair*(
-convect([up1,up2,up3],-dt,up1)*v1/dt
-convect([up1,up2,up3],-dt,up2)*v2/dt
-convect([up1,up2,up3],-dt,up3)*v3/dt
))+
on(1, u1=u1entree, u2=u2entree, u3=u3entree)+on(3, u1=0, u2=0, u3=0);
matrix<real> Loc = bilinNS(Vh,Vh);
real[int] rhs = bilinNS(0,Vh);
Mat A;
**buildMat(cylindre,getARGV("-split",1), A, Pk, mpiCommWorld)**
A = Loc;
//From here on I am copying the example stokes-fieldsplit-3d-PETSc.edp, I do not really know what is happening
Vh<real> def(sol) = [1.0,1.0,1.0,2.0];
string[int] names(2);
names[0] = "velocity";
names[1] = "pressure";
set(A, sparams = "-ksp_type fgmres -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type lower -pc_fieldsplit_detect_saddle_point -fieldsplit_velocity_sub_pc_type lu " + " -fieldsplit_pressure_sub_pc_type lu -fieldsplit_velocity_sub_pc_factor_mat_solver_type mumps " + " -fieldsplit_pressure_sub_pc_factor_mat_solver_type mumps -ksp_monitor -ksp_view " + " -fieldsplit_velocity_ksp_type gmres -fieldsplit_velocity_ksp_max_it 5 -fieldsplit_pressure_ksp_type gmres -fieldsplit_pressure_ksp_max_it 5 -ksp_rtol 1e-6", fields = sol[], names = names)
sol[] = 0.0;
sol[] = A^-1 * rhs;
plotD(cylindre,sol,cmm="sol distribuee");