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=W*t; y=0; label=inlet;};
border b2(t=0., 1.){x=W; y=L*t; label=wall2;};

border b3(t=0., 1.){x=W-W

*t; y=L; label=outlet;};*

border b4(t=0., 1.){x=0; y=L-Lt; label=wall1;};

border b4(t=0., 1.){x=0; y=L-L

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)(

du1*v1
+ du2*v2

+ dt

*mu1*(dx(du1)

*dx(v1) + dy(du1)*mu1*(F/(R

*dy(v1))*

+ 2dt+ 2

*T))*(u1*(dx(dphi)

*dx(v1) + dy(dphi)*(du1

*dy(v1)) + du1*(dx(phi)*dx(v1) + dy(phi)*V/(W^2))*dy(v1)))*

- dt(6- dt

*x*(W-x)*dy(v1))*

+ dtmu2*(dx(du2)

+ dt

*dx(v2) + dy(du2)*(F/(R

*dy(v2))*

+ dtmu2+ dt

*T))*(u2*(dx(dphi)

*dx(v2) + dy(dphi)*(du2

*dy(v2)) + du2*(dx(phi)*dx(v2) + dy(phi)*V/(W^2))*dy(v2)))*

- dt(6- dt

*x*(W-x)*dy(v2))*

+ dtF*(2

+ dt

*mu1-2*mu3)

*(dx(du1)*(dx(du2)

*dx(v3) + dy(du1)*(mu2-mu3)*dy(v3))*

+ dtF+ dt

*dx(v3) + dy(du2)*T))

*dy(v3))*

+ dt(F^2/(R+ dt

*((4*mu1+2

*mu3)*(u1*(dx(dphi)

*dx(v3) + dy(dphi)*(dx(u1)

*dy(v3)) + du1*(dx(phi)*dx(v3) + dy(phi)*(dx(dphi)*dy(v3))) + (mu2+mu3)*(u2*dx(v3) + dy(dphi)*

+ dtmu1*dy(v3)) + du2*(dx(phi)*dx(v3) + dy(phi)*

+ u2v2*dy(v3))))*

)

+ int2d(dom4)(

u1v1)

+ int2d(dom4)(

u1

+ u2

+ dt

*dx(v1) + dy(u1)*mu1*(F/(R

*dy(v1))*

+ 2dt+ 2

*T))*(u1

*u1*(dx(phi)*dx(v1) + dy(phi)*V/(W^2))*dy(v1))*

- dt(6- dt

*x*(W-x)*dy(v1))*

+ dtmu2*(dx(u2)

+ dt

*dx(v2) + dy(u2)*(F/(R

*dy(v2))*

+ dtmu2+ dt

*T))*(u2

*u2*(dx(phi)*dx(v2) + dy(phi)*V/(W^2))*dy(v2))*

- dt(6- dt

*x*(W-x)*dy(v2))*

+ dtF*(2

+ dt

*mu1-2*mu3)

*(dx(u1)*(dx(u2)

*dx(v3) + dy(u1)*(mu2-mu3)*dy(v3))*

+ dtF+ dt

*dx(v3) + dy(u2)*T))

*dy(v3))*

+ dt(F^2/(R+ dt

*((4*mu1+2

*mu3)*(v1

*u1*(dx(phi)*dx(v3) + dy(phi)*V/(W^2))*dy(v3)) + (mu2+mu3)*

)

+ int1d(dom4, outlet)(dt(6*u2*(dx(phi)*dx(v3) + dy(phi)*

-uu2v2*dy(v3)))*

)

+ int2d(dom4)(

-uu1v1)

+ int2d(dom4)(

-uu1

-uu2

)

+ int1d(dom4, outlet)(dt

*x*(W-x)*u1) + dt*(6

*V/(W^2))*(v2

*x*(W-x)*u2))*

- int1d(dom4, wall2)(dt(1/(2

- int1d(dom4, wall2)(dt

*F))*F

*v1*(4*k0pb*u1

*sinh((F/(R*T))

*(-phi-Eminus-(R*T/(2

*F))*(4

*log(abs(u1)))))*

+ dtv3+ dt

*F*k0pb

*u1*sinh((F/(R

*T))*(-phi-Eminus-(R

*T/(2*F))

*log(abs(u1))))))*

+ int1d(dom4, wall1)(dt(1/(2

+ int1d(dom4, wall1)(dt

*F))*F

*v1*(4*k0pb*u1*(u2/1.5)

*sinh((F/(R*T))

*(3-phi-Eplus+(R*T/(2

*F))*F

*log(abs(u1/u2)))))*

-dt(2/F)-dt

*v2*(4*k0pb*u1*(u2/1.5)

*sinh((F/(R*T))

*(3-phi-Eplus+(R*T/(2

*F))*(4

*log(abs(u1/u2)))))*

-dtv3-dt

*F*k0pb

*u1*(u2/1.5)

*sinh((F/(R*T))

*(3-phi-Eplus+(R*T/(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,
```

}