I’ve run into a situation recently where I’m solving a nonlinear PDE; to do this, I am writing a Newton solver, e.g. solving a system J p = F. In this case, F is a vector for my residual and J is my Jacobian matrix. When creating these objects, I must use tgv = 0 for F otherwise my approximate solution becomes unbounded. Even when I use the exact solution of my problem as the initial guess for Newton’s method, my approximate solution will continue to grow after each iteration. (Note however, that the tgv for J does not matter, as any value will work as long as tgv for F is 0.)
I’ve tried reading some of the documentation as well as other discussions on here about the tgv, but I guess I’m still not sure theoretically why the tgv value of 0 works in this case. Does anyone have any thoughts on this?
Hopefully this is clear enough. Please let me know if any other information would be helpful. Thank you!
Hi,
Setting tgv=0 for the right-hand side means that you put zero as Dirichlet value.
I suspect that you do not manage correctly the boundary condition related to your Newton method.
Consider that you have to solve J(u_{k+1}-u_k)+F(u_k)=0, with J the jacobian matrix of F at u_k. There are two ways to do it: either you solve it in the unknown u_{k+1}, or in the unknown \delta u_{k+1}\equiv u_{k+1}-u_k. Let me consider the second choice, that seems to correspond to what you say.
Assume that you want to set the Dirichlet boundary condition u=u_b on the boundary.
Then if you solve the problem in the unknown \delta u_{k+1}, the boundary condition that you have to set is \delta u_{k+1}=u_b-u_k on the boundary. If by chance u_k has the right BC (u_k=u_b), the right-hand side becomes 0. This is maybe what you do by setting tgv=0. On the contrary if you set \delta u_{k+1}=u_b on the boundary, with tgv non zero, this is wrong (except if u_b=0).
The correct way should be to set the BC \delta u_{k+1}=u_b-u_k on the boundary, with tgv=-1 for matrix and right-hand side.
Hi @fb77 , thank you very much for the explanation. Indeed, as you have described, I am trying to solve a system J (\delta u_{k+1)} = -F(u_k).
I tried setting the boundary to u_b - u_k and then using tgv=-1 as you suggested, however my approximate solution still becomes unbounded after each Newton iteration. As before, the only time my algorithm converges is when the tgv for the right-hand side is 0. I’m not sure if perhaps I have misunderstood what you said or if there is some bug in my code. I am copying my code here, in case that is helpful.
int n = 100;
real k1 = 0.1;
real k2 = -0.16;
int Ftgv = -1;
int Jtgv = -1;
//define the boundary: R in [0,1] and z in [−1,1]
border Gamma1(t=0., 1.){x = t; y = -1;}
border Gamma2(t=0., 1.){x = 1; y = -1+2t;}
border Gamma3(t=0., 1.){x = 1-t; y = 1;}
border Gamma4(t=0., 1.){x = 0; y = 1-2t;}
mesh M = buildmesh(Gamma1(n) + Gamma2(n) + Gamma3(n) + Gamma4(n));
int[int] lab = labels(M);
//define finite element space
fespace Vh(M, P2);
Vh psi, phi, psiExact, error;
//exact solution
func Exact = -60./(9. * x^2 + (y+2.)^2);
psiExact = Exact;
//define variational forms
varf vJ(w,phi) = int2d(M)(dx(w)dx(phi) + dy(w)dy(phi) + dx(w)phi/x - wpsiphi(2k1+3k2x^2psi)) + on(lab, w=psiExact-psi);
varf vF(w,phi) = int2d(M)(dx(psi)dx(phi) + dy(psi)dy(phi) + dx(psi)phi/x - (psi^2(k1+k2x^2psi))phi) + on(lab, w=psiExact-psi);
//initial guess
psi = Exact(1. + x*(x-1.)(z+1.)(z-1.));
plot(psiExact, wait=true, cmm=“Exact Solution”, value=true);
//Newton’s Method
int iter = 20;
real tol = 1.0e-12;
for(int i=1; i < iter; ++i){
matrix J = vJ(Vh, Vh, tgv = Jtgv);
real[int] F = vF(0, Vh, tgv = Ftgv);
real[int] inc = J^-1 * F;
psi -= inc;
plot(psi, wait=true, cmm=“Solution at iteration #”+i, value=true);
if(F.linfty < tol) break;
}
error = abs(psi - Exact);
real L2error = error.l2;
cout << "L2error = " << L2error << endl;
You are absolutely correct, my sincerest apologies for this silly mistake. After changing the boundary condition, my code is now working as it should be. Additionally, the L^2 error is much more reasonable than before.
Thank you for taking the time to help me!