Dear all,

I have attached a demo code for incompressible two phase flow using hybrid-mixed FEM (RT0) for pressure and DG-P0 for saturation.

Currently, I use the standard DG updates for the saturation conservation law equation.

I believe this default version could be significantly improved, but I don’t know what’s the best option available now in ngsolve. I would like to get ideas about speeding up this code. Many thanks as always!

Best,

Guosheng

Attachment: quarterFiveSpot2D_compile.ipynb

Hi Guosheng,

I tried to modify your ForwardEuler by replacing the Matrix-Vector multiplication, but this did not really help.

```
def ForwardEuler(sh, dt):
rhsS.data = aS.mat * sh
#rhsS.data = invm*rhsS
fesS.SolveM(rhsS)
return sh-dt*rhsS
```

Also to “inline” it did not bring a lot

```
#gfuS.vec.data = ForwardEuler(gfuS.vec, dt)
aS.Apply(gfuS.vec,rhsS)
gfuS.vec.data -= dt*invm*rhsS
```

However, the Compile flag made a huge difference! It is important to compile everything, not only a part of it

```
if kdg>0:
aS += (-uh*fw(ss)*grad(tt)).Compile(true_compile, wait=True)*dx
...
aS += (uh*n*Fhat*tt).Compile(true_compile, wait=True)*dx(element_boundary=True)
```

With this I got a speed up of a factor 3!

It is also possible to speed up your “fast” part even more. Instead of calling mat.Invert every time you can “update” the inverse:

```
a.Assemble()
inv = a.mat.Inverse(fesP.FreeDofs(True),inverse="sparsecholesky")
while time < tend-dt/2 :
if stepS%nsubcycles == 0:
a.Assemble()
inv.Update()
...
gfuP.vec.data = inv * rhsP
```

Then the symbolic factorization is done only once and not every step. This improved the “fast” part from 2.4 to 1.5 seconds on my machine.

Best

Michael

Hi Michael,

Thanks about the two tricks, they are helpful.

I also get roughly a factor of 3 speedup with the correct use of the Compile flag, and a factor of 2 for the sparse Cholesky inverse update.

Changing element_boundary formulation to skeleton-based formulation further helps to speed-up the explicit DG solver slightly.

I am now thinking about using implicit DG solver for saturation. But Newton solver fails since .Other() only works for proxyfunctions, not gridfunctions. Do you have any suggestions on this issue?

Best,

Guosheng