Subcycling in partitioned heat conduction: OpenFOAM

[MakisH]

Follow-up of #281. Needs the latest preCICE develop (precice/precice#1892).

---

(outdated) This is actually exploding at the moment:

---[precice]  iteration: 1 of 100 (min 1), time-window: 2, time: 0.18 of 1, time-window-size: 0.1, max-time-step-size: 0.02, ongoing: yes, time-window-complete: no, 
---[preciceAdapter] The solver's timestep is smaller than the coupling timestep. Subcycling...
Time = 0.19

DICPCG:  Solving for T, Initial residual = 0.00414973, Final residual = 8.40349e-13, No Iterations 78
ExecutionTime = 1.11 s  ClockTime = 2 s

---[precice]  iteration: 1 of 100 (min 1), time-window: 2, time: 0.19 of 1, time-window-size: 0.1, max-time-step-size: 0.01, ongoing: yes, time-window-complete: no, 
---[preciceAdapter] The solver's timestep is smaller than the coupling timestep. Subcycling...
Time = 0.2

DICPCG:  Solving for T, Initial residual = 0.0041563, Final residual = 8.41794e-13, No Iterations 78
ExecutionTime = 1.15 s  ClockTime = 2 s

---[precice]  relative convergence measure: relative two-norm diff of data "Heat-Flux" = 1.02e+00, limit = 1.00e-10, normalization = 2.02e+01, conv = false
---[precice]  relative convergence measure: relative two-norm diff of data "Temperature" = 6.99e-02, limit = 1.00e-10, normalization = 3.48e+01, conv = false
/usr/include/c++/11/bits/stl_vector.h:1063: std::vector<_Tp, _Alloc>::const_reference std::vector<_Tp, _Alloc>::operator[](std::vector<_Tp, _Alloc>::size_type) const [with _Tp = precice::time::Stample; _Alloc = std::allocator<precice::time::Stample>; std::vector<_Tp, _Alloc>::const_reference = const precice::time::Stample&; std::vector<_Tp, _Alloc>::size_type = long unsigned int]: Assertion '__n < this->size()' failed.
[tickly:403477] *** Process received signal ***
[tickly:403477] Signal: Aborted (6)
[tickly:403477] Signal code:  (-6)
[tickly:403477] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7fc719a42520]
[tickly:403477] [ 1] /lib/x86_64-linux-gnu/libc.so.6(pthread_kill+0x12c)[0x7fc719a969fc]
[tickly:403477] [ 2] /lib/x86_64-linux-gnu/libc.so.6(raise+0x16)[0x7fc719a42476]
[tickly:403477] [ 3] /lib/x86_64-linux-gnu/libc.so.6(abort+0xd3)[0x7fc719a287f3]
[tickly:403477] [ 4] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x7e0ca)[0x7fc71587e0ca]
[tickly:403477] [ 5] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x765ab1)[0x7fc715f65ab1]
[tickly:403477] [ 6] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x763221)[0x7fc715f63221]
[tickly:403477] [ 7] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x762afd)[0x7fc715f62afd]
[tickly:403477] [ 8] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x280dce)[0x7fc715a80dce]
[tickly:403477] [ 9] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x718cb)[0x7fc7158718cb]
[tickly:403477] [10] /home/gc/repos/precice/precice/build/libprecice.so.3(+0xa50e3)[0x7fc7158a50e3]
[tickly:403477] [11] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x255e30)[0x7fc715a55e30]
[tickly:403477] [12] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x291705)[0x7fc715a91705]
[tickly:403477] [13] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x24fa27)[0x7fc715a4fa27]
[tickly:403477] [14] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x6c4d87)[0x7fc715ec4d87]
[tickly:403477] [15] /home/gc/repos/precice/precice/build/libprecice.so.3(+0x6b024d)[0x7fc715eb024d]
[tickly:403477] [16] /home/gc/repos/precice/precice/build/libprecice.so.3(_ZN7precice11Participant7advanceEd+0x35)[0x7fc715e821f9]
[tickly:403477] [17] /home/gc/OpenFOAM/gc-v2012/platforms/linux64GccDPInt32Opt/lib/libpreciceAdapterFunctionObject.so(_ZN14preciceAdapter7Adapter7advanceEv+0x1c)[0x7fc7165987ec]
[tickly:403477] [18] /home/gc/OpenFOAM/gc-v2012/platforms/linux64GccDPInt32Opt/lib/libpreciceAdapterFunctionObject.so(_ZN14preciceAdapter7Adapter7executeEv+0xb3)[0x7fc7165b11a3]
[tickly:403477] [19] /home/gc/OpenFOAM/gc-v2012/platforms/linux64GccDPInt32Opt/lib/libpreciceAdapterFunctionObject.so(_ZN4Foam15functionObjects28preciceAdapterFunctionObject7executeEv+0x11)[0x7fc7165dfca1]
[tickly:403477] [20] /usr/lib/openfoam/openfoam2012/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam18functionObjectList7executeEv+0x1db)[0x7fc71a5bc7cb]
[tickly:403477] [21] /usr/lib/openfoam/openfoam2012/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam4Time3runEv+0x105)[0x7fc71a5cc315]
[tickly:403477] [22] /usr/lib/openfoam/openfoam2012/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam4Time4loopEv+0x17)[0x7fc71a5c7817]
[tickly:403477] [23] heatTransfer(+0x114d7)[0x55a9b7a4a4d7]
[tickly:403477] [24] /lib/x86_64-linux-gnu/libc.so.6(+0x29d90)[0x7fc719a29d90]
[tickly:403477] [25] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0x80)[0x7fc719a29e40]
[tickly:403477] [26] heatTransfer(+0x12985)[0x55a9b7a4b985]
[tickly:403477] *** End of error message ***
Aborted (core dumped)

It works without subcycling.

[MakisH]

@fsimonis trace logs:

openfoam-dirichlet.log
openfoam-neumann.log

[MakisH]

With precice/precice#1892, this seems to work for me.

[MakisH]

@davidscn @BenjaminRodenberg this is ready to be reviewed and works with the latest develop of the OpenFOAM adapter and preCICE.

[davidscn]

Works in principle. The error is for both OpenFOAM-OpenFOAM cases (t=0.1 and t=0.01) almost the same though.
with subcycling:
l2 error (sqrt(sum(err^2)/n)): 0.0944029
without subcycling
l2 error (sqrt(sum(err^2)/n)): 0.0944041
Not sure what to expect here in terms of accuracy gain, I could compare with a higher-order (in space) FE code if desired. Technically, the changes work of course.

[BenjaminRodenberg]

Works in principle. The error is for both OpenFOAM-OpenFOAM cases (t=0.1 and t=0.01) almost the same though. with subcycling: l2 error (sqrt(sum(err^2)/n)): 0.0944029 without subcycling l2 error (sqrt(sum(err^2)/n)): 0.0944041 Not sure what to expect here in terms of accuracy gain, I could compare with a higher-order (in space) FE code if desired. Technically, the changes work of course.

Can you swap the time stepping rates for the solvers? (don't mind. Both solvers are using 0.01)

Generally the error looks pretty large to me. For the FEniCS case I get the exact solution. But maybe this also has something to do with the different space discretization and therefore the space discretization error is shadowing the time discretization error in your case for any setup.

[BenjaminRodenberg]

The interesting part here is to compare substeps="true" and substeps="false" w.r.t the error. But since the analytical solution is only linear in time and we are doing linear interpolation in any case, it does not really make a difference. If you switch to a manufactured solution with quadratic or higher-degree time dependence, you should be able to see a difference in the error, if you set substeps="false".

Here, I guess it is already good news that everything works and subcycling does not introduce any error. If you want to compare to the v2.5 implementation of this case: There you usually see some error, if you use subcycling.

[MakisH]

If you want to compare to the v2.5 implementation of this case: There you usually see some error, if you use subcycling.

Previously, subcycling was not even working, as we discovered at some point. I would suggest merging this and creating a separate issue for investigating the errors, as I am afraid this could take us into a very deep dive.

Unless you think that something is terribly wrong and we should investigate before releasing the v3-variant of the OpenFOAM adapter.

[davidscn]

Previously, subcycling was not even working, as we discovered at some point. I would suggest merging this and creating a separate issue for investigating the errors, as I am afraid this could take us into a very deep dive.

No it's ok. I guess the space resolution is not particularly high in our case and I would have to take a deeper look into the error calculation. This is good to go here and I would only open an issue if you have some suspicion about the correctness of the results.

[MakisH]

Do I understand correctly that FEniCS defines a 10x10 mesh? What kind of elements are these?

tutorials/partitioned-heat-conduction/fenics/problem_setup.py

Line 36 in f174aa4

nx = ny = 11

In OpenFOAM we define a 100x100 mesh:

tutorials/partitioned-heat-conduction/openfoam-dirichlet/system/blockMeshDict

Line 25 in f174aa4

hex (0 1 2 3 4 5 6 7) (100 100 1) simpleGrading (1 1 1)

Are we using the same convergence criteria? I see some 1e-14 in FEniCS (not sure what for):

tutorials/partitioned-heat-conduction/fenics/problem_setup.py

Line 19 in f174aa4

tol = 1E-14

In OpenFOAM, the solver uses 1e-6:

tutorials/partitioned-heat-conduction/openfoam-dirichlet/system/fvSolution

Line 15 in f174aa4

tolerance 1e-06;

[BenjaminRodenberg]

The mesh does not matter for this setup, if you do everything correctly and do not introduce any errors (adapter, incorrect BCs, preCICE). I'm using second order Lagrange elements in the FEniCS case. This is not the case because they are absolutely needed for the original problem as a second degree polynomial solution means first degree elements would be fine. We need second degree elements because otherwise the reconstructed heat flux (gradient of temperature) would be incorrect (grad of piecewise linear has jumps). In the corresponding chapter of the FEniCS tutorials book you can also find some more information on the monolithic setup and the method of manufactured solutions (https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf).

@MakisH the tol you are mentioning above is just the tolerance w.r.t checks whether a certain point is inside the geometry or not. FEniCS uses this to determine if mesh points lie on the boundary.

[MakisH]

So, @BenjaminRodenberg @davidscn should we merge this already, or not?

[BenjaminRodenberg]

I would merge this PR. The subcycling implementation looks ok.

I guess there was already an error/inaccuracy (compared to the FEniCS case) in how you are treating the boundary conditions before. I would suggest to open an issue and try to fix this to reach the same accuracy as for the FEniCS case without subcycling in a separate PR. Not sure whether this is actually possible with OpenFOAM: In FEniCS, I needed 2nd order FEM to do the trick with the flux. I have no idea about the capabilities of OpenFOAM and FV here.