Integrating Nonlinear Terms

Hello,

I am learning to write FreeFem code for the Cahn-Hilliard equation with Neumann boundary conditions. I have simplified one of the nonlinear terms, but it seems that the int2d function requires a bilinear form. I believe this is where my problem lies. How should I write the integration for a nonlinear form?


Bao_h1_Neu3.edp (3.0 KB)

Lina

You cannot write a bilinear varf for a nonlinear term.
Instead you have to write a linear varf

varf bphi(unused,zz) = int2d(Th)(Fhat()*zz/epsilon);
real[int] b2 = bphi(0,Vh);

Here the value of phi is taken via the definition of Fhat()

Also, the definition of A should be
matrix A = [[Amumu,Amuphi],[Aphimu,Aphiphi]];

It leads to a blow up, probably because your nonlinearity is taken explicitly.
Instead you should use an implicit scheme (put \phi^{n+1} instead of \phi^n in F_\phi), but this has to be resolved by an iteration scheme like the Newton method, as is done in
Cahn-Hilliard.edp (3.2 KB)
This is a simplification to P1 from https://community.freefem.org/t/can-you-tell-me-what-is-the-error-on-this-codes-progrrame-running-still-exact-and-approximate-solution-not-matching/3220/16 where Discontinuous Galerkin is used with P1dc

1 Like

Thank you very much for your guidance. I have reviewed the code Cahn-Hilliard.edp you shared. However, I encountered an issue during the Newton iteration process. It seems that the Jacobian matrix varf vhJ in the code should be derived as varf vdJ by differentiating with respect to the variables. However, I am unable to derive varf vhJ. Could you kindly provide some guidance on this?

I agree with your equations. You have to solve
\int 2\frac{v^{n+1/2}-v^n}{\delta t}\phi+\int\nabla z^{n+1/2}\cdot\nabla\phi-\int g^{n+1/2}\phi=0,
\int -z^{n+1/2}\psi+\int a\nabla v^{n+1/2}\cdot\nabla\psi+\int F'(v^{n+1/2})\psi=0.
The Newton method defines a sequence (v_k,z_k) such that as k\to\infty, one has v_k\to v^{n+1/2}, z_k\to z^{n+1/2}. The sequence is defined through
\int 2\frac{v_{k+1}-v^n}{\delta t}\phi+\int\nabla z_{k+1}\cdot\nabla\phi-\int g^{n+1/2}\phi=0,
\int -z_{k+1}\psi+\int a\nabla v_{k+1}\cdot\nabla\psi+\int \bigl(F'(v_k)+F''(v_k)(v_{k+1}-v_k)\bigr)\psi=0.
Denoting \delta v=v_{k+1}-v_k, \delta z=z_{k+1}-z_k it can be written
\int 2\frac{\delta v}{\delta t}\phi+\int\nabla \delta z\cdot\nabla\phi=-\Bigl(\int 2\frac{v_k-v^n}{\delta t}\phi+\int\nabla z_{k}\cdot\nabla\phi-\int g^{n+1/2}\phi\Bigr),
\int -\delta z\psi+\int a\nabla \delta v\cdot\nabla\psi+\int \bigl(F''(v_k)\delta v\bigr)\psi=-\Bigl(\int -z_{k}\psi+\int a\nabla v_{k}\cdot\nabla\psi+\int F'(v_k)\psi\Bigr).
In this system the linear operator on the left-hand side is the jacobian matrix (defined in the code by vhJ), and the right-hand side is minus the residual.
The instruction auxv[]=vdJ(0,Vh2) computes the residual.
Then vv[]=H^-1*auxv[] computes vv=-(\delta v,\delta z).
Finally oldv[]-=vv[] (or equivalently oldv[]=oldv[]-vv[]) computes (v^{k+1},z^{k+1})=(v^k,z^k)+(\delta v,\delta z).

1 Like

Thank you so much for your guidance. Previously, when using Newton’s iteration, I always computed the Jacobian matrix through differentiation and then determined the solution’s progression direction. Your approach has been truly enlightening to me, and I deeply appreciate it!

Indeed a trick for Newton’s method is as follows.

When you have to solve F(X)=0 (X and F(X) are vectors), the Newton method writes
X_{k+1}-X_k=-dF(X_k)^{-1}F(X_k).
You can write it as
F(X_k)+dF(X_k)(X_{k+1}-X_k)=0\qquad (1)
In other words, instead of solving F(X_{k+1})=0, you replace F by its first order expansion around X_k, it gives (1).

This is what I wrote in my previous message: starting from the equations on (v^{n+1/2},z^{n+1/2}) I derive the equations on (v_{k+1},z_{k+1}) by making a first order expansion around (v_k,z_k).