[Muerio]

Hello everyone,

I'm trying to update the values of a custom field I added to the multicomponentFluid solver in OpenFOAM 11.
The field should only depend on the values from a fvModel and shouldn't be affected by the flow itself.

Therefore I tried to extend the thermophysicalPredictor like this (copied from component fields Y_ and removed terms for transport):

void Foam::solvers::electrochemFluid::thermophysicalPredictor()
{
    tmp<fv::convectionScheme<scalar>> mvConvection(
        fv::convectionScheme<scalar>::New(mesh, fields, phi,
                                          mesh.schemes().div("div(phi,Yi_h)")));

    reaction->correct();

    volScalarField &myField = myField_;
    fvScalarMatrix myFieldEqn(fvm::ddt(myField) == fvModels().source(myField));


    myFieldEqn.relax();

    fvConstraints().constrain(myFieldEqn);

    myFieldEqn.solve("myField");

    fvConstraints().constrain(myField);

    forAll(Y, i)
    {
        ...

For testing purposes I use a simple fvModel that assigns the value of 1.0 to all cells:

forAll(cells, i){eqn.source()[cells[i]] = 1.0;}

All boundaries are defined as zeroGradient:

internalField   uniform 0;

boundaryField
{
    ".*"
    {
        type            zeroGradient;
    }
}

It compiles and starts without issue, but once it reaches the "myFieldEqn.solve" line, the solver crashes with an Arithmetic exception in Foam::divide. I guess it divides by 0 while trying to solve?!

Is there something obvious I missed?
What could I try to resolve this issue?

Thank you in advance for any help!