Hi Edward,
Thank you for you answer, it did help. The code is now compiling (with a
warning-> porousmediumsubproblem.hh:324:19: warning: variable ‘values’
set but not used [-Wunused-but-set-variable]) but doesn't run :-(!
The linear solver doesn't converge so there is still a problem somewhere.
To keep it simple, I just added a constant temperature in the
initialAtPos function in the porousmediumsubproblem and a
freeflowsubproblem but it doesn't work. Any idea why it doesn't converge ?
I am attaching the .hh files that I am using that are a slightly
modified version of the exercises/exercise-coupling-ff-pm/models.
Thank you for your help,
Suzon
On 02/04/2024 16:25, Coltman, Edward wrote:
Hello Suzon,
I'm not sure if this is really where you are working, but assuming:
- the error you are working on is coming from the
"evaluateInterfaceFluxes" function,
- in porousmediumsubproblem.hh header file,
- and you have added a call to evaluate the face-wise flux at
somepoint that looks like this:
>
NumEqVector flux =
couplingManager().couplingData().massCouplingCondition(element,
fvGeometry, elemVolVars, scvf)
* scvf.area() *
elemVolVars[scvf.insideScvIdx()].extrusionFactor() *
FluidSystem::molarMass(1) * -1.0 * 86400.0;
> (taken from the exercise solution)
Then there will be a compatibility problem between the call to the
"massCouplingCondition(...)" function, and the type NumEqVector.
NumEqVector will have three entries (FieldVector<[...],3>, including
the added temperature), and the returned vector from the
massCouplingCondition(...) function will only have two entries
(FieldVector<[...],2>, only evaluating the mass fluxes).
To fix this, you can use the correct type (for this I would recommend
to simply use the auto type).
If this is not where your error is coming from, it would be very
helpful if you could send us a bit more output from the compiler
error. This is just a best guess based on the line number and error
you have mentioned.
Best,
Ned
------------------------------------------------------------------------
*From:* DuMux <[email protected]> on behalf of
Suzon Jammes <[email protected]>
*Sent:* Tuesday, April 2, 2024 4:01 PM
*To:* DuMux User Mailing List
*Subject:* Re: [DuMux] Non isothermal coupled exercise
Hi Christoph,
I was able to use non isothermal condition with another example using
TwoPNCNI model so I guess that the problem comes from the coupling part.
I assume that indeed the temperature should be taken into at the
coupling interface but I am not sure how to do that.
Should I modify something in the couplingManager ?
Thanks for your help,
—————————————————————————————
Suzon Jammes
Researcher & Consultant
M&U sas, Geology by Research
www.mandu-geology.fr
[email protected] / +33(0)768692934
Ce message est confidentiel. Son contenu ne représente en aucun cas
un engagement de la part de M&U SAS sous reserve de tout accord conclu
par écrit avec la société. Toute publication, utilisation
ou diffusion, meme partielle, doit être autorisée préalablement.
Si vous n'êtes pas destinataire de ce message, merci d'en
avertir immédiatement l'expéditeur ou par email :[email protected].
This message is confidential. Its contents do not constitute
a commitment by M&U SAS except
where provided for in a written agreement with the company.
Any unauthorised disclosure, use
or dissemination, either whole or partial, is prohibited. If you
are not the intended recipient, please notify the sender
immediately or with the following email: [email protected]
On 30 Mar 2024, at 16:23, Christoph Grüninger <[email protected]> wrote:
Hi Suzon,
I did not look into the code, so I might guess wrong.
I suspect that by adding temperature, you got a third primary
variable. The error points out that you are using a vector with two
entries (representing two primary variables) at a place where you
should use a vector with three entries (for all three primary vars).
I hope this helps you figuring out what went wrong.
Bye
Christoph
Am 29.03.24 um 11:14 schrieb Suzon Jammes:
Hi all,
I am now working from the coupled exercise
(exercises/exercise-coupling-ff-pm/models) with a 2pnc model. I
wanted to make it non-isothermal so I change TwoPNC by TwoPNCNI and
NavierStockesNC by NavierStockesNCNI in the properties file but I
get the following error :
error: conversion from ‘FieldVector<[...],2>’ to non-scalar type
‘FieldVector<[...],3>’ requested
263 | elemVolVars[scvf.insideScvIdx()].extrusionFactor() *
FluidSystem::molarMass(1) * -1.0 * 86400.0;
I don't understand what does it mean. Does anybody know how to fix
this problem ?
Thank you
Suzon Jammes
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
--
L'enjeu est de bâtir la France de nos enfants,
pas de ressasser la France de notre enfance.
[Emmanuel Macron, 2022]
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \ingroup NavierStokesTests
* \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
*/
#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/freeflow/navierstokes/staggered/problem.hh>
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
namespace Dumux {
/*!
* \ingroup NavierStokesTests
* \brief Test problem for the one-phase compositional (Navier-) Stokes problem.
*
* Horizontal flow from left to right with a parabolic velocity profile.
*/
template <class TypeTag>
class FreeFlowSubProblem : public NavierStokesStaggeredProblem<TypeTag>
{
using ParentType = NavierStokesStaggeredProblem<TypeTag>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
using FluidState = GetPropType<TypeTag, Properties::FluidState>;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
public:
FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Freeflow"),
eps_(1e-6),
couplingManager_(couplingManager)
{
velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
temperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "SpatialParams.Temperature");
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary segment.
*
* \param element The finite element
* \param scvf The sub control volume face
*/
BoundaryTypes boundaryTypes(const Element& element,
const SubControlVolumeFace& scvf) const
{
BoundaryTypes values;
const auto& globalPos = scvf.center();
if(onLeftBoundary_(globalPos))
{
values.setDirichlet(Indices::conti0EqIdx + 1);
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
else if(onRightBoundary_(globalPos))
{
values.setDirichlet(Indices::pressureIdx);
values.setOutflow(Indices::conti0EqIdx + 1);
}
else
{
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
values.setNeumann(Indices::conti0EqIdx);
values.setNeumann(Indices::conti0EqIdx + 1);
}
if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
values.setCouplingNeumann(Indices::conti0EqIdx);
values.setCouplingNeumann(Indices::conti0EqIdx + 1);
values.setCouplingNeumann(Indices::momentumYBalanceIdx);
values.setBeaversJoseph(Indices::momentumXBalanceIdx);
}
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Dirichlet control volume.
*
* \param element The element
* \param scvf The subcontrolvolume face
*/
PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
{
PrimaryVariables values(0.0);
values = initialAtPos(pos);
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables
* \param elemFaceVars The element face variables
* \param scvf The boundary sub control volume face
*/
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFaceVariables& elemFaceVars,
const SubControlVolumeFace& scvf) const
{
PrimaryVariables values(0.0);
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
values[Indices::conti0EqIdx] = massFlux[0];
values[Indices::conti0EqIdx + 1] = massFlux[1];
}
return values;
}
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param globalPos The global position
*/
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
// This is only an approximation of the actual hydorostatic pressure gradient.
// Air is compressible (the density depends on pressure), thus the actual
// vertical pressure profile is non-linear.
// This discrepancy can lead to spurious flows at the outlet if the inlet
// velocity is chosen too small.
FluidState fluidState;
updateFluidStateForBC_(fluidState, pressure_);
const Scalar density = FluidSystem::density(fluidState, 0);
values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->gridGeometry().bBoxMax()[1] - globalPos[1]);
values[Indices::temperatureIdx] = temperature_;
// gravity has negative sign
values[Indices::conti0EqIdx + 1] = moleFraction_;
values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->gridGeometry().bBoxMin()[1])
* (this->gridGeometry().bBoxMax()[1] - globalPos[1])
/ (height_() * height_());
return values;
}
/*!
* \brief Return the sources within the domain.
*
* \param globalPos The global position
*/
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }
/*!
* \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
*/
Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
{ return couplingManager().couplingData().darcyPermeability(element, scvf); }
/*!
* \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
*/
Scalar alphaBJ(const SubControlVolumeFace& scvf) const
{ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); }
/*!
* \brief Set the coupling manager
*/
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManager_ = cm; }
/*!
* \brief Get the coupling manager
*/
const CouplingManager& couplingManager() const
{ return *couplingManager_; }
void setTimeLoop(TimeLoopPtr timeLoop)
{ timeLoop_ = timeLoop; }
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
bool onRightBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
//! \brief updates the fluid state to obtain required quantities for IC/BC
void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
{
fluidState.setTemperature(this->spatialParams().temperatureAtPos({}));
fluidState.setPressure(0, pressure);
fluidState.setSaturation(0, 1.0);
fluidState.setMoleFraction(0, 1, moleFraction_);
fluidState.setMoleFraction(0, 0, 1.0 - moleFraction_);
typename FluidSystem::ParameterCache paramCache;
paramCache.updatePhase(fluidState, 0);
const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
fluidState.setDensity(0, density);
const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
fluidState.setMolarDensity(0, molarDensity);
const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
fluidState.setEnthalpy(0, enthalpy);
}
// the height of the free-flow domain
const Scalar height_() const
{ return this->gridGeometry().bBoxMax()[1] - this->gridGeometry().bBoxMin()[1]; }
Scalar eps_;
Scalar velocity_;
Scalar pressure_;
Scalar moleFraction_;
Scalar temperature_;
TimeLoopPtr timeLoop_;
std::shared_ptr<CouplingManager> couplingManager_;
};
} //end namespace Dumux
#endif
[TimeLoop]
DtInitial = 1 # s
EpisodeLength = -360 # s # 0.25 days
TEnd = 256000 # s # 2 days
[Freeflow.Grid]
LowerLeft = 0 1
UpperRight = 1 2
Cells = 16 16
[PorousMedium.Grid]
UpperRight = 1 1
Cells = 16 16
[Freeflow.Problem]
Name = freeflow
EnableGravity = true
EnableInertiaTerms = true
MoleFraction = 0.0 # -
Pressure = 1e5 # Pa
Velocity = 1 # m/s
[PorousMedium.Problem]
Name = porousmedium
EnableGravity = true
Saturation = 0.1 # -
MoleFraction = 0.1 # -
Pressure = 1e5 # Pa
[SpatialParams]
Permeability = 2.65e-10 # m^2
Porosity = 0.4 # -
AlphaBeaversJoseph = 1.0 # -
# EXNUMBER >= 1
VanGenuchtenN = 8.0
VanGenuchtenAlpha = 6.5e-4
Swr = 0.005
Snr = 0.01
PcLowSw = Swr * 5.0
PcHighSw = 1.0 - Snr * 5.0
Temperature = 293.15
[Problem]
Name = models_coupling
ExportStorage = true
PlotStorage = true
ExportFluxes = true
[Newton]
MaxSteps = 12
MaxRelativeShift = 1e-7
TargetSteps = 5
[Vtk]
AddVelocity = 1
[Assembly]
NumericDifference.BaseEpsilon = 1e-8
[Component]
SolidDensity = 2700 # solid density of granite
SolidThermalConductivity = 2.8 # solid thermal conducitivity of granite
SolidHeatCapacity = 790 # solid heat capacity of granite
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief A simple Darcy test problem (cell-centered finite volume method).
*/
#ifndef DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
#define DUMUX_POROUSMEDIUMFLOW_SUBPROBLEM_HH
#include <dumux/common/properties.hh>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/io/gnuplotinterface.hh>
#include <dumux/common/metadata.hh>
#include <dumux/porousmediumflow/problem.hh>
namespace Dumux {
/*!
* \brief The porous medium flow sub problem
*/
template <class TypeTag>
class PorousMediumSubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
// copy some indices for convenience
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
// grid and world dimension
static constexpr int dim = GridView::dimension;
static constexpr int dimworld = GridView::dimensionworld;
// primary variable indices
static constexpr int conti0EqIdx = Indices::conti0EqIdx;
static constexpr int pressureIdx = Indices::pressureIdx;
static constexpr int phaseIdx = 0;
static constexpr int temperatureIdx = Indices::temperatureIdx;
// TODO: dumux-course-task 2.A
// set the `transportCompIdx` to `Indices::switchIdx`.
static constexpr int transportCompIdx = Indices::switchIdx;
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
public:
PorousMediumSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "PorousMedium"),
eps_(1e-7),
couplingManager_(couplingManager)
{
moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
temperature_= getParamFromGroup<Scalar>(this->paramGroup(), "SpatialParams.Temperature");
// initialize storage output file (.csv)
exportStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.ExportStorage", false);
plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
if (exportStorage_)
{
storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
initializeStorageOutput();
}
// initialize flux output file (.json)
exportFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.ExportFluxes", false);
if (exportFluxes_)
{
fluxFileName_ = "flux" + getParam<std::string>("Problem.Name");
simulationKey_ = getParam<std::string>("Problem.Name", "case1");
initializeFluxOutput();
}
}
void initializeStorageOutput()
{
storageFile_.open(storageFileName_);
storageFile_ << "#Time[s]" << ";"
<< "WaterMass[kg]" << ";"
<< "WaterMassLoss[kg]" << ";"
<< "EvaporationRate[mm/d]"
<< std::endl;
}
void initializeFluxOutput()
{
Dumux::MetaData::Collector collector;
if (Dumux::MetaData::jsonFileExists(fluxFileName_))
Dumux::MetaData::readJsonFile(collector, fluxFileName_);
Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
}
template<class SolutionVector, class GridVariables>
void startStorageEvaluation(const SolutionVector& curSol,
const GridVariables& gridVariables)
{
// TODO: dumux-course-task 2.B
// Initialize `initialWaterContent_` and assign that to `lastWaterMass_`.
initialWaterContent_ = evaluateWaterMassStorageTerm(curSol, gridVariables);
lastWaterMass_ = initialWaterContent_;
}
void startFluxEvaluation()
{
// Get Interface Data
std::vector<GlobalPosition> faceLocation;
std::vector<Scalar> faceArea;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
for (auto&& scvf : scvfs(fvGeometry))
{
if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
continue;
faceLocation.push_back(scvf.ipGlobal());
faceArea.push_back(scvf.area());
}
}
std::vector<Scalar> faceEvaporation(faceArea.size(), 0.0);
Dumux::MetaData::Collector collector;
if (Dumux::MetaData::jsonFileExists(fluxFileName_))
Dumux::MetaData::readJsonFile(collector, fluxFileName_);
collector[simulationKey_]["Time"]= {0.0};
collector[simulationKey_]["TimeStepIdx"]= {0};
collector[simulationKey_]["FaceCenterCoordinates"] = faceLocation;
collector[simulationKey_]["FaceLength"] = faceArea;
collector[simulationKey_]["FaceEvaporation"] = {faceEvaporation};
collector[simulationKey_]["FluxOverFullSurface"] = {0.0};
Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
}
template<class SolutionVector, class GridVariables>
void postTimeStep(const SolutionVector& curSol,
const GridVariables& gridVariables)
{
evaluateWaterMassStorageTerm(curSol, gridVariables);
evaluateInterfaceFluxes(curSol, gridVariables);
if (exportStorage_ && plotStorage_)
{
gnuplotStorage_.resetPlot();
gnuplotStorage_.setDatafileSeparator(';');
gnuplotStorage_.setXlabel("time [d]");
gnuplotStorage_.setXRange(0.0, getParam<Scalar>("TimeLoop.TEnd"));
gnuplotStorage_.setYlabel("evaporation rate [mm/d]");
gnuplotStorage_.setOption("set yrange [0.0:]");
gnuplotStorage_.setOption("set y2label 'cumulative mass loss'");
gnuplotStorage_.setOption("set y2range [0.0:0.5]");
gnuplotStorage_.setOption("set y2range [0.0:0.5]");
gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
gnuplotStorage_.plot("StorageOverTime");
}
}
template<class SolutionVector, class GridVariables>
Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol,
const GridVariables& gridVariables)
{
// compute the mass in the entire domain
Scalar waterMass = 0.0;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, curSol);
for (auto&& scv : scvs(fvGeometry))
{
// const auto& volVars = elemVolVars[scv];
// TODO: dumux-course-task 2.B
// Insert calculation of the water mass here
const auto& volVars = elemVolVars[scv];
waterMass += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx) * volVars.density(phaseIdx)
* volVars.saturation(phaseIdx) * volVars.porosity()
* scv.volume() * volVars.extrusionFactor();
}
}
Scalar cumMassLoss = initialWaterContent_ - waterMass;
Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400
/ (this->gridGeometry().bBoxMax()[0] - this->gridGeometry().bBoxMin()[0])
/ timeLoop_->timeStepSize();
lastWaterMass_ = waterMass;
if (exportStorage_)
{
storageFile_ << timeLoop_->time() << ";"
<< waterMass << ";"
<< cumMassLoss << ";"
<< evaporationRate
<< std::endl;
}
return waterMass;
}
template<class SolutionVector, class GridVariables>
void evaluateInterfaceFluxes(const SolutionVector& curSol,
const GridVariables& gridVariables)
{
std::vector<Scalar> faceEvaporation;
for (const auto& element : elements(this->gridGeometry().gridView()))
{
auto fvGeometry = localView(this->gridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, curSol);
for (auto&& scvf : scvfs(fvGeometry))
{
if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
continue;
// TODO: dumux-course-task 2.B
// Use "massCouplingCondition" from the couplingManager here
auto flux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf)
* scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor() * FluidSystem::molarMass(1) * -1.0 * 86400.0;
faceEvaporation.push_back(flux[transportCompIdx]);
}
}
Scalar fluxOverFullSurface = std::accumulate(faceEvaporation.begin(), faceEvaporation.end(), 0.0);
if (exportFluxes_)
{
Dumux::MetaData::Collector collector;
if (Dumux::MetaData::jsonFileExists(fluxFileName_))
Dumux::MetaData::readJsonFile(collector, fluxFileName_);
collector[simulationKey_]["Time"].push_back(time());
collector[simulationKey_]["TimeStepIdx"].push_back(timeStepIndex());
collector[simulationKey_]["FaceEvaporation"].push_back(faceEvaporation);
collector[simulationKey_]["FluxOverFullSurface"].push_back(fluxOverFullSurface);
Dumux::MetaData::writeJsonFile(collector, fluxFileName_);
std::cout << "Total flux per local method = " << fluxOverFullSurface << "\n";
}
}
/*!
* \brief Specifies which kind of boundary condition should be
* used for which equation on a given boundary control volume.
*
* \param element The element
* \param scvf The boundary sub control volume face
*/
BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
{
BoundaryTypes values;
values.setAllNeumann();
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values.setAllCouplingNeumann();
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Neumann control volume.
*
* \param element The element for which the Neumann boundary condition is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables
* \param scvf The boundary sub control volume face
*
* For this method, the \a values variable stores primary variables.
*/
template<class ElementVolumeVariables>
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
{
auto values = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
}
return values;
}
/*!
* \brief Evaluate the source term for all phases within a given
* sub-control-volume.
*
* \param element The element for which the source term is set
* \param fvGeomentry The fvGeometry
* \param elemVolVars The element volume variables
* \param scv The subcontrolvolume
*/
template<class ElementVolumeVariables>
NumEqVector source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolume& scv) const
{ return NumEqVector(0.0); }
/*!
* \brief Evaluate the initial value for a control volume.
*
* \param element The element
*
* For this method, the \a priVars parameter stores primary
* variables.
*/
PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
static const Scalar freeflowPressure = getParamFromGroup<Scalar>("Freeflow", "Problem.Pressure");
PrimaryVariables values(0.0);
// TODO: dumux-course-task 2.A
// Declare here which phases are present.
values.setState(2/*secondPhaseOnly*/);
values[transportCompIdx] = moleFraction_;
values[temperatureIdx] = temperature_;
values[Indices::pressureIdx] = freeflowPressure;
// values[transportCompIdx] = moleFraction_;
return values;
}
//! Set the coupling manager
void setCouplingManager(std::shared_ptr<CouplingManager> cm)
{ couplingManager_ = cm; }
//! Get the coupling manager
const CouplingManager& couplingManager() const
{ return *couplingManager_; }
void setTimeLoop(TimeLoopPtr timeLoop)
{ timeLoop_ = timeLoop; }
Scalar time() const
{ return timeLoop_->time(); }
int timeStepIndex() const
{ return timeLoop_->timeStepIndex(); }
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->gridGeometry().bBoxMin()[0] + eps_; }
bool onRightBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; }
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
Scalar eps_;
Scalar moleFraction_;
Scalar temperature_;
TimeLoopPtr timeLoop_;
std::shared_ptr<CouplingManager> couplingManager_;
// Storage calculation and output
Scalar initialWaterContent_ = 0.0;
Scalar lastWaterMass_ = 0.0;
std::string storageFileName_;
std::ofstream storageFile_;
bool exportStorage_;
bool plotStorage_;
Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
// Flux output
std::string fluxFileName_;
std::string simulationKey_;
bool exportFluxes_;
};
} //end namespace Dumux
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Properties file for coupled models example.
*/
#ifndef DUMUX_EXERCISE_COUPLED_MODELS_PROPERTIES_HH
#define DUMUX_EXERCISE_COUPLED_MODELS_PROPERTIES_HH
// Both Domains
#include <dune/grid/yaspgrid.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include <dumux/io/gnuplotinterface.hh>
#include <dumux/material/fluidsystems/1padapter.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
// Porous medium flow domain
// TODO: dumux-course-task 2.A
// Include 2pnc model here
#include <dumux/porousmediumflow/2pnc/model.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
// TODO: dumux-course-task 2.A
// Include spatial params for a 2-phase system
#include "../2pspatialparams.hh"
#include "porousmediumsubproblem.hh"
// Free-flow
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
#include "freeflowsubproblem.hh"
namespace Dumux::Properties {
// Create new type tags
namespace TTag {
// TODO: dumux-course-task 2.A
// Change to property of the `FluidSystem` such that `H2OAir` is used directly.
struct PorousMediumOnePNC { using InheritsFrom = std::tuple<TwoPNCNI, CCTpfaModel>; };
struct FreeflowNC { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
} // end namespace TTag
// Set the coupling manager
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::FreeflowNC>
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::PorousMediumOnePNC>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::PorousMediumOnePNC>
{
using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeflowNC, Properties::TTag::FreeflowNC, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::PorousMediumOnePNC> { using type = Dumux::PorousMediumSubProblem<TypeTag>; };
template<class TypeTag>
struct Problem<TypeTag, TTag::FreeflowNC> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
// The fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::PorousMediumOnePNC>
{
using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
// using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
using type = H2OAir;
};
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::FreeflowNC>
{
using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
};
// Use moles
template<class TypeTag>
struct UseMoles<TypeTag, TTag::PorousMediumOnePNC> { static constexpr bool value = true; };
template<class TypeTag>
struct UseMoles<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
// Do not replace one equation with a total mass balance
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::PorousMediumOnePNC> { static constexpr int value = 3; };
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::FreeflowNC> { static constexpr int value = 3; };
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::PorousMediumOnePNC> { using type = Dune::YaspGrid<2>; };
template<class TypeTag>
struct Grid<TypeTag, TTag::FreeflowNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
//! Use a model with constant tortuosity for the effective diffusivity
template<class TypeTag>
struct EffectiveDiffusivityModel<TypeTag, TTag::PorousMediumOnePNC>
{ using type = DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>; };
// TODO: dumux-course-task 2.A
// Define new formulation for primary variables here.
template<class TypeTag>
struct Formulation<TypeTag, TTag::PorousMediumOnePNC>
{ static constexpr auto value = TwoPFormulation::p1s0; };
// Set the spatial parameters type
template<class TypeTag>
// TODO: dumux-course-task 2.A
// Adapt the spatial params here.
struct SpatialParams<TypeTag, TTag::PorousMediumOnePNC> {
using type = TwoPSpatialParams<GetPropType<TypeTag, GridGeometry>, GetPropType<TypeTag, Scalar>>;
};
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
} // end namespace Dumux::Properties
#endif
_______________________________________________
DuMux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
