solver.edp

// File that solve the Navier Stokes stationary in time using a fixed-point method for the
//different domain: "LAA_CF.msh", "LAA_CW.msh", "LAA_WS.msh", "LAA_CS.msh"
// In the end, indicators for the probability of formation of a clot are defined, the aim is to classify the domains for the most probable to the less probable

// 1. geometry
load "gmsh"
mesh Th = gmshload("LAA_CF.msh"); //load of the mesh

// 2. finite element spaces
fespace Vh(Th, [P2, P2, P1]); //FE space for velocity ux,uy and pressure p
fespace Ah(Th,P2); //Fe space for the total velocity u=sqrt(ux^2+uy^2)

// 3. set data of the problem
real rho=0.00106; //density [kg/cm^3]
real mu=3.5e-5; //dynamic viscosity [kg/(cm*s)]
real nu=mu/rho; //kinematic viscosity [cm^2/s]
real u=35; //characteristic velocity (the one at inlet) [cm/s]
real L=1; //characteristic length (the one at inlet) [cm]
real Re = u*L/nu; //Reynolds number

// 4. define field of the FE spaces
//  [ux0h, uy0h, p0h] is the solution at the previous iteration
// [uxh, uyh, ph] is the solution at the actual iteration
//  [vxh, vyh, qh] are the test functions

Vh [ux0h, uy0h, p0h],
   [uxh, uyh, ph],
   [vxh, vyh, qh];

// uh is the total velocity
Ah uh;

// 5. variational forms
//macro functions
macro Grad(f) [dx(f), dy(f)] //
macro Grad2(fx,fy) [dx(fx), dy(fx), dx(fy), dy(fy)] //
macro Div(fx,fy) (dx(fx) + dy(fy)) //
macro UGrad(bx,by, f) [bx,by]' * Grad(f) //
macro UGrad2(bx,by, fx,fy) [UGrad(bx,by,fx), UGrad(bx,by,fy)] //

//variational formulation of the dimensionless form of the problem, using the fixed-point method
problem fixedpoint([uxh, uyh, ph], [vxh, vyh, qh]) =
    int2d(Th)(UGrad2(ux0h, uy0h, uxh, uyh)' * [vxh, vyh]
              + (Grad2(uxh, uyh)' * Grad2(vxh, vyh)) / Re
              - ph * Div(vxh, vyh)
              + Div(uxh, uyh) * qh)
    + on(2, uxh=0.0, uyh=-1.0)
    + on(1, uxh=0.0, uyh=0.0);

//variational fomrulation of the residual of the problem
varf residual([ux,uy,p],[vx,vy,q]) =
    int2d(Th)(UGrad2(uxh, uyh, uxh, uyh)' * [vx, vy]
              + (Grad2(uxh, uyh)' * Grad2(vx, vy)) / Re
              - ph * Div(vx, vy)
              + Div(uxh, uyh) * q)
                + on(2, ux=0, uy=0)
    + on(1, ux=0, uy=0);

// 6. functions
func int step() {
    // input: [ux0h, uy0h, p0h] -> old solution
    // output: [uxh, uyh, ph] -> new solution
    fixedpoint;
    return 0;
}

// calculate the L2-norm of the residual
real[int] res(Vh.ndof);
func real error() {
    res = residual(0, Vh);
    return sqrt(res' * res);
}

// 7. nonlinear solver
int nbiter = 80; //max iterations
real eps = 1e-5; //threshold for the error

[ux0h, uy0h, p0h] = [0.0, 0.0, 0.0]; //initial solution
fixedpoint; //at least one iteration is done
[ux0h, uy0h, p0h] = [uxh, uyh, ph]; //update the solution

int iter = 0;
for (iter = 0; iter < nbiter; ++iter) {
    int errc = step();
    if (errc != 0) { //check if the variational form of the problem is solved
        cout << "ERROR: iteration " << iter << " failed!" << endl;
        break;
    }

    real err=error(); //calculate the L2-norm of the residual
    if (err < eps) //if the residual is less than the threshold
        break; //the fixed-point method has reached convergence
    [ux0h, uy0h, p0h] = [uxh, uyh, ph]; //update the solution
}

cout << "# iter: " << iter << endl; //iteration for which has reached oncergence

// 8. plot solution
plot(ph, nbiso=100,fill=1, wait=1, value=1, cmm="Pressure");
plot(uxh, nbiso=100, fill=1, wait=1, value=1,cmm="Velocity X");
plot(uyh, nbiso=100, fill=1, wait=1,value=1,cmm="Velocity Y");
plot([uxh, uyh], fill=1,value=1, cmm="Velocity");

// 9. Indicators Thrombous
uh=sqrt(uxh^2+uyh^2); //define the total velocity

//define the indicators
//ind1: area of still blood / area total
//ind2: mean of the blood / area total
real ind1= int2d(Th)((abs(uh)<0.05556))/Th.area;
real ind2=int2d(Th)(uh)/Th.area;

cout<<"The first indicator is "<<ind1<<endl;
cout<<"The second indicator is "<<ind2<<endl;