Yan and I noticed that all FV terms are treated incorrectly in the Generic FV classes or not present at all. We are fixing these step by step and checking the revised files into the SVN. We strongly suggest that the Leoben team has a look at this too, since otherwise the same error will appear in CSM++ 1.0. I am not sure how the ETH team is using source/sink terms for their pressure-enthalpy simulations, but this could have an effect too? Here are some thoughts:

Take a standard linear advection equation with a conserved variable

(1)where $\phi$ is the conserved variable $\mathbf{u}$ is the velocity, and $q$ is the source term. Standard finite volume discretisation yields on a control volume $V_i$ over time increment $\Delta t$

(2)where $\mathcal{F}$ is the usual upwind weighted finite volume discretisation of the flux term (in implicit or explicit manner). Now, obviously, this becomes

(3)If you have a look at the function `ComposeSolution()` in `ExplicitNodeCenteredTransport` you will see the following

for ( stl_index nidx=0U; nidx<net.FiniteVolumes(); nidx++ ) { // reading the transported variable value fT solution = net.N(nidx)->Read( pmem, NodeCenteredFiniteVolumeTransport<fT,dim>::ad1_key.index ); // adding potential (volumetric) source or sink terms due to a divergence of the flow (+ dt sum_j^e 1/3 V_e q_j) if ( with_flux_balance_correction ) { if ( this->FLUX_BALANCE[nidx] != static_cast<fT>(0.) ) RESULT[nidx] += solution * -this->FLUX_BALANCE[nidx]; } // adding externally assigned source or sink terms RESULT[nidx] += solution * net.N(nidx)->Read( pmem, this->src_key.index ); // subtracting the flux time-interval product RESULT[nidx] = solution - (time_interval / this->FVPOREVOL[nidx]) * RESULT[nidx]; }

Here, `nidx` corresponds to the ID of a control volume, `RESULT` to $\mathcal{F}$, `solution` to $\phi^t$, `time_interval` to $\Delta t$, and `net.N(nidx)->Read( pmem, this->src_key.index );` to $q$. This implies that we are solving

This is obviously wrong for several reasons (source has the wrong sign, is multiplied with $\phi^t$ and divided by the control volume $V_i$.

The function `Compose2PhaseSolution()` in `TwoPhaseExplicitNodeCenteredFVTransport` does even worse:

// 2. Nodal fluid SOURCE/SINK terms const fT nodal_src = net.N(nidx)->Read( pmem, this->src_key.index ); // if there is a source term if ( nodal_src != static_cast<fT>(0.) ) { // alarm if source term exceeds second stability criterion, Guignot, p. 166 (first criterion does // not apply since the source term should be positive in certain situations) if ( (nodal_src / RESULT[nidx]) * time_interval < -2. ) { #ifndef NDEBUG cerr <<"\nTwoPhaseExplicitNodeCenteredFVTransport<fT,dim>::Compose2PhaseSolution: "; cerr <<" cell "<< nidx <<", source term: "<< nodal_src <<" violated stability criterion."; #endif // explicit solution for source term at delta t_max (Guignot, p. 166) RESULT[nidx] += nodal_src * (-2. * RESULT[nidx] / nodal_src); } // overshoot criterion else if ( (nodal_src * time_interval + RESULT[nidx]) > 1. ) { #ifndef NDEBUG cerr <<"\nTwoPhaseExplicitNodeCenteredFVTransport<"<< typeid(fT).name() <<","<< dim <<"m>::Compose2PhaseSolution: "; cerr <<" cell "<< nidx <<", source term: "<< nodal_src <<" violated stability criterion."; #endif // explicit solution for source term at delta t_max (Guignot, p. 166) RESULT[nidx] += nodal_src * ((1. - RESULT[nidx]) / nodal_src); } // explicit solution for source terms due to capillary flow O.K. else RESULT[nidx] += nodal_src * time_interval; } // 3. ACCUMULATION: subtracting the flux time-interval product RESULT[nidx] = net.N(nidx)->Read( pmem, this->ad1_key.index ) - (time_interval / this->FVPOREVOL[nidx]) * RESULT[nidx];

Here, we end up with the following equations, depending on which statement is reached in the `if else` condition:

or

(6)or

(7)Again, the source term has the wrong sign, is multiplied by $\Delta t/V_i$ or $\Delta t^2/V_i$ or does not appear at all (instead, $\mathcal{F}$ is used).

As far as we can see, source terms are not at all considered in the implicity finite volume schemes. A variable `"nodal source"` is read in but never used. As mentioned above, we will work on these to fix it and provide new files via the SVN. Have a look at our recent plea of making more use of the SVN