3
$\begingroup$

In a Newton iterative FEM PDE solver I am writing, each iterative step involves updating the current solution function u by adding to it the increment function du, which is the solution of some linearised PDE whose coefficients depend on the current u. Both u and du are InterpolatingFunction.

Now, as the number of iterates n grows, it becomes quite inefficient to evaluate the sum of n interpolating functions.

How can two InterpolatingFunction over a finite element mesh be added together to get a single InterpolatingFunction, whilst keeping all the original unstructured mesh properties / interpolation order?

ClearAll["Global`*"]
<< NDSolve`FEM`
rad = 5; th = Pi/4;
reg = Disk[{0,0}, rad, {-th,th}];
mesh = ToElementMesh[reg];

Say the two functions to be added are

f1 = NDSolveValue[{Laplacian[u[x,y],{x,y}] == 1^2 * u[x,y]
       + NeumannValue[1, True]}, u, Element[{x,y}, mesh]];
f2 = NDSolveValue[{Laplacian[u[x,y],{x,y}] == 2^2 * u[x,y]
       + NeumannValue[1, True]}, u, Element[{x,y}, mesh]];

Attempt 0.

sum0 = Function[{x,y}, f1[x,y] + f2[x,y] // Evaluate];
Do[  f1[0,0], 100000] // AbsoluteTiming  (* {0.58909, Null} *)
Do[sum0[0,0], 100000] // AbsoluteTiming  (* {1.64215, Null} *)

Using Function obviously doesn't lump the sum into a single InterpolatingFunction.

Attempt 1.

sum1 = NDSolveValue[s[x,y] == f1[x,y] + f2[x,y], s, Element[{x,y}, mesh]];

NDSolveValue::derivs: No derivatives of dependent variables were found in the equations. NDSolveValue is designed to solve differential or differential algebraic equations. Use NSolve or FindRoot to numerically solve algebraic equations.

See (23416); it's a pity you can't solve "zeroth order differential equations" anymore.

Attempt 2.

sum2 = FunctionInterpolation[f1[x,y] + f2[x,y], Element[{x,y}, mesh]];

FunctionInterpolation::range: "Argument {x,y}[Element]mesh is not in the form of a range specification, {x, xmin, xmax}."

Using the form {x, xmin, xmax}, ... is not an option because the reg is not rectangular and mesh is unstructured, and I want the sum to inherit these properties.

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
4
$\begingroup$

You can do something like this (assuming of course that the InterpolatingFunction f1 and f2 are generating using the same mesh):

f12 = ElementMeshInterpolation[{mesh}, f1["ValuesOnGrid"] + f2["ValuesOnGrid"]]

Where ElementMeshInterpolation is from the 

NDSolve`FEM`

package.

f1[1, 1] + f2[1, 1]
(* -1.84159 *)

f12[1, 1]
(* -1.84159 *)

Do[f1[0, 0], 100000] // AbsoluteTiming
(* {0.223145, Null} *)

Do[f12[0, 0], 100000] // AbsoluteTiming
(* {0.213769, Null} *)
Improve this answer
$\endgroup$
2
$\begingroup$

Attempt 3.

It turns out you can trick Mathematica into not complaining about lack of derivatives by adding a dependent derivative in a separate equation:

sum3 = NDSolveValue[
         {Laplacian[dummy[x, y], {x, y}] == 0,
          s[x, y] == f1[x, y] + f2[x, y]
         }, s, Element[{x, y}, mesh]]

Of course this is slower than chuy's method.

Improve this answer
$\endgroup$