Implementing FreeFem on a problem including equation in the non weak form format

General Discussion
1

I have a system of partial differential equations (PDEs) that I’ve discretized in time using a specific scheme. Now, I want to discretize it in space using finite elements. However, I’m facing an issue: my last equation doesn’t conform to the weak form. The weak formulation of my problem is as follows:

Let’s denote our domain as Ω, and W and P as test functions from appropriate function spaces. Additionally, Μ^{k+1}, U^{k+1}, and R^{k+1} are our unknowns at time step k+1, and we have the values of all variables up to time step k.

Here’s the weak formulation with all three equations:

<U^{k+1}, W> = -<∇Μ^{k+1}, W>
<Μ^{k+1}, P> = <R^{k+1}, P> + <∇U^{k+1}, ∇P>
R^{k+1} = ∫_Ω A(U^k) * U^{k+1} dΩ

However, the last equation doesn’t adhere to the weak form. My goal is to incorporate it into the second equation in a form that satisfies the weak formulation.

I’m looking for guidance on how to implement this in FreeFEM++. Specifically, I need to address the issue where the last equation doesn’t fit the weak form, and I want to integrate it into the second equation while ensuring the overall problem remains consistent with the weak formulation.

Any help or insights would be greatly appreciated. Thank you

I attempted to integrate the third equation into the second equation, but FreeFEM is encountering an error and is unable to solve the problem.

2

Am I correct in understanding that R is a scalar? If so, it seems natural to include the last equation as a Lagrange multiplier with an additional scalar-valued test function.

1 Like
3

Thank you for your response. R is a scalar function of t and is independent of both x and y. I simplify my system in fact. Could you provide me an example of how I can do it in Freefem++?

4

I don’t know if there is an easier way but if you are going to go rogut
anyway generally the examples in the documentation IIRC use
matrix and RHS extraction from an ff “varf” and then augmenting the
matrix with Lagrange constraint. IIRC that is in the documentation and
comes up on some keyword search.

1 Like
5

Thank you for your response. I’ll certainly give the method you explained a try. However, in the meantime, I would greatly appreciate it if you could let me know if there is a simpler approach available. If there’s an easier method mentioned in the documentation or if I can find it through a keyword search, please do inform me. Your assistance would be highly valued. Thank you again!

6

There is an example that uses a Lagrange multiplier to enforce a constraint here.

1 Like
7

Great. Thanks a lot.

8

I still have problems solving my system in freefem++. I want to make sure that it is possible to do it in freefem++ and if not I will try to use another software like Fenics or MatLab. This is my big system in weak form except for the last equation :

image
image985×649 150 KB

. considering that \bar{r}^{k+1/2}=0.5(r^k+r^{k+1}) I want to replace last equation in 2nd and 4th and then solve the system. r^k,E1 also are given. Do you think it is possible to do this system by extracting the matrices and using the idea of a Lagrange multiplier?

9

I guess I’m not real sure on the notation - is eqn 4 for r linear in the unknowns phi and n?
That is, is phi bar known? Offhand it doesn’t look much different from say a chemical
potential if you are looking for examples. You just need to find the contribution of the
term for each dof ( apparently phi and n ) to the integral over Omega to make your eqn?
Regardless of how complicated it actually is you can always linearize although maybe
the question is if there are easier to user solvers.

10

Thanks for your answer. In fact, the system is linear. With the scheme I used, I did linearization. The unknowns are \phi,\bar{\mu_{\phi}},n,\bar{\mu_n} at time level k+1 and all is known up to time step k. In the last equation r is linear in the unknowns \phi and n (at time step k+1). the rest of the notations are known. just for you to know the notations are as follows:

image
image950×352 56.2 KB

All the things in the time step k+1 are unknowns. I just have a problem with the last equation which is not in the weak form. I implemented the code in Freefem++ but when I replace \bar{r} in the second and fourth equation with the last integral it will be a complicated relation with two nested integrals that it seems Freefem can not handle it.

11

I am quite sure that this can be done in FreeFEM, but it may require a somewhat advanced level of knowledge to actually implement. If you can come up with a simple MWE, perhaps more advanced users on the forum can help you.

1 Like
12

“This is my FreeFem code, and I’m encountering an error in the macro definitions of EXR1 and EXR2. The primary section of the code begins with ‘////////////////////////////////////////////////Second Order scheme() ,’ followed by the problem definition > problem PHITMUMUN(phiT, muphibar, n, munbar, phiTb, muphib, nb, munb) .
Prior to this section, I have defined the finite element spaces required for solving the main problem and have also computed two initial values necessary for the primary scheme. Subsequently, I implement the problem-solving process using a for loop for each time step. I would greatly appreciate any assistance with this issue.”

//initials
int n1=8;
load "lapack";
//border
border b (t =-1,1) { x =t ; y =-1 ; label =2;}
border c (t =-1 ,1) { x =1 ; y =t ; label =3;}
border d1 (t =1,-1) { x =t ; y =1; label =4;}
border e1(t =1,-1) { x =-1;y =t ; label =5;}
mesh Th100 = buildmesh(b(n1)+c(n1)+d1(n1)+e1(n1));
//plot ( Th100 , ps=" wholedomainphi.eps ",cmm="square",value = true,wait = true);
//savemesh(Th100,"Th110.msh"); // Save the Mesh
//mphi = M*phi^2, Mn = M
real p0=300,epsil=0.02,error = 0.001,dt = 0.00001,T =0.0001,delta = 0.01,Gamma = 0.25,chi0=0.05,M=10,threshold = 0.000001,X0=0.25;
real numadapt = 0,c0 = 0.25,alfa = 3;
fespace Vh100(Th100,P2);
Vh100 h=hTriangle;
macro df(u) 4*Gamma*u^3-6*Gamma*u^2+2*Gamma*u      //
macro f(u) Gamma*u^2*(1-u)^2   //      //
macro P(s) p0*s*delta*(s>=0)+0*(s<0)  //
macro  MM(u) M*u^2        //
macro chi(u,v) -X0*u*v  //
macro chiu(u,v) -X0*v//
macro chiv(u,v) -X0*u//
macro E1(u,v) int2d(Th100)(Gamma*u^2*(1-u)^2+chi(u,v))//
macro be1(u,v) (df(u)+chiu(u,v))/sqrt(E1(u,v)+c0)  //
macro be2(u,v) ( chiv(u,v) )/sqrt(E1(u,v)+c0)   //
macro EXR1(u,v,z) int2d(Th100)((be1(u,v))*z)//
macro EXR2(u,v,z) int2d(Th100)((be2(u,v))*z)//

Vh100 phiTp,phiT,phiTb,muphi,muphipp,muphib,muphip,muphibar,n,np,npp=1,nb,mun,munp,munpp,munbar,munb,phitilda,ntilda,muphitilda,muntilda,Ptilda,MMtilda;
Vh100  rp,rtildap;
Vh100 phiTpp =tanh((0.15-sqrt(x^2+y^2))/(sqrt(2)*epsil));
Vh100 rpp = sqrt(E1(phiTpp,npp)+c0);
plot(Th100,phiTpp,wait = true,ps ="phiTpp.eps",cmm = "phiTpp",fill = true,value = true);
plot(Th100,npp,wait = true,ps ="npp.eps",cmm = "npp",fill = true,value = true);
//muphipp
problem MUphiPP(muphipp,muphib)=int2d(Th100)(muphipp*muphib)
-int2d(Th100)((df(phiTpp))*muphib)-int2d(Th100)((chiu(phiTpp,npp))*muphib)
-int2d(Th100)(epsil^2* (dx(phiTpp)*dx(muphib)+dy(phiTpp)*dy(muphib)));
//
//munpp
problem MUnPP(munpp,munb)=int2d(Th100)(munpp*munb)-int2d(Th100)(chiv(phiTpp,npp)*munb)
-int2d(Th100)((1/delta)*npp*munb);
MUphiPP;
MUnPP;
plot(Th100,muphipp,wait = true,ps ="muphipp.eps",cmm = "muphipp",fill = true,value = true);
plot(Th100,munpp,wait = true,ps ="munpp.eps",cmm = "munpp",fill = true,value = true);

//find the second initial of the system(The one which is according to the final system);
//I used perturbation for the second initial means purturb the first initial very small and take the result as a second initial but I will replace it with an appropriate method later on
phiTp = phiTpp + 0.00001;
muphip = muphipp + 0.00001;
np = npp+0.00001;
munp = munpp+0.00001;;
////////////////////////////////////////////////Second Order scheme
problem PHITMUMUN(phiT,muphibar,n,munbar,phiTb,muphib,nb,munb) = 
//phiT which is volume fraction of tumor cells equation
int2d(Th100)(phiT*phiTb)
-int2d(Th100)(phiTp*phiTb)
+int2d(Th100)(dt*MMtilda*(dx(muphibar)*dx(phiTb)+dy(muphibar)*dy(phiTb)))
-int2d(Th100)(dt*Ptilda*munbar*phiTb)
+int2d(Th100)(dt*Ptilda*muphibar*phiTb)
//muphi equation
+int2d(Th100)(muphibar*muphib)
-int2d(Th100)(0.25*epsil^2*(dx(phiTpp)*dx(muphib)+dy(phiTpp)*dy(muphib)))
-int2d(Th100)(0.75*epsil^2*(dx(phiT)*dx(muphib)+dy(phiT)*dy(muphib)))
-int2d(Th100)(rp*(be1(phitilda,ntilda))*muphib)
-int2d(Th100)(0.25*(EXR1(phitilda,ntilda,phiT))*(be1(phitilda,ntilda))*muphib)
+int2d(Th100)(0.25*(EXR1(phitilda,ntilda,phiTp))*(be1(phitilda,ntilda))*muphib)
-int2d(Th100)(0.25*(EXR2(phitilda,ntilda,n))*(be1(phitilda,ntilda))*muphib)
+int2d(Th100)(0.25*(EXR2(phitilda,ntilda,np))*(be1(phitilda,ntilda))*muphib)
//nutrition rich watter equation
+int2d(Th100)(n*nb)
-int2d(Th100)(np*nb)
+int2d(Th100)(M*dt*(dx(munbar)*dx(nb)+dy(munbar)*dy(nb)))
+int2d(Th100)(dt*Ptilda*munbar*nb)
-int2d(Th100)(dt*Ptilda*muphibar*nb)
//mun
+int2d(Th100)(munbar*munb)
-int2d(Th100)(rp*(be2(phitilda,ntilda))*munb)
-int2d(Th100)(0.25*(EXR1(phitilda,ntilda,phiT))*(be2(phitilda,ntilda))*munb)
+int2d(Th100)(0.25*(EXR1(phitilda,ntilda,phiTp))*(be2(phitilda,ntilda))*munb)
-int2d(Th100)(0.25*(EXR2(phitilda,ntilda,n))*(be2(phitilda,ntilda))*munb)
+int2d(Th100)(0.25*(EXR2(phitilda,ntilda,np))*(be2(phitilda,ntilda))*munb);

for(int i=0;i<=T/dt;i+=1){
phitilda = 1.5*phiTp-0.5*phiTpp;
Ptilda = 1.5*(P(phiTp))-0.5*(P(phiTpp));
MMtilda = 1.5*(MM(phiTp))-0.5*(MM(phiTpp));
rp = sqrt(E1(phiTp,np)+c0);
rtildap = sqrt(E1(phitilda,ntilda)+c0);
PHITMUMUN;
real AA = dx(phiT)^2+dy(phiT)^2+alfa*dx(n)^2+alfa*dy(n)^2;
if(AA>threshold){
numadapt=numadapt+1;
Th100 = adaptmesh(Th100,AA);
phitilda = 1.5*phiTp-0.5*phiTpp;
ntilda = 1.5*np - 0.5*npp;
Ptilda = 1.5*(P(phiTp))-0.5*(P(phiTpp));
MMtilda = 1.5*(MM(phiTp))-0.5*(MM(phiTpp));
rp = sqrt(E1(phiTp,np)+c0);
PHITMUMUN
;};
muphi[] = 2*muphibar[]-muphip[];
mun[] = 2*munbar[]-munp[];
phiTpp[]= phiTp[];
phiTp[] = phiT[];
muphipp[] = muphip[];
muphip[] = muphi[];
npp[] = np[];
np[] =n[]; 
munpp[] = munp[];
munp[] = mun[];
cout<<"number of iteration="<<i<<endl;
cout<<"Metric="<<AA<<endl;
//***if(i==10 | i==50 | i==100 | i == 500){
plot(Th100,phiTp,wait = true,ps ="phiTp.eps",cmm = "phiTp",fill = true,value = true);
//plot(phiDp,wait = true,ps ="phiDp.eps",cmm = "phiDp",fill = true,value = true);
//plot(mup,wait = true,ps ="mup.eps",cmm = "mup",fill = true,value = true);
//plot(Th100,pp, wait=true,ps="pp.ps",cmm="pp",fill = true,value = true);
//plot(Th100,np, wait=true,ps="np.ps",cmm="np",fill = true,value = true);
////////****cout<<"phip="<<phip[]<<endl;
//plot(mup,wait = true,ps ="mup.eps",cmm = "mup",fill = true,value = true);
/////***cout<<"mup="<<mup[]<<endl;
//plot(pp, wait=true,ps="pp.ps",cmm="pp",fill = true,value = true);
///////***cout<<"pp="<<pp[]<<endl;
//};
};
plot(Th100,phiTp,wait = true,ps ="phiT.eps",cmm = "phiT",fill = true,value = true);
cout<<"phiT="<<phiTp[]<<endl;
plot(Th100,muphi,wait = true,ps ="muphi.eps",cmm = "muphi",fill = true,value = true);
cout<<"muphip="<<muphip[]<<endl;
plot(Th100,np, wait=true,ps="n.ps",cmm="n",fill = true,value = true);
cout<<"np="<<np[]<<endl;
plot(Th100,munp, wait=true,ps="mun.ps",cmm="mun",fill = true,value = true);
cout<<"munp="<<munp[]<<endl;
cout<<"size of mesh =" <<h[].max<<endl;
cout<<"number of adaptation="<<numadapt<<endl;
ofstream file("phiT40.txt");
file << phiT[];
13

Please upload your code as an attachment or enclose it in ``` to format it correctly.

14

I’ve made the necessary edits. Thank you for bringing it to my attention.