Dear all experts and respected professors, i am wondering why the following simple equation does not giving good results>
Here, is my codes:
Heat_Elliptic_combined.edp (2.7 KB)
Thanks in advance!.
There are several mistakes in your code:
t+=dteps*z*phi see the FF documentation https://doc.freefem.org/pdf/FreeFEM-documentation.pdf, note at the bottom p.218z, one possibility is to decouple the problem (at each timestep): first solve the equation on u, then once u and its average can be computed, solve the equation for z. However this works only if there is no dependence in z in the equation on u (this is true for your example).z in the equation on u) and you want to solve the problem in a coupled manner (solving u and z simultaneously), you can use the Lagrange multiplier method. This means that you take the average of u as additional scalar unknown.Thanks you so much for your responses. It helps me a lot.
Can this Lagrange multiplier method is applicable if the equation is looks like the following:
Many many Thanks sir!.
A priori yes, but the stability of the scheme is quite much unclear. You have to use \bar\phi instead of \bar\phi_0 to be sure that the right-hand side of (2.1c) has vanishing integral at the discrete level and for all times.
There is a condition missing in your system to fix a possible additive constant in \xi when solving (2.1c) with Neumann BC. But it may be unnecessary since adding a constant to \mu via (2.1b) will not modify (2.1a) nor (2.1d). The most simple is to set \int\xi=0 (this is what is taken in the previous code).
So, i will take average of phi instead of average of phi_{0} in the code sir??.
Actually, i have to solve equation (2.1a) -(2.2c) by coupling.
Yes. For the continuous problem one has \int\phi(t)=\int\phi_0, but for the numerical approximation it is not true exactly. You have to take \int\phi(t).
Okay thank you sir.
One more thing sir that this method can work in DG setting??. (That is, i have to use method in DG setting by replacing all function space in P1dc or P2dc).
I don’t know. You can try.
Yes sir. Thank you so much❤.
It seems that you can somehow simplify your system by setting \mu_1=f(\phi)-\varepsilon\Delta\phi. Then you get
\phi_t+\nabla\cdot(\phi u)-\varepsilon\Delta \mu_1+\varepsilon\theta(\phi-\bar\phi)=0,
\mu_1=f(\phi)-\varepsilon\Delta\phi,
-\Delta\xi=\theta(\phi-\bar\phi),
with BC \partial_n\phi=\partial_n\mu_1=\partial_n\xi=0.
In this form \xi can be computed afterwards. The coupling is only in the unknowns \phi and \mu_1.
This approach is possible as long as you solve the u separately.
If you need to compute all variables \phi, \mu, \xi, u simultaneously (for stability reasons?), then the previous decoupling is not possible because \nabla \mu (and therefore \nabla\xi) is needed in (2.1d).
I think that the resolution of u separately can be stable because it is involved in (2.1a) only in the advection term, which is low order with respect to the \Delta\mu term.
Yes sir. But one paper is there where they solved all the variables simultaneously by using pressure correction sir. They provided the stablity also.
Here, is the paper>
DG PPC of CHDS.pdf (2.0 MB)
Ok, do what you think is best!
Thanks sir. DG setting is also working sir.
Many many thanks for your support.
Respected sir, i have implemented your idea of Lagrange multiplier to solve first three equations (2.1a), (2.1b), and (2.1c) in the following equations:
Here, is a piece of code sir>>
Blockquote
real theta=1.;
real a=0.1;
real ac=1./a;
// Finite Element Spaces
fespace Vh(Th,P2);
Vh u,w,z,phi,psi,rho;
fespace Vh3(Th,[P2,P2,P2]);
Vh3 [wu,ww,wz];
real uAverage;
varf averagephi(u,w,z,phi,psi,rho)=int2d(Th)(phi/Th.area);
real[int] avphiform=averagephi(0,Vh3);
varf averagepsi(u,w,z,phi,psi,rho)=int2d(Th)(psi/Th.area);
real[int] avpsiform=averagepsi(0,Vh3);
varf averagerho(u,w,z,phi,psi,rho)=int2d(Th)(rho/Th.area);
real[int] avrhoform=averagerho(0,Vh3);
// To solve First three equations (2.1a), (2.1b), and (2.1c) simultaneously
varf pbuwz(u,w,z,phi,psi,rho)=int2d(Th)(u*phi/dt+a*Grad(w)'*Grad(phi)) // Due to (2.1a)
-int2d(Th)(w*psi)+int2d(Th)((ac*(ukold^3-ukold)-stab*ukold)*psi)// Non-linear term for (2.1b)
+int2d(Th)(stab*u*psi)//stabilization
+int2d(Th)(a*Grad(u)'*Grad(psi)) // for (2.1b)
+int2d(Th)(z*psi) // for equaution of z in (2.1c)
-int2d(Th)(u0*v10*dx(phi)+u0*v20*dy(phi)) // advection term in (2.1a)
-int2d(Th)(u0*phi/dt+f(t+dt/2.)*phi)// RHS of first equation in (2.1a)
+int2d(Th)(Grad(z)'*Grad(rho)) // for equation of z in (2.1c)
+int2d(Th)(1.e-8*z*rho) // Due to (2.1c)
-int2d(Th)(u*rho); // Due to (2.1c)
// solve
matrix A=pbuwz(Vh3,Vh3,Vh3);
real[int] rhs=pbuwz(0,0,Vh3);
rhs=-rhs;
matrix B=[[A,avrhoform],[avphiform',-1.],[avpsiform',-1.]];
set(B,solver="UMFPACK");
real[int] rhsB(rhs.n+1);
rhsB=[rhs,0., 0.];
real[int] solB=B^-1*rhsB;
[wu[],ww[],uAverage]=solB;
u=wu;
z=wz;
w=ww;
Blockquote
Here, u0, [v10, v20](initial velocity) are knows due to initial datas.
I have used the continuous setting of the paper>>
DG PPC of CHDS.pdf (2.0 MB)
I have not given nonlinear loop and time loop here sir (ukold comes for nonlinearity) (as i know you know better than me sir).
I have tried to solve as in the paper in continuous setting:
But some errors are coming in matrix A.
Here, is the errors>>
I will be very grateful sir if i get some help on this.
Thanks in advance sir!.
It has to be
matrix A=pbuwz(Vh3,Vh3);
real[int] rhs=pbuwz(0,Vh3);
rhs=-rhs;
matrix B=[[A,avrhoform],[avphiform',-1.]];
set(B,solver="UMFPACK");
real[int] rhsB(rhs.n+1);
rhsB=[rhs,0.];
real[int] solB=B^-1*rhsB;
[wu[],uAverage]=solB;
u=wu;
w=ww;
z=wz;
Thank you so much sir for your great help. Now, code is working fine sir.
**But Can you explain in short why **
matrix A=pbuwz(Vh3,Vh3) instead of matrix A=pbuwz(Vh3,Vh3,Vh3);
still there were three unkowns in the bilinear form of pbuwz. I am not getting this sir??.
The declaration is varf pbuwz(u,w,z,phi,psi,rho)
Then A=pbuwz(Vh3,Vh3) means that:
the first Vh3 is for the set of unknowns (u,w,z)\in Vh3,
the second Vh3 is for the test functions (phi,psi,rho)\in Vh3
Okay sir. Very nice explanation.