1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
|
// run with MPI: ff-mpirun -np 4 script.edp
// NBPROC 4
load "PETSc" // PETSc plugin
macro partitioner()metis// EOM // metis, scotch, or parmetis
macro dimension()3// EOM // 2D or 3D
include "macro_ddm.idp" // additional DDM functions
macro def(i)[i, i#B, i#C, i#D]// EOM// vector field definition
macro init(i)[i, i, i, i]// EOM // vector field initialization
macro grad(u)[dx(u), dy(u), dz(u)]//// two-dimensional gradient
real Sqrt = sqrt(2.);
macro div(u)(dx(u) + dy(u#B) + dz(u#C))// EOM
func Pk = [P2, P2, P2, P1]; // finite element space
int s = getARGV("-split", 1); // refinement factor
if(verbosity > 0 && mpirank == 0) {
cout << " --- " << mpirank << "/" << mpisize;
cout << " - stokes-3d-PETSc.edp - input parameters: refinement factor = " << s << endl;
}
meshN Th = buildlayers(square(1, 1), 1);
fespace Wh(Th, Pk); // local finite element space
int[int] arrayIntersection; // ranks of neighboring subdomains
int[int][int] restrictionIntersection(0); // local-to-neighbors renumbering
real[int] D; // partition of unity
{
mesh ThGlobal2d = square(getARGV("-global", 12), getARGV("-global", 12), [x, y]); // global mesh
ThGlobal2d = trunc(ThGlobal2d, (x <= 0.5) || (y <= 0.5), label = 5);
ThGlobal2d = trunc(ThGlobal2d, (y >= 0.25) || (x >= 0.25), label = 5);
mesh Th2d = movemesh(ThGlobal2d, [-x, y]);
ThGlobal2d = ThGlobal2d + Th2d;
meshN ThBorder, ThGlobal = buildlayers(ThGlobal2d, getARGV("-global", 12) / 2, zbound = [0, 0.4]);
build(Th, ThBorder, ThGlobal, 10, s, 1, D, arrayIntersection, restrictionIntersection, Wh, Pk, mpiCommWorld, false)
}
varf vPb([u, uB, uC, p], [v, vB, vC, q]) = intN(Th)(grad(u)' * grad(v) + grad(uB)' * grad(vB) + grad(uC)' * grad(vC) - div(u) * q - div(v) * p + 1e-10 * p * q) + on(0, 1, 3, 5, u = 0, uB = 0, uC = 0) + on(2, u = 1000*y*(0.5-y)*z*(0.4-z), uB = 0, uC = 0);
matrix<real> A = vPb(Wh, Wh, tgv = -1);
real[int] rhs = vPb(0, Wh, tgv = -1);
dmatrix Mat(A, arrayIntersection, restrictionIntersection, D, bs = 1);
set(Mat, sparams = "-pc_type lu -pc_factor_mat_solver_package mumps");
Wh<real> def(b);
b[] = Mat^-1 * rhs;
plotMPI(Th, b[], "Global solution PETSc", Pk, def, real, 3, 1);
|