Hello freefem community,
I am trying to solve a semilinear elliptic problem where the forcing term f(u) depends on the solution. Here I am using fixed point method implemented for nonlinear Schwarz method to solve the original problem. Unfortunately, my code is running but does not show any outputs. I guess there is an error in the fixed point method’s loop coding. Can someone comment on my problem and repair the error? With thanks.
Here is my code attempt:
int n=0; //iteration
int i=0; //number of meshes
func g= x+y; //DIRICHELET boundary condition
real x0=0, x1=9./16., y0=0, y1=1; //SUB-DOMAIN1 Dimensions
real xx0=7./16., xx1=1, yy0=0, yy1=1; //SUB-DOMAIN2 Dimensions
// Size of both meshes
real[int] h1(5) , h2(5); //size of both meshes
real Globalnorm=1.;
//SUB-DOMAIN1 borders
border AB1(t=x0, x1){x=t; y=y0; label=1;};
border BC1(t=y0, y1){x=x1; y=t; label=2;};
border CD1(t=x1, x0){x=t; y=y1; label=3;};
border DA1(t=y1, y0){x=x0; y=t; label=4;};
//SUB-DOMAIN2 borders
border AB2(t=xx0, xx1){x=t; y=yy0; label=1;};
border BC2(t=yy0, yy1){x=xx1; y=t; label=2;};
border CD2(t=xx1, xx0){x=t; y=yy1; label=3;};
border DA2(t=yy1, yy0){x=xx0; y=t; label=4;};
mesh Th1=square(2, 2, [x0+(x1-x0)*x, y0+(y1-y0)*y]); //initial Rectangular mesh of SUB-DOMAIN1
mesh Th2=square(2, 2, [xx0+(xx1-xx0)*x, yy0+(yy1-yy0)*y]); //initial Rectangular mesh of SUB-DOMAIN2
//4 MESHES OF BOTH SUB-DOMAINS (Two choices of meshing)
for (int i=0; i<=4; i++){
h1[i]=(1./8.)(1./(2.^i));
h2[i]=(1./4.)(1./(2.^i));
Th1 = adaptmesh(Th1, h1[i], IsMetric=1);
Th2 = adaptmesh(Th2, h2[i], IsMetric=1);
//Finite element spaces
fespace Vh1(Th1, P1);
Vh1 uh1, vh1, f1, R1, uold1=0, diff1;
fespace Vh2(Th2, P1);
Vh2 uh2, vh2, f2, R2, uold2=0, diff2;
// plot to see the 2 meshes.
plot(Th1, Th2, wait=true);
//solve both nonlinear subproblems
R1=-1/(1+0.25*uold1^2);
f1=uold1*V1;
problem Schw1(uh1, vh1, solver=CG, eps=1.e-6)
=int2d(Th1)(dx(uh1)*dx(vh1)+dy(uh1)dy(vh1)) //bilinear term
+int2d(Th1)(f1vh1) //right hand side(nonlinear term )
+on(1,3,4, uh1=g)+on(2, uh1=uold2); //external+internal boundaries
R2=-1/(1+0.25*uold2^2);
f2=uold2*V2;
problem Schw2(uh2, vh2, solver=CG, eps=1.e-6)
=int2d(Th2)(dx(uh2)*dx(vh2)+dy(uh2)dy(vh2)) //bilinear term
+int2d(Th2)(f2vh2) //right hand side (nonlinear term)
+on(1,2,3, uh2=g)+on(4, uh2=uh1); //external+internal boundarie
for (int n=1; n<=35; n++){
while ( Globalnorm >= 1e -6 ) { //fixed point loop for both subproblems
Schw1;
Schw2;
V1 = -1/(1+0.25*uh1^2); // actualization
V2 = -1/(1+0.25*uh2^2); // actualization
uold1=uh1;
uold2=uh2;
};
diff1 = abs(uh1-uold1);
norm1=diff1.max;
diff2 = abs(uh2-uold2);
norm2=diff2.max;
Globalnorm=max(norm1, norm2);
cout<<" --@ iteration no. --“<<n<<endl;
cout<<“Difference error on subdomain 1 at iteration”<<”(“<<n<<”) =“<<norm1<<endl;
cout<<“Difference error on subdomain 2 at iteration”<<”(“<<n<<”) ="<<norm2<<endl;
};
plot(uh1,uh2, wait=true);
};