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<mailto: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<mailto: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
