Commit d8e101d3 authored by Katharina Heck's avatar Katharina Heck
Browse files

add the four examples

parents
Pipeline #1683 failed with stages
# build system clutter
build-*
Testing
# auto-saved files
*~
# hidden files
.cproject
.project
# left overs from git rebase
*.orig
*.rej
# latex clutter
*.pdf
*.aux
*.blg
*.log
*.bbl
*.dvi
*.idx
*.out
*.tdo
*.toc
*.synctex.gz
# always consider files containing source code regardless of their name
!*.cc
!*.hh
!*.c
!*.h
!*.sh
!*.py
# always consider reference solutions
!*reference.vtu
# We require version CMake version 3.1 to prevent issues
# with dune_enable_all_packages and older CMake versions.
cmake_minimum_required(VERSION 3.1)
project(Heck2020a CXX)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
${PROJECT_BINARY_DIR})
endif()
#find dune-common and set the module path
find_package(dune-common REQUIRED)
list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
${dune-common_MODULE_PATH})
#include the dune macros
include(DuneMacros)
# start a dune project with information from dune.module
dune_project()
dune_enable_all_packages()
add_subdirectory(doc)
add_subdirectory(cmake/modules)
add_subdirectory(appl)
add_subdirectory(dumux)
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
prefix=@prefix@
exec_prefix=@exec_prefix@
libdir=@libdir@
includedir=@includedir@
CXX=@CXX@
CC=@CC@
DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: Heck2020a module
URL: http://dune-project.org/
Requires: dumux
Libs: -L${libdir}
Cflags: -I${includedir}
Preparing the Sources
=========================
Additional to the software mentioned in README you'll need the
following programs installed on your system:
cmake >= 2.8.12
Getting started
---------------
If these preliminaries are met, you should run
dunecontrol all
which will find all installed dune modules as well as all dune modules
(not installed) which sources reside in a subdirectory of the current
directory. Note that if dune is not installed properly you will either
have to add the directory where the dunecontrol script resides (probably
./dune-common/bin) to your path or specify the relative path of the script.
Most probably you'll have to provide additional information to dunecontrol
(e. g. compilers, configure options) and/or make options.
The most convenient way is to use options files in this case. The files
define four variables:
CMAKE_FLAGS flags passed to cmake (during configure)
An example options file might look like this:
#use this options to configure and make if no other options are given
CMAKE_FLAGS=" \
-DCMAKE_CXX_COMPILER=g++-5 \
-DCMAKE_CXX_FLAGS='-Wall -pedantic' \
-DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-5 and set compiler flags
If you save this information into example.opts you can pass the opts file to
dunecontrol via the --opts option, e. g.
dunecontrol --opts=example.opts all
More info
---------
See
dunecontrol --help
for further options.
The full build system is described in the dune-common/doc/buildsystem (Git version) or under share/doc/dune-common/buildsystem if you installed DUNE!
add_subdirectory(example1)
add_subdirectory(example2)
add_subdirectory(example3)
add_subdirectory(example4)
add_input_file_links()
dune_add_test(NAME test_stokes1p2cdarcy2p2c_radiation
SOURCES main.cc
CMAKE_GUARD HAVE_UMFPACK)
set(CMAKE_BUILD_TYPE Debug)
// -*- 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 2 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_DARCY_2P2C_RADIATION_SUBPROBLEM_HH
#define DUMUX_DARCY_2P2C_RADIATION_SUBPROBLEM_HH
#include <dune/grid/yaspgrid.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/porousmediumflow/2p2c/model.hh>
#include <dumux/porousmediumflow/problem.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/io/gnuplotinterface.hh>
#include "spatialparams.hh"
#include "radiation.hh"
namespace Dumux
{
template <class TypeTag>
class DarcySubProblem;
namespace Properties
{
// Create new type tags
namespace TTag {
struct DarcyTwoPTwoCTypeTag { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
} // end namespace TTag
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
// the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
//! Set the default formulation to pw-Sn: This can be over written in the problem.
template<class TypeTag>
struct Formulation<TypeTag, TTag::DarcyTwoPTwoCTypeTag>
{ static constexpr auto value = TwoPFormulation::p1s0; };
//// The gas component balance (air) is replaced by the total mass balance
template<class TypeTag>
struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { static constexpr int value = 5; };
// Set the grid type
template<class TypeTag>
struct Grid<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
template<class TypeTag>
struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { static constexpr bool value = true; };
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = TwoPTwoCSpatialParams<TypeTag>; };
template<class TypeTag>
struct ThermalConductivityModel<TypeTag, TTag::DarcyTwoPTwoCTypeTag> { using type = ThermalConductivitySomerton<GetPropType<TypeTag, Properties::Scalar>>; };
}
template <class TypeTag>
class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = GetPropType<TypeTag, Properties::GridView>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
// copy some indices for convenience
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
enum {
// primary variable indices
conti0EqIdx = Indices::conti0EqIdx,
contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx,
pressureIdx = Indices::pressureIdx,
switchIdx = Indices::switchIdx
};
enum {energyEqIdx = Indices::energyEqIdx}; //water temp
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
using Radiation = typename Dumux::Radiation<Scalar, FVGridGeometry>;
public:
DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<CouplingManager> couplingManager)
: ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
{
pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
temperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Temperature");
initialPhasePresence_ = getParamFromGroup<int>(this->paramGroup(), "Problem.InitPhasePresence");
time_ = 0.0;
initializationTime_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.InitializationTime");
}
void setTime(Scalar time)
{ time_ = time; }
const Scalar time() const
{return time_; }
/*!
* \name Simulation steering
*/
// \{
template<class SolutionVector, class GridVariables>
void postTimeStep(const SolutionVector& curSol,
const GridVariables& gridVariables,
const Scalar time)
{
if (time_>=initializationTime_)
{
// compute the mass in the entire domain
Scalar massWater = 0.0;
Scalar numScvf = 0.0;
Scalar averageTemperature = 0.0;
Scalar evaporation = 0.0;
Scalar energyFlux = 0.0;
Scalar netRadiation = 0.0;
// bulk elements
for (const auto& element : elements(this->fvGridGeometry().gridView()))
{
auto fvGeometry = localView(this->fvGridGeometry());
fvGeometry.bindElement(element);
auto elemVolVars = localView(gridVariables.curGridVolVars());
elemVolVars.bindElement(element, fvGeometry, curSol);
for (auto&& scv : scvs(fvGeometry))
{
const auto& volVars = elemVolVars[scv];
for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
{
massWater += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx)*volVars.density(phaseIdx)
* scv.volume() * volVars.saturation(phaseIdx) * volVars.porosity() * volVars.extrusionFactor();
}
}
for (auto&& scvf : scvfs(fvGeometry))
{
if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
continue;
// NOTE: binding the coupling context is necessary
couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
const auto flux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
const auto fluxEnergy = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf);
const auto radiation = Radiation::radationEquilibrium(element, elemVolVars[scvf.insideScvIdx()], scvf, couplingManager().couplingData().stokesTemperature(element, scvf), time);
evaporation += flux[0] * scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor();
energyFlux += fluxEnergy* scvf.area() * elemVolVars[scvf.insideScvIdx()].extrusionFactor();
averageTemperature += elemVolVars[scvf.insideScvIdx()].temperature();
numScvf += 1;
netRadiation = radiation;
}
}
// convert to kg/s if using mole fractions
evaporation = evaporation * FluidSystem::molarMass(FluidSystem::H2OIdx);
averageTemperature = averageTemperature/numScvf;
std::cout<<std::setprecision(12)<<"mass of water is: " << massWater << std::endl;
std::cout<<"evaporation from pm "<<evaporation<<std::endl;
std::cout<<"energyFlux from pm "<<energyFlux<<std::endl;
std::cout<<"radiation "<<netRadiation<<std::endl;
//do a gnuplot
x_.push_back(time); // in seconds
y_.push_back(evaporation);
gnuplot_.resetPlot();
gnuplot_.setXRange(0,std::max(time, 86400.0));
gnuplot_.setYRange(0, 1e-5);
gnuplot_.setXlabel("time [s]");
gnuplot_.setYlabel("kg/s");
gnuplot_.addDataSetToPlot(x_, y_, "evaporation");
gnuplot_.plot("evaporation");
//do a gnuplot
y2_.push_back(massWater);
gnuplot2_.resetPlot();
gnuplot2_.setXRange(0,std::max(time, 86400.0));
gnuplot2_.setYRange(0, 20);
gnuplot2_.setXlabel("time [s]");
gnuplot2_.setYlabel("kg");
gnuplot2_.addDataSetToPlot(x_, y2_, "water mass");
gnuplot2_.plot("watermass");
//do a gnuplot
y3_.push_back(averageTemperature);
gnuplot3_.resetPlot();
gnuplot3_.setXRange(0,std::max(time, 86400.0));
gnuplot3_.setYRange(280,300);
gnuplot3_.setXlabel("time [s]");
gnuplot3_.setYlabel("K");
gnuplot3_.addDataSetToPlot(x_, y3_, "temperatureAverage");
gnuplot3_.plot("temperature");
//do a gnuplot
y4_.push_back(netRadiation);
gnuplot4_.resetPlot();
gnuplot4_.setXRange(0,std::max(time, 86400.0));
gnuplot4_.setYRange(0,600);
gnuplot4_.setXlabel("time [s]");
gnuplot4_.setYlabel("W/m^2");
gnuplot4_.addDataSetToPlot(x_, y4_, "net radiation");
gnuplot4_.plot("radiation");
}
}
/*!
* \brief Returns true if a restart file should be written to
* disk.
*/
bool shouldWriteRestartFile() const
{ return false; }
/*!
* \name Problem parameters
*/
// \{
bool shouldWriteOutput() const //define output
{ return true; }
// \}
/*!
* \name Boundary conditions
*/
// \{
/*!
* \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;
const auto& globalPos = scvf.center();
values.setAllNeumann();
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values.setAllCouplingNeumann();
// else if (onLowerBoundary_(globalPos))
// values.setAllDirichlet();
else
values.setAllNeumann();
return values;
}
/*!
* \brief Evaluate the boundary conditions for a Dirichlet control volume.
*
* \param element The element for which the Dirichlet boundary condition is set
* \param scvf The boundary subcontrolvolumeface
*
* For this method, the \a values parameter stores primary variables.
*/
PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
{
PrimaryVariables values(0.0);
values = initialAtPos(scvf.center());
// values[switchIdx]= 0.8;
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
*
*/
NumEqVector neumann(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf) const
{
NumEqVector values(0.0);
const auto& volVars = elemVolVars[scvf.insideScvIdx()];
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf) &&time_>=initializationTime_)
{
const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
for(int i = 0; i< massFlux.size(); ++i)
values[i] = massFlux[i];
values[energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf);
values[energyEqIdx] += -1*Radiation::radationEquilibrium(element, elemVolVars[scvf.insideScvIdx()], scvf, couplingManager().couplingData().stokesTemperature(element, scvf),time());
}
else
{
Scalar temperatureInside = volVars.temperature();
Scalar plexiglassThermalConductivity =
getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PlexiglassThermalConductivity");
Scalar plexiglassThickness = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PlexiglassThickness");
values[energyEqIdx] = plexiglassThermalConductivity
* (temperatureInside - temperature_)
/ plexiglassThickness;
}
return values;
}
// \}
/*!
* \name Volume terms
*/
// \{
/*!
* \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
*
* For this method, the \a values variable stores the rate mass
* of a component is generated or annihilate per volume
* unit. Positive values mean that mass is created, negative ones
* mean that it vanishes.
*/
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.
*
* For this method, the \a priVars parameter stores primary
* variables.
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
values.setState(initialPhasePresence_);
values[pressureIdx] = pressure_ + 1. * this->spatialParams().gravity(globalPos)[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
values[switchIdx] = initialSw_;
values[energyEqIdx] = temperature_; //20
return values;
}
// \}
/*!
* \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_; }
private:
bool onLeftBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
bool onRightBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
bool onLowerBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
bool onUpperBoundary_(const GlobalPosition &globalPos) const
{ return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
Scalar pressure_;
Scalar initialSw_;
Scalar temperature_;
int initialPhasePresence_;
Scalar eps_;
Scalar time_;
Scalar initializationTime_;
std::shared_ptr<CouplingManager> couplingManager_;
std::vector<double> x_;
std::vector<double> y_;
std::vector<double> y2_;
std::vector<double> y3_;
std::vector<double> y4_;
Dumux::GnuplotInterface<double> gnuplot_;
Dumux::GnuplotInterface<double> gnuplot2_;
Dumux::GnuplotInterface<double> gnuplot3_;
Dumux::GnuplotInterface<double> gnuplot4_;
};
} //end namespace
#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH