[antonin.m]

Hi everyone,

I'm currently working on a hydrological model that aims to unify Richards' equations for sub-surface infiltration flows and surface runoff. We've decided to model this runoff using a single surface mesh layer, introducing a “flag” variable that tells us whether we're in the subsurface or in the runoff layer (value of -1 in runoff, +1 in subsurface). By analogy for the different variables in the system, a single linear system must solve the runoff and subsurface equations at each iteration. The linear system can be written as follows:

{
  fvScalarMatrix psiEqn
  (
      Crel*fvm::ddt(psi)
  ==
      (
          Krel.is_isotropic()
        ? fvm::laplacian(Krel.isotropic(),   psi, "laplacian(Krel,psi)")
        : fvm::laplacian(Krel.anisotropic(), psi, "laplacian(KrelT,psi)")
      )
    +
      (
          Krel.is_isotropic()
        ? (fvc::grad(Krel.isotropic(), "grad(Krel)") & vuz)
        : (fvc::div((Krel.anisotropic() & vuz), "div(Krel&z)"))
      )
    - AET
  );

  psiEqn.solve();
}

The solution was therefore to do the case disjunction on flag for Crel, Krel, AET, here's an example for the calculation of Krel with Manning's equations that we wish to use at the surface:

template<class HydConductivityFieldType, class VanGenuchtenFieldType, class ManningFieldType, class FlagFieldType, class IdFieldType, class KadimFieldType>
tmp<HydConductivityFieldType> calc_Krel_Impl
(
    const HydConductivityFieldType& K,
    const volScalarField& psi,
    const volScalarField& H,
    const VanGenuchtenFieldType& n,
    const ManningFieldType& m,
    const volScalarField& thtil,
    const FlagFieldType& flag,
    const IdFieldType& Id,
    const KadimFieldType& Kadim,
    const fvMesh& mesh
)
{
    return
    (
        K * pos(flag) * // subsurface
            (
                pos0(psi)
                + neg(psi)
                * sqrt(thtil)
                * sqr(1 - pow((1-pow(thtil,(n/(n - 1)))), (1 - (1/n))))
            )

        + neg(flag) * // run-off layer
            ( Kadim * // this part calculates Kxx and Kyy components in the run-off tensor.
                (
                    pos0(psi)
                    * pow(psi, 5/3)/(m*sqrt(0.0016)) / faceArea
                )
            + Id // ... while this part keeps the Kzz equal to 1.
            )
    );
}

In order to realize our analogy between our two equations, one of which is expressed in [L.T-1] and [T-1], for the surface runoff layer I need to divide the quantities Crel, Krel and AET by the cell area parallel to the vertical axis (the thickness being 1, this actually corresponds to the cell height).

I tried several implementations based on what I read online, such as

point p(0.,0.,0.); //closest point to the run-off downslope cell
label cellNo = mesh.findCell(p);
scalar faceArea = 0.0;
const cell& cFace = mesh.cells()[cellNo];
faceArea = mag(mesh.Sf()[cFace[0]]);

but this one and other variations, for example, can't be used for parallel calculations, creating Segmentation Faults, which I understand because the mesh is divided and not all processors have access to the cell. This wasn't a problem when I implemented the boundary condition myself, as I was looping through the cells on the patch, but here, looping through all the cells again in each iteration would be a real problem for performance.
I realize that all methods based on mesh.magSf() combined with a particular cell search don't work in parallel for this reason. By using the mesh.magSf() surface field directly in the calculations (Crel / mesh.magSf()), I'm having operator problems between the different types.
I haven't found any other way around this problem without directly writing Crel = Crel / 5. in the code, because I know its precise height but I'd like to avoid having to recompile my module every time I modify the case study. By the way, the 0.0016 value corresponds to the slope at the surface and I encountered the same issue weeks ago, worked on it and gave up but I'm not willing to reproduce this and that's why I finally post this.
I'm therefore open to any hint, even the obvious ones that I might have missed, as I'm more or less new to OpenFoam.

Many thanks in advance. Best regards,
Antonin.