Mesh orientation effect for the temperature propagation

[i-l-s-h-a-t]

Describe the bug
The temperature propagation is mesh-oriented for a hexahedral inclined mesh (tstep 1):

There are 3 meshes used: hexahedral rectangular, hexahedral inclined, and tetrahedral.
Neither of them work for the purpose:

• hexahedral inclined - has the mesh-orientation effect
• hexahedral rectangular - is not aligned to materials boundary, and it was made just to check the mesh orientation effect
• tetrahedral - doesn't converge

Results for the tetrahedral mesh (only timesteps before the simulation divergence are available):

Results for the hexahedral rectangular mesh:

To Reproduce
Steps to reproduce the behavior:

1. Run the models (the link is provided below)
2. Load the resulting vtk files into Paraview and observe the issue

Expected behavior

1. The temperature propagation front is perpendicular to Z-axis

Platform (please complete the following information):
Machine: RHEL 8.10
Compiler: gcc 13.1.0
GEOS Version: 1.1.0
GEOS Version sha: develop 1e617be , tpls: 9687680c61fa3378f7b2ce7eea8327736809fb63

Models and logs
You can find them in this repository: https://gitlab.com/i-l-s-h-a-t/geos_models/-/tree/main/THM/triaxial/thermal_issue/mesh_effect?ref_type=heads

Additional context
The idea of the model is to simulate a triaxial experiment for the heterogenous rock.
Model uses 2 materials (without contact mechanics), currently they have the same properties. It uses singlePhase singleComponent physics with TPFA discretization.

Next boundary conditions are defined: no load from the top and side boundaries, lower boundary has a roller condition except the central part where it is fixed. Cooling is imposed at the top side by setting the temperature lower than the initial.

I thought that HybridFVM might help and found there are number of features not implemented for the case:

1. There is no iterative linear solver for the SinglePhaseHybridFVM:
   PoroelasticitySolver: MGR strategy is not implemented for thermal SinglePhasePoromechanics/SinglePhaseHybridFVM
2. After switching to the direct linear solver: The thermal option is not supported by SinglePhaseHybridFVM

Checking the code for CompositionalMultiphaseHybridFVM showed there is no such an exception so maybe adding a dummy component and phase is a solution, but I have not checked yet. Anyway, it would be good to resolve the issue for the SinglePhase physics.

This issue could be related to this #3739.

[paveltomin]

what do you expect for TPFA???

[i-l-s-h-a-t]

what do you expect for TPFA???

I know that TPFA is not good for anisotropic cases. I didn't experiment with TPFA for distorted meshes with a thermal option, but I have a feeling that for homogeneous isotropic media the mesh orientation effect should not be that much pronounced...

Here is an illustration with a temperature cutoff:

[OmarDuran]

While anisotropy is a key factor, non-K-orthogonality also arises from the misalignment of canonical axes composing any tensor (including isotropic) with specific regions of the mesh. This effect is well-documented; for instance, see the reference 10.1063/5.0282130 and the citations therein. Specifically, Figure 20 illustrates how mesh distortion follows a linear transformation, similar to the behavior observed in your mesh.