Consistent treatment of flow and heat sources in models

[jhabriel]

As pointed out in https://github.com/pmgbergen/porepy/discussions/674#discussion-4187333 flow and heat sources must have a consistent treatment in their respective mass and energy equations.

The proposed changes consist of keeping the standard definition of an integrated source term (i.e., as given in the governing equation) and including the scaling with time-step explicitly.

We shall assume that all semi-discrete conservation laws are written such that the time-step $\tau$ is multiplying the whole equation, i.e.:

$(\phi^{n} - \phi^{n-1}) + \tau \nabla\cdot \mathbf{q} - \tau f = 0$

To introduce consistent changes in all current models, the following modifications are necessary:

• Model: SlightlyCompressibleFlow

Proposed change:

# LINE99 from slightly_compressible_flow_model.py
# Compute accumulation term using BE
accumulation_term = (
    accumulation_term.mass * (p - p.previous_timestep())
)

# Retrieve incompressible flow equations  
incomp_eq =  self._eq_manager.equations["subdomain_flow"]

#  Update flow equation
self._eq_manager.equations["subdomain_flow"] = (
    self._ad.time_step * incomp_eq
    + accumulation_term
)

• Model ContactMechanicsBiot

Proposed change:

#LINE821 from contact_mechanics_biot_model.py    
eq = (
    ad.time_step
    * (  # Time scaling of flux terms, both inter-dimensional and from
        # the higher dimension
        div_scalar * flux
        - ad.mortar_projections_scalar.mortar_to_secondary_int
        * ad.interface_flux
    )
    + accumulation_all
    + ad.subdomain_projections_scalar.cell_prolongation([self._nd_subdomain()])
    * biot_accumulation_primary
    - ad.subdomain_projections_scalar.cell_prolongation(g_frac)
    * accumulation_fracs
    - ad.time_step * flow_source  # note explicit scaling by time step
)

• Model THM

Proposed change:

#LINE764 from thm_model.py
eq = (
    ad.time_step
    * (
        # Heat fluxes internal to the subdomain, and from lower-dimensional
        # neighboring subdomains (see self._heat_flux())
        div_scalar * heat_flux
        # Conductive heat flux from higher-dimensional neighbors
        - ad.mortar_projections_scalar.mortar_to_secondary_int
        * ad.conductive_interface_flux
        # Advective heat flux from higher-dimensional neighbors
        - ad.mortar_projections_scalar.mortar_to_secondary_int
        * ad.advective_interface_flux
    )
    + accumulation_all
    + ad.subdomain_projections_scalar.cell_prolongation([self._nd_subdomain()])
    * biot_accumulation_primary
    - ad.time_step * heat_source  # note explicit scaling by time step
)

[keileg]

I added this to the coming model update, it should, if not fixed, then at least clearer what to do there.

[keileg]

@jhabriel: What do you think of this one?

[jhabriel]

Now that we have proper mass (and not volumetric) balances, it should be understood that an external source is the rate of addition of fluid mass.

On that note, I think it is fair to close the issue.