Mesh becoming overly fine during adaptation

jna1g18 · July 18, 2023, 11:55am

Hi all, I’m quite new to freefem so please try to bear with me. I’ve been trying to implement a mesh adapation into my code for a nonlinear system of three equations in order to make it more efficient. My code uses a Newton scheme to deal with the nonlinearity, which seems to work fine, but I cannot get the mesh adaptation to behave properly; my problem is that with each timestep, the adapted mesh becomes increasingly fine to point where it makes my code less efficient than it was without any mesh adaptation and I see no reason why it should be doing this. The mesh adaptation does not seem to be following the contours of u1, u2 and phi. See below the initial mesh:

At timestep t=0.1, the adapted mesh looks good and seems to match u1, u2 and phi:

But then it becomes unnecessarily fine for t=0.2 and t=0.3 etc.:

Why is this happening? Why does the mesh become so fine? The mesh does not seem to be adapting to the functions that I have input into the adaptmesh() function. Below I have attached my code for reference (apologies, it is a little messy). Please let me know if my question is not clear.

//Time step
real dt = 0.1;

//Model parameters
real mu1 = 0.001;
real mu2 = 0.001;
real mu3 = 0.001;
real V = 0.2;
real D = 0.01;
real Vr = 0.0015;
real F = 96485.3;
real T = 298;
real R = 8.31;
real k0pb = 0.00000021;
real k0pb02 = 0.00000025;
real Eminus = 0.13;
real Eplus = 1.56;

//Initial adaptation error
real Eadapt = 0.1;

//Error threshold and initial error for Newton loop
real eps = 1e-6;
real err = 1;

//Construction of initial mesh
real L = 0.02;
real W = 0.012;
int meshSize2 = 40;
int meshsize3 = meshSize2*(L/W);

int wall1 = 1;
int wall2 = 2;
int inlet = 3;
int outlet = 4;

mesh dom4;

border b1(t=0., 1.){x=Wt; y=0; label=inlet;};
border b2(t=0., 1.){x=W; y=Lt; label=wall2;};
border b3(t=0., 1.){x=W-Wt; y=L; label=outlet;};
border b4(t=0., 1.){x=0; y=L-Lt; label=wall1;};

dom4 = buildmesh(b1(meshSize2) + b2(meshsize3) + b3(meshSize2) + b4(meshsize3));
plot(dom4, wait=true);

//Defining function spaces, variables and test functions (note, test functions are denoted as v followed by a number e.g. v1)
fespace Vh(dom4, P2);
Vh du1, v1, u1, uu1, uin1, uin01, uout01;

fespace Vh2(dom4, P2);
Vh2 du2, v2, u2, uu2, uin2, uin02, uout02;

fespace Vh3(dom4, P2);
Vh3 dphi, v3, phi;

//Setting up problem (the unknowns are the corrections du1, du2 and dphi to the variables u1, u2 and phi)
problem dHeat (du1, du2, dphi, v1, v2, v3)
=int2d(dom4)(
du1v1
+ du2v2
+ dtmu1(dx(du1)dx(v1) + dy(du1)dy(v1))
+ 2dtmu1*(F/(RT))(u1*(dx(dphi)dx(v1) + dy(dphi)dy(v1)) + du1(dx(phi)dx(v1) + dy(phi)dy(v1)))
- dt(6V/(W^2))x(W-x)(du1dy(v1))
+ dtmu2*(dx(du2)dx(v2) + dy(du2)dy(v2))
+ dtmu2(F/(RT))(u2*(dx(dphi)dx(v2) + dy(dphi)dy(v2)) + du2(dx(phi)dx(v2) + dy(phi)dy(v2)))
- dt(6V/(W^2))x(W-x)(du2dy(v2))
+ dtF*(2mu1-2mu3)(dx(du1)dx(v3) + dy(du1)dy(v3))
+ dtF(mu2-mu3)(dx(du2)dx(v3) + dy(du2)dy(v3))
+ dt(F^2/(RT))((4mu1+2mu3)(u1*(dx(dphi)dx(v3) + dy(dphi)dy(v3)) + du1(dx(phi)dx(v3) + dy(phi)dy(v3))) + (mu2+mu3)(u2(dx(dphi)dx(v3) + dy(dphi)dy(v3)) + du2(dx(phi)dx(v3) + dy(phi)dy(v3))))
)
+ int2d(dom4)(
u1v1
+ u2v2
+ dtmu1(dx(u1)dx(v1) + dy(u1)dy(v1))
+ 2dtmu1*(F/(RT))u1(dx(phi)dx(v1) + dy(phi)dy(v1))
- dt(6V/(W^2))x(W-x)(u1dy(v1))
+ dtmu2*(dx(u2)dx(v2) + dy(u2)dy(v2))
+ dtmu2(F/(RT))u2(dx(phi)dx(v2) + dy(phi)dy(v2))
- dt(6V/(W^2))x(W-x)(u2dy(v2))
+ dtF*(2mu1-2mu3)(dx(u1)dx(v3) + dy(u1)dy(v3))
+ dtF(mu2-mu3)(dx(u2)dx(v3) + dy(u2)dy(v3))
+ dt(F^2/(RT))((4mu1+2mu3)u1(dx(phi)dx(v3) + dy(phi)dy(v3)) + (mu2+mu3)u2(dx(phi)dx(v3) + dy(phi)dy(v3)))
)
+ int2d(dom4)(
-uu1v1
-uu2v2
)
+ int1d(dom4, outlet)(dt(6V/(W^2))x(W-x)(v1u1) + dt(6V/(W^2))x(W-x)(v2u2))
- int1d(dom4, wall2)(dt(1/(2F))v1(4Fk0pbu1sinh((F/(RT))(-phi-Eminus-(RT/(2F))log(abs(u1)))))
+ dtv3(4Fk0pbu1sinh((F/(RT))(-phi-Eminus-(RT/(2F))log(abs(u1))))))
+ int1d(dom4, wall1)(dt(1/(2F))v1(4Fk0pbu1*(u2/1.5)sinh((F/(RT))(3-phi-Eplus+(RT/(2F))log(abs(u1/u2)))))
-dt(2/F)v2(4Fk0pbu1*(u2/1.5)sinh((F/(RT))(3-phi-Eplus+(RT/(2F))log(abs(u1/u2)))))
-dtv3(4Fk0pbu1(u2/1.5)sinh((F/(RT))(3-phi-Eplus+(RT/(2*F))*log(abs(u1/u2))))))

+ int1d(dom4, outlet)(dt*(6*V/(W^2))*x*(W-x)*(v1*du1) + dt*(6*V/(W^2))*x*(W-x)*(v2*du2))
- int1d(dom4, wall2)(dt*(1/(2*F))*v1*(du1*4*F*k0pb*(sinh((F/(R*T))*(-phi-Eminus-(R*T/(2*F))*log(abs(u1)))) - 0.5*cosh((F/(R*T))*(-phi-Eminus-(R*T/(2*F))*log(abs(u1))))) + dphi*(-4*F^2*k0pb/(R*T))*u1*cosh((F/(R*T))*(-phi-Eminus-(R*T/(2*F))*log(abs(u1)))))
+ dt*v3*(du1*4*F*k0pb*(sinh((F/(R*T))*(-phi-Eminus-(R*T/(2*F))*log(abs(u1)))) - 0.5*cosh((F/(R*T))*(-phi-Eminus-(R*T/(2*F))*log(abs(u1))))) + dphi*(-4*F^2*k0pb/(R*T))*u1*cosh((F/(R*T))*(-phi-Eminus-(R*T/(2*F))*log(abs(u1))))))
+ int1d(dom4, wall1)(dt*(1/(2*F))*v1*(du1*4*F*k0pb02*(u2/1.5)*(sinh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))) + 0.5*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))))
+ du2*4*F*k0pb02*(u1/1.5)*(sinh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))) - 0.5*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2))))) + dphi*(-4*F^2*k0pb02/(R*T))*u1*(u2/1.5)*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))))
-dt*(2/F)*v2*(du1*4*F*k0pb02*(u2/1.5)*(sinh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))) + 0.5*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))))
+ du2*4*F*k0pb02*(u1/1.5)*(sinh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))) - 0.5*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2))))) + dphi*(-4*F^2*k0pb02/(R*T))*u1*(u2/1.5)*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))))
-dt*v3*(du1*4*F*k0pb02*(u2/1.5)*(sinh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))) + 0.5*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))))
+ du2*4*F*k0pb02*(u1/1.5)*(sinh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2)))) - 0.5*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2))))) + dphi*(-4*F^2*k0pb02/(R*T))*u1*(u2/1.5)*cosh((F/(R*T))*(3-phi-Eplus+(R*T/(2*F))*log(abs(u1/u2))))))
+ on(inlet, du1=uin1-u1)
+ on(inlet, du2=uin2-u2)
;

//initialisation

real t=0;

uu1=0.7; //old value of u1 (value at the previous time step)
uin01=0.7; //initial inlet value of u1
uout01=0.7; //intial outlet value of u1
u1=0.7; //Initial guess for Newton iteration (note, this is just the value at t=0)

uu2=1.5; //old value of u2 (value at the previous time step)
uin02=1.5; //initial inlet value of u2
uout02=1.5; //intial outlet value of u2
u2=1.5; //Initial guess for Newton iteration (note, this is just the value at t=0)

phi=1; //starting guess for phi in Newton iteration

for (int m=0; m<=5/dt; m++){ //Time loop

t = t+dt;


uin1 =   uin01 + dt*(V*D*W/Vr)*(uout01-uin01);  //forward difference scheme applied to equation for inlet value of u1
uin01 = uin1;   

uin2 =   uin02 + dt*(V*D*W/Vr)*(uout02-uin02);  //forward difference scheme applied to equation for inlet value of u2
uin02 = uin2;

for(int G=0; G<4; G++){  //Mesh adaptation loop
    for (int i=0; i<=50; i++){  //Newton iteration loop
        dHeat;  
        u1[]+=du1[];    //add the corrections to the current approximations for the unknowns   i.e. u1=u1+du1 etc.
        u2[]+=du2[];
        phi[]+=dphi[];
        real Lu1= u1[].linfty, Lu2=u2[].linfty, Lphi=phi[].linfty;
        err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dphi[].linfty/Lphi;
        cout<<i<<"err="<<err<<" eps="<<eps<<endl;

        if(err<eps) break;  //If the error is less than the threshold value, end the Newton loop
        if(i>3 && err>10) break;  //If the error begins to diverge, end the Newton loop
}
dom4 = adaptmesh(dom4, u1, u2, phi, err=Eadapt, nbvx=40000, hmin=0.0001);   //Adapt the mesh to the newly obtained solutions
u1=u1;  
u2=u2;      //Interpolate the newly found solutions upon the new version of the mesh 
phi=phi;
Eadapt=Eadapt/2; 
}


uu1 = u1;   //update solution at previous timestep
uout01 = int1d(dom4, outlet)(u1)/W;  //Update the old outlet value

uu2 = u2;
uout02 = int1d(dom4, outlet)(u2)/W;

//Plot the results for this current timestep
plot(phi, value=true, wait=true); 
    

plot(u2, value=true, wait=true); //fill=true,
    

plot(u1, value=true, wait=true); //fill=true,

}

aszaboa · July 18, 2023, 3:00pm

I played around with 2D mesh adaptation a few years ago. For me, it turned out not to be worth it.
Since your problem is not that big (mesh limit max 40000) you might just do better with creating a sufficiently fine mesh by hand.

I think the adaptation is based of the change (derivative) of the function you supply. As your solution seems to be sufficiently smooth, the program just refines everywhere.

You can play around with the error parameter, max mesh size, and specify different functions to adapt to (e.g., solution^2 or something).

marchywka · July 18, 2023, 3:22pm

The documentation on adaptmesh should help but probably you want
it to be fine where things vary quickly. I guess if you have an expoential,
it will vary more quickly where it is large but you can play with that
and also the number of veritcies allowed.

jna1g18 · July 18, 2023, 6:45pm

Hi, thanks for taking the time to reply. I’ll take onboard what you’ve told me! The code I posted is a solution to a toy model; in the real model, the diffusion coefficients (parameters mu1, mu2, mu3) are significantly smaller and, consequently, a boundary layer forms at both wall1 and wall2, over which the variables u1 and u2 vary quite rapidly. Do you think that a mesh adaptation would be good to employ in this case?

jna1g18 · July 18, 2023, 6:48pm

Hi, thanks for the response! I’ve had quite a close look at the adaptmesh section and I tried a few different things, such as lowering nbvx and increasing hmin; I had some success, however it merely shifted the issues onto a later timestep unfortunately.

marchywka · July 18, 2023, 8:22pm

I’m working on a phase boundary problem and used something like this,
“c” is the phase dof and the others are things like charge distribution.

` Thnew=adaptmesh(Thnew, 1.0Grad22(c)+10Grad22(pmn)+10Grad22(cln)+1e5Grad22(v),hmin=0.2*.001szy ,hmax=.03szy,abserror=1, nbvx=1e5);

or closer,

but the reason I’m posting is to show this lol,

which gives me more control over how the mesh is organized. In this case,
I’m defining a local etch rate based on concentrations solved by FF iteration.
that move this phase boundary down from the liquid above. You could probably
do something similar by manipulating isoline implemented by freefem
but there are a lot of complexities when actually moving a surface, see
for example the freefem movemesh situation. I just trash the existing mesh,
create an array of points at multiples of the local etch amount and use freefem
triangulate on the points. I guess regular meshes could be pathological
and with enough points the result should not be sensitive to the mesh but
there are issues with adaptation like conserving amount of stuff etc.
This also lets me play with interpolation issues. fwiw.

`

jna1g18 · July 19, 2023, 2:30pm

This looks interesting, I’ll certainly have a close look at this!

cmd · July 20, 2023, 11:09am

@aszaboa I would encourage you to try again with adaptation. It is very powerful, and I personally haven’t encountered problems on “smooth” solutions (quite the opposite, in fact). If the mesh becomes too coarse for your liking, you can adjust the err or hmax parameters.

cmd · July 20, 2023, 11:11am

Can you re-post the file that is giving you issues as an attachment? The code did not paste correctly.

jna1g18 · July 20, 2023, 11:28am

Hi Chris, hopefully this works. Let me know if it doesn’t or if there is a better way to upload my code.

ToPostOnForum.edp (7.7 KB)

cmd · July 20, 2023, 11:44am

It is just a silly mistake: you are overwriting the value of Eadapt. Here’s a fix:
forumfix.edp (7.5 KB)

jna1g18 · July 20, 2023, 3:12pm

Ah wow, thank you! I’m not sure how I missed this XD

Topic		Replies	Views
Mesh adaptation General Discussion	2	212	March 10, 2023
Adaptive FEM about adaptmesh General Discussion	0	208	March 23, 2023
Implementing a non-constant stabilized parameter General Discussion	3	325	July 27, 2022
How to calculate adaptive meshing time? General Discussion	0	171	August 30, 2022
Integrals of functions defined on different meshes General Discussion	0	217	July 8, 2022

Mesh becoming overly fine during adaptation

Related Topics