Point values of unknown fespace functions in variational form

General Discussion
1

Hello,

I am iteratively solving a flow problem with a nonlinear forcing term which is dependent on the value of the solution at a boundary. I use Newton’s method to solve, but as I have not yet figured out how to implement the linearized forcing term in FreeFEM, I have not been able to achieve quadratic convergence.

With some simplifications to keep this clear, suppose my nonlinear problem is N(u) with Jacobian J(u, du) and with a forcing term F(u)=ub^2, where ub is the value of u along the boundary x=x0. At each iteration, I form and solve Ax=b, where the residual (b) is computed by evaluating the forcing term (which is defined along the boundary) at every point in the domain.

mesh th; fespace Xh(th, P2); Xh u, du, v; // preliminary stuff

varf b(du, v) = int2d(th)( N(u) * v + u(x0, y)^2 * v ); // varf for residual vector
real[int] bb = b(0, Xh); // build residual vector

I would like to use the same approach to build the sparse matrix (A) with a linearized version of the forcing term. Ideally, I could do this as:

varf A(du, v) = int2d(th)( J(u, du) * v + 2. * u(x0, y) * du(x0, y) * v ); // varf for jacobian matrix

However, I get an error with this. It seems I cannot evaluate the point values of the unknown fespace function du in the varf. Has anyone encountered this before? Any workarounds?

Cheers!
-Chris

2

Since I have not yet received a response, I have decided to give a more detailed working example of my problem. This is a somewhat minimal example using the steady Burgers’ equation on a square mesh. It is designed to demonstrate my issue and not to express any real physics. The core issue here is how to implement the linearization of the forcing term (if this is even possible) since its global strength is determined by its value along a boundary.

// A simple example of my problem using Burgers' equation
mesh th = square(25, 25);fespace Xh(th,[P2, P2]); // preliminaries
Xh [ux, uy], [dux, duy], [vx, vy]; // definitions

real nu = 0.005; // viscosity parameter
real f; // forcing strength

// Nonlinear operator (chosen here to be Burgers' equation)
macro N(u,v)
  (u#x*dx(u#x)*v#x + u#y*dy(u#x)*v#x + nu*(dx(u#x)*dx(v#x) + dy(u#x)*dy(v#x))
   + u#x*dx(u#y)*v#y + u#y*dy(u#y)*v#y + nu*(dx(u#y)*dx(v#y) + dy(u#y)*dy(v#y))) // End macro
// Nonlinear forcing term dependent on the value at a boundary
macro F(u,v)
  (f*(u#y^2 - u#y(1,y)^2)*v#x) // End macro
// Jacobian of nonlinear operator
macro dN(du,u,v)
  (du#x*dx(u#x)*v#x + du#y*dy(u#x)*v#x + u#x*dx(du#x)*v#x + u#y*dy(du#x)*v#x + nu*(dx(du#x)*dx(v#x) + dy(du#x)*dy(v#x))
   + du#x*dx(u#y)*v#y + du#y*dy(u#y)*v#y + u#x*dx(du#y)*v#y + u#y*dy(du#y)*v#y + nu*(dx(du#y)*dx(v#y) + dy(du#y)*dy(v#y))) // End macro
// Jacobian of nonlinear forcing term
macro dF(du,u,v)
  (f*2.*du#y*u#y*v#x) // End macro, not sure how to implement this
  // NOTE: the linearization dF is incomplete as it does not include the boundary term
  
// Boundary conditions
macro BC(du,u) (on(1,3,dux=-ux,duy=-uy) + on(4,dux=y*(1-y)-ux,duy=y*(1-y)-uy)) // End macro

// Declare variational forms
varf b([dux, duy], [vx, vy]) = int2d(th)(-(N(u,v) + F(u,v))) + BC(du,u);
varf A([dux, duy], [vx, vy], solver=sparsesolver) = int2d(th)(dN(du,u,v) + dF(du,u,v)) + BC(du,u);

// Newton iteration on linearized system with f=0
ux[] = 0.;
f = 0.;
real error = 1.0;
int i = 0;
cout << " ------ UNFORCED SYSTEM ------" << endl;
while (error > 1e-9){
  real[int] bb = b(0,Xh); // evaluate residual
  matrix AA = A(Xh,Xh); // build matrix
  dux[] = AA^-1*bb; // solve for correction
  ux[] += dux[]; // update solution
  error = dux[].l2/ux[].l2; // compute relative change
  i++;
  cout << " Iteration = " << i << ", error = " << error << endl;
}
cout << " Converged with " << i << " iterations." << endl;



// Newton iteration on linearized system with f=1
ux[] = 0.;
f = 1.;
error = 1.0;
i = 0;
cout << "------ FORCED SYSTEM ------" << endl;
while (error > 1e-9){
  real[int] bb = b(0,Xh); // evaluate residual
  matrix AA = A(Xh,Xh); // build matrix
  dux[] = AA^-1*bb; // solve for correction
  ux[] += dux[]; // update solution
  error = dux[].l2/ux[].l2; // compute relative change
  i++;
  cout << " Iteration = " << i << ", error = " << error << endl;
}
cout << " Converged with " << i << " iterations." << endl;
3

Special thanks to @prj for helping with this problem. It is now resolved and I will leave the solution below in case anyone else has a similar question. The key here was to manipulate the Jacobian matrix using a low-rank update. This update uses a projection onto a “1D” fespace that shares the boundary governing the forcing and then projects the result back to the full fespace. This approach restricts the operator to the boundary degrees of freedom for the unknowns but not for the test functions. I hope someone finds this useful! Cheers.

   // A simple example of my problem using Burgers' equation
func Pk = [P2, P2];
macro def(u) [u#x, u#y] // End macro
    
mesh th = square(25, 25), th2 = square(1, 25, [1 + 0.01*x, y]); // define meshes

fespace Xh(th, Pk), Xh2(th2, Pk, periodic = [[2, y], [4, y]]); // define spaces
Xh def(u), def(du), def(v); // definitions

matrix Ih  = interpolate(Xh, Xh2); // from skinny to normal fespace
matrix Ih2 = interpolate(Xh2, Xh);  // from normal to skinny fespace
matrix II  = Ih*Ih2;
real f, tol = 1e-10, nu = 0.008; // forcing strength, tolerance, viscosity parameter

// Nonlinear operator (chosen here to be Burgers' equation)
macro N(u,v)  (u#x*dx(u#x)*v#x + u#y*dy(u#x)*v#x + nu*(dx(u#x)*dx(v#x) + dy(u#x)*dy(v#x)) + u#x*dx(u#y)*v#y + u#y*dy(u#y)*v#y + nu*(dx(u#y)*dx(v#y) + dy(u#y)*dy(v#y))) // End macro
// Nonlinear forcing term dependent on the value at a boundary
macro F(u,v)  ((u#y^2 - f*u#y(1, y)^2)*v#x) // End macro
// NOTE: the boundary term is proportional to -f
// Jacobian of nonlinear operator
macro dN(du,u,v)  (du#x*dx(u#x)*v#x + du#y*dy(u#x)*v#x + u#x*dx(du#x)*v#x + u#y*dy(du#x)*v#x + nu*(dx(du#x)*dx(v#x) + dy(du#x)*dy(v#x))
   + du#x*dx(u#y)*v#y + du#y*dy(u#y)*v#y + u#x*dx(du#y)*v#y + u#y*dy(du#y)*v#y +     nu*(dx(du#y)*dx(v#y) + dy(du#y)*dy(v#y))) // End macro
// Jacobian of nonlinear forcing term
macro dF(du,u,v)  (2.*du#y*u#y*v#x) // End macro, not sure how to implement this
// NOTE: the linearization dF is incomplete as it does not include the boundary term

// Boundary conditions
macro BC(du,u) (on(1, 3, dux = -ux, duy = -uy) + on(4, dux = y*(1 - y) - ux, duy = y*(1 - y) - uy)) // End macro
// Declare variational forms
varf b(def(du), def(v)) = int2d(th)(-(N(u, v) + F(u, v))) + BC(du, u);
varf A(def(du), def(v)) = int2d(th)(dN(du, u, v) + dF(du, u, v)) + BC(du, u);
varf C(def(du), def(v)) = int2d(th)(-2.*f*uy(1, y)*duy*vx); // linearized forcing term

// Newton iteration on linearized system with f=0
ux[] = 0.; f = 0.; real error = 1.0; int i = 0;
cout << " ------ UNFORCED SYSTEM ------" << endl;
while (error > tol){
  real[int] bb = b(0, Xh); // evaluate residual
  matrix AA = A(Xh, Xh, solver = sparsesolver); // build matrix
  dux[] = AA^-1*bb; // solve for correction
  ux[] += dux[]; // update solution
  error = dux[].l2/ux[].l2; // compute relative change
  i++;  cout << " Iteration = " << i << ", error = " << error << endl;
}cout << " Converged with " << i << " iterations." << endl;

// Newton iteration on linearized system with f=5
ux[] = 0.; f = 5.; error = 1.0; i = 0;
cout << "------ FORCED SYSTEM (boundary term not linearized) ------" << endl;
while (error > tol){
  real[int] bb = b(0, Xh); // evaluate residual
  matrix AA = A(Xh, Xh, solver = sparsesolver); // build matrix
  dux[] = AA^-1*bb; // solve for correction
  ux[] += dux[]; // update solution
  error = dux[].l2/ux[].l2; // compute relative change
  i++;  cout << " Iteration = " << i << ", error = " << error << endl;
}cout << " Converged with " << i << " iterations." << endl;

// Newton iteration on linearized system with f=5
ux[] = 0.; f = 5.; error = 1.0; i = 0;
cout << "------ FORCED SYSTEM (fully linearized) ------" << endl;
while (error > tol){
  real[int] bb = b(0, Xh); // evaluate residual
  matrix AA = A(Xh, Xh); // build matrix
  matrix CC = C(Xh, Xh); // build matrix
  CC = CC*II;
  AA += CC;
  set(AA, solver = sparsesolver);
  dux[] = AA^-1*bb; // solve for correction
  ux[] += dux[]; // update solution
  error = dux[].l2/ux[].l2; // compute relative change
  i++;  cout << " Iteration = " << i << ", error = " << error << endl;
}cout << " Converged with " << i << " iterations." << endl;