Dear Mona,
The discretization error depends of course on the discretization scheme
and discretization length (for the ODE just Δt). For DuMux currently
only backwards Euler (implicit, default in DuMux) and forwards Euler
(explicit) are readily available. The only real control you have is
through the time step size, by choosing a lower maximum time step size.
In your example the rates are quite high, requiring small time steps for
good accuracy.
Note that with numeric vs analytic differentiation we refer to how the
discretized equations are solved. For a linear ODE the two produce the
same result and the discretized equation is solved exactly.
Best regards,
Mathis
On 11/08/2023 12:50, Giraud, Mona wrote:
Dear Mathis,
Many thanks for this speedy reply!
If I understood correctly, there are no other ways of getting more
accurate results (with numeric differentiation ) than by manually
setting a lower maximal time step?
As you correctly guessed, I was indeed using "MaxTimeStepSize = 1 # [s]".
a) no indeed, I am (currently) using the numeric differentiation with
(the default) forward differences method.
With best regards,
Mona
------------------------------------------------------------------------
*From:* DuMux <dumux-boun...@listserv.uni-stuttgart.de> on behalf of
Mathis Kelm <mathis.k...@iws.uni-stuttgart.de>
*Sent:* Friday, August 11, 2023 12:14:32 PM
*To:* dumux@listserv.uni-stuttgart.de
*Subject:* Re: [DuMux] implementation of solDependentPointSource
Dear Mona Giraud,
if you want a global source term you should be able to simply use the
`source(..)` function [1] and make use of the elemVolVars and scv
parameters just as you do now in `bioChemicalReactions(..)`, see e.g. [2].
Regarding the result you obtain, what time step are you using?
Discretizing your equation dS/dt = r S with an implicit Euler scheme
yields C^(n+1) - C^(n) = Δt r C^(n+1) which rearranges to C^(n+1)(1 -
Δt r) = C^(n) and C^(n+1) = C^(n) / (1 - Δt r) which for Δt = 1s would
give the observed output.
Regarding other questions:
a) The `add[...]Derivatives` functions are only needed for analytic
differentiation which I assume you are not using?
c) I don't think you need to overload any other functions. As far as I
can see analytic derivatives of point sources would fall under
addSourceDerivatives too, not that you need them.
d) As I understand it the `addPointSources` function lets you register
source functions at different points, for `SolDependentPointSource`
passing a `ValueFunction` that is evaluated in its `update` function.
In `scvPointSources` first this update is called and then the
problem's `pointSource()` which again receives all the local
information to compute a solution dependent source. You always need to
add your point sources but can choose whether to pass a value function
for the update call or to overload `pointSource` in your problem.
Best regards,
Mathis
[1]:
https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/releases/3.0/dumux/common/fvproblem.hh#L294
[2]:
https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/releases/3.0/test/porousmediumflow/2pnc/implicit/fuelcell/problem.hh#L181
On 11/08/2023 11:07, Giraud, Mona wrote:
Dear members of the DuMux community,
I am contacting you because I would like to get some pointers on how
to implement a source term [1] valid in the whole domain, [2] which
would depend on the solution vector.
From what I understood, I need to use the "PointSource" function
rather than the source() function.
After looking at "dumux\test\porousmediumflow\richardsnc\implicit", I
added the following:
*1) in the main.cc:*
before "timeLoop->start();":
problem->computePointSourceMap()
*2) in the problem.hh: *
using PointSource = GetPropType<TypeTag, Properties::PointSource>;
using Problem = GetPropType<TypeTag, Properties::Problem>;
[...]
void addPointSources(std::vector<PointSource>& pointSources) const
{
// loop over all vertices of the 1D domain to get one source
per subcontrol volume.
for (const auto& vertex : Dune::vertices(
this->fvGridGeometry().gridView() ) )
{
auto globalPos = vertex.geometry().center();//1D space
pointSources.emplace_back(globalPos,
[this](const Problem &problem,
const Element &element,
const FVElementGeometry &fvGeometry,
const ElementVolumeVariables &elemVolVars,
const SubControlVolume &scv)
{
return bioChemicalReaction(elemVolVars,scv);
}
);
}
}
PrimaryVariables biochemicalReactions(const ElementVolumeVariables
&elemVolVars,
const SubControlVolume &scv)
{
const auto& volVars = elemVolVars[scv];
double volume = scv.volume(); //to have m3 scv.
const auto massOrMoleDensity = [](const auto& volVars, const int
compIdx, const bool isFluid)
{
return (useMoles ? volVars.molarDensity(compIdx) :
volVars.density(compIdx) );
};
const auto massOrMoleFraction= [](const auto& volVars, const int
phaseIdx, const int compIdx)
{
return (useMoles ? volVars.moleFraction(phaseIdx, compIdx) :
volVars.massFraction(phaseIdx, compIdx) )
};
//mol comp/m3 scv
double C_S = massOrMoleDensity(volVars, h2OIdx) //mol phase/ m3 phase
* massOrMoleFraction(volVars,h2OIdx, mucilIdx) //mol comp/ mol phase
* volVars.saturation(h2OIdx) //m3 phase / m3 pores
* volVars.porosity(); // m3 pores/m3 scv
double rate = - 0.5; //1/s
double dwater_dt = 0.;
// reaction rate in mol/s for all components
// of the subcontrol volume
double dC_S_dt = rate * C_S * volume; // example of a simple reaction
return PrimaryVariables({dwater_dt,dC_S_dt});
}
*3) in properties.hh:*
template<class TypeTag>
struct PointSource<TypeTag, TTag::RichardsNCTT> { using type =
SolDependentPointSource<TypeTag>; };
####
For this simplified toy model, [1] there are no flux BCs, [2] the
(molar) density of the phase is always equal to that of pure water,
[3] there is a simple equation for the reaction to make the
comparison with a reference analytical solution easier,
(the biochemical reactions implemented later will be more complex).
Therefore the only change in the concentration of the component comes
from the reaction = [d C_S/dt = -0.5 * C_S].
I have the issue that the solver yields at each time step:
[ C_S_end = C_S_init/(1 - rate ) ]
a) Do I need to overload the "addSourceDerivatives()" function and/or
change the solver settings to get the correct result?
b) If yes, do you have one or several examples using the
"addSourceDerivatives()"?
c) Are there other functions which need to be overloaded/changed
(like "pointSourceDerivative()")?
d) If I understood correctly, using "addPointSources()" or
"pointSource()" yield the same results? cf. comments in the
function "dumux\common\fvproblem.hh:scvPointSources()"
NB: I use version 3.0 of DuMux
Many thanks for your help!
With best regards,
Mona Giraud
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
Forschungszentrum Juelich GmbH
52425 Juelich
Sitz der Gesellschaft: Juelich
Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498
Vorsitzender des Aufsichtsrats: MinDir Stefan Müller
Geschaeftsfuehrung: Prof. Dr. Astrid Lambrecht (Vorsitzende),
Karsten Beneke (stellv. Vorsitzender), Dr. Ir. Pieter Jansens, Prof.
Dr. Frauke Melchior
------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
--
********************************************************
Mathis Kelm
Universitaet Stuttgart
Institut fuer Wasser- und Umweltsystemmodellierung
Lehrstuhl fuer Hydromechanik und Hydrosystemmodellierung
Pfaffenwaldring 61, 70569 Stuttgart
Tel.: 0711/685-60146
mathis.k...@iws.uni-stuttgart.de
https://www.iws.uni-stuttgart.de/lh2/
********************************************************
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
--
********************************************************
Mathis Kelm
Universitaet Stuttgart
Institut fuer Wasser- und Umweltsystemmodellierung
Lehrstuhl fuer Hydromechanik und Hydrosystemmodellierung
Pfaffenwaldring 61, 70569 Stuttgart
Tel.: 0711/685-60146
mathis.k...@iws.uni-stuttgart.de
https://www.iws.uni-stuttgart.de/lh2/
********************************************************
_______________________________________________
DuMux mailing list
DuMux@listserv.uni-stuttgart.de
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
