WIP: Import from FEniCS and simplify code. by BenjaminRodenberg · Pull Request #1 · NiklasVin/tutorials

BenjaminRodenberg · 2023-11-16T09:45:16Z

No description provided.

BenjaminRodenberg · 2023-11-16T09:45:46Z

partitioned-heat-conduction/fenics/utils/utils.py

    # apparently you need 2k+2 samples for degree k
-    x_dist = float(dt) / (2 * degree + 2)
-    nodes = []
+    nodes = np.linspace(0,dt,2*3+3)  # <-- But these are 2k+3 samples


This contradicts the statement above. Which one is true?

The comment about the number of samples needed was a bit imprecise. You need at least 2k+2 samples.
I will check, if you can reduce the samples to at least 2k, as I changed this code segment here:
https://github.com/NiklasVin/tutorials/commit/00007d7be9a4d8334089c2f91b3a25e49ed6b0ec#diff-0b22c9acd2661e0e0fe1d0350fdd7736431cd134ab6f1e5b088d90b0f8c7b2b8L34-R35
If you used the constructor directly instead of using the return values of splrep, the constructor had for some reason a problem less than 2k+2 samples were given. Maybe an implementation error from my side, but I am on it :)

partitioned-heat-conduction/fenics/heat.py

BenjaminRodenberg

Noteworthy observation: If I try to use a different time stepping scheme via tsm = BackwardEuler() or tsm = RadauIIA(1) I get an error. Nothing we have to fix now, but important to know that we cannot use arbitrary time stepping methods for some reason. @NiklasVin do you observe the same problem? Maybe this is also due to my FEniCS installation and some compatibility problem.

BenjaminRodenberg

I think I found the problem 🎉

partitioned-heat-conduction/fenics/heat.py

BenjaminRodenberg · 2023-11-16T12:33:06Z

partitioned-heat-conduction/fenics/heat.py

        for i in range(tsm.num_stages):
            bc.append(DirichletBC(Vbig.sub(i), du_dt[i], remaining_boundary))
-            F += vs[i] * coupling_expressions[i] * dolfin.ds
+            F += vs[i] * coupling_expressions[-1] * dolfin.ds


It looks like we need to use the heat flux corresponding to the end of the current time step for the Neumann boundary condition. I currently do not have a good explanation why this is the case, but if we think about the way how Neumann boundary conditions work in general and compare this to the way Dirichlet work it intuitively makes sense to me. However, there is certainly a better explanation including more hard facts and less intuition 😏

This problem was masked by the bug in the creation of the coupling expressions because due to the pointer-nature of the problem, all coupling expressions were pointing to the coupling expressions corresponding to the end of the current time step.

Another thing we should keep in mind here: If we are using a time stepping scheme, where c[-1] does not correspond to the end of the window, things will probably become tricky. Instead of referring to the coupling_expression via the index, I think we need a different strategy for creating coupling_expressions for the Neumann participant. What's also important: coupling_expression[0] is not needed for the Neumann participant. So why bother creating it in the first place?

Maybe I was a bit quick here: One thing we should definitely test is what happens, if the heat flux changes over time. Currently, we have a constant heat flux and this might hide some possible errors.

In d410f66 I realized that the heat flux was not initialized (no problem with implicit Euler). This results in zero heat flux for t=0, which also leads to an incorrect interpolant. I think this is the much more reasonable fix.

It looks like we need to use the heat flux corresponding to the end of the current time step for the Neumann boundary condition.

Yes, we need to set the heat flux to the end of the time step of the respective stages.
Otherwise it would e.g. not comply with the following statement from the Irksome paper:

We note that weakly-enforced boundary conditions [...] require [...] just the evaluation of any time-dependent data at the correct stage times.

So for the weak form for the stage $k_i$ at $t+c_i\cdot\delta t$.

I tested this as well with $g_{tri}$ (from the Quasi-Newton waveform iteration paper) to check if the convergence of the solution is equal to the order of the time stepping scheme. With the introduction of a time dependent heat flux, wrongly imposed Neumann Boundaries would impact the convergence of the solution.
Until now, I tested it with LobattoIIIC(2) and the results look promising:
Neumann BC at $t+c_i\cdot\delta t$

$\Delta t$ $\delta t$ L2 Dirichlet L2 Neumann

0.1 0.05 1,15E-05 7,09E-05

0.05 0.025 2,59E-06 1,74E-05

0.025 0.0125 6,51E-07 4,32E-06

So the convergence order approximately is, as expected, two.

As I was not quite sure, if you referred to the time of each stage by "end of current time step" or the actual end of the current time step, and as a sanity check, if I understood the idea behind the k-ansatz correctly, I also tested, how the errors look like if you only use data at the boundaries retrieved only from $t+\delta t$. And, thankfully, the results have no convergence order of two (otherwise I would, again, be clueless why this is happening😅):

$\Delta t$ $\delta t$ L2 Dirichlet L2 Neumann

0.1 0.05 2,73E-04 3,98E-04

0.05 0.025 1,59E-04 2,08E-04

0.025 0.0125 8,96E-05 1,07E-04

These are exactly the tests I meant with

One thing we should definitely test is what happens, if the heat flux changes over time.

Thank you! Looks good.

Last but not least, a comment about:

If we are using a time stepping scheme, where c[-1] does not correspond to the end of the window, things will probably become tricky.

I think this should be fine. As

Yes, we need to set the heat flux to the end of the time step of the respective stages. Otherwise it would e.g. not comply with the following statement from the Irksome paper:

We note that weakly-enforced boundary conditions [...] require [...] just the evaluation of any time-dependent data at the correct stage times.

Additionally, tested with the GaussLegendre(1), the time dependent heat flux and the Neumann BC only used at $t+\frac{\delta t}{2}$, the order seems fine.
So it is nice to see, that actually no special treatment for such time stepping methods are required. I hope for the best, that this will also be true for multi step methods.

partitioned-heat-conduction/fenics/heat.py

NiklasVin · 2023-11-16T14:42:03Z

Noteworthy observation: If I try to use a different time stepping scheme via tsm = BackwardEuler() or tsm = RadauIIA(1) I get an error.

Is this still valid for the most recent commit? Using 1b9d9cb both time stepping schemes work fine for me.

With the change, the function shouldn't get into trouble if it created a BSpline with a degree >3

BenjaminRodenberg · 2023-11-17T10:30:02Z

Noteworthy observation: If I try to use a different time stepping scheme via tsm = BackwardEuler() or tsm = RadauIIA(1) I get an error.

Is this still valid for the most recent commit? Using 1b9d9cb both time stepping schemes work fine for me.

Yes. I still get the following error:

Traceback (most recent call last):
  File "heat.py", line 108, in <module>
    tsm = BackwardEuler()
  File "/home/benjamin/Programming/tutorials/partitioned-heat-conduction/fenics/utils/ButcherTableaux.py", line 137, in __init__
    super(BackwardEuler, self).__init__(1)
  File "/home/benjamin/Programming/tutorials/partitioned-heat-conduction/fenics/utils/ButcherTableaux.py", line 127, in __init__
    L = FIAT.GaussRadau(U, num_stages - 1)
AttributeError: module 'FIAT' has no attribute 'GaussRadau'

I think this is probably a problem with the specific version of FEniCS or one of the dependencies I'm using. I'm using the following versions:

pip3 list
...
fenics                       2019.1.0
fenics-basix                 0.4.2
fenics-dijitso               2019.1.0
fenics-dolfin                2019.2.0.dev0
fenics-dolfinx               0.4.1
fenics-ffc                   2019.1.0.post0
fenics-ffcx                  0.4.2
fenics-fiat                  2019.1.0
fenics-ufl                   2019.1.0
...

This might also be a problem with my local FEniCS installation. So let's not worry about this issue here. If the other time stepping schemes work on your side that's good. If we want to test things properly and exclude side effects originating from the runtime system we can use the docker image provided by the FEniCS project and make the time stepping scheme configurable through some command line argument. But this is not really relevant yet.

Import from FEniCS and simplify code.

b1f1639