plotproperties.cc 21.44 KiB
// -*- 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 MaterialTests
* \brief Plot properties of components and fluids.
*/
#include "config.h"
#include <array>
#include <algorithm>
#include <cstring>
#include <limits>
#include <vector>
#include <dumux/common/typetraits/isvalid.hh>
#include <dumux/common/typetraits/typetraits.hh>
#include <dumux/common/parameters.hh>
#include <dumux/io/gnuplotinterface.hh>
#include <dumux/material/components/air.hh>
#include <dumux/material/components/ammonia.hh>
#include <dumux/material/components/benzene.hh>
#include <dumux/material/components/brine.hh>
#include <dumux/material/components/calcite.hh>
#include <dumux/material/components/calciumion.hh>
#include <dumux/material/components/cao.hh>
#include <dumux/material/components/cao2h2.hh>
#include <dumux/material/components/carbonateion.hh>
#include <dumux/material/components/ch4.hh>
#include <dumux/material/components/chlorideion.hh>
#include <dumux/material/components/constant.hh>
#include <dumux/material/components/co2.hh>
#include <dumux/material/components/glucose.hh>
#include <dumux/material/components/granite.hh>
#include <dumux/material/components/h2.hh>
#include <dumux/material/components/h2o.hh>
#include <dumux/material/components/heavyoil.hh>
#include <dumux/material/components/mesitylene.hh>
#include <dumux/material/components/n2.hh>
#include <dumux/material/components/nacl.hh>
#include <dumux/material/components/o2.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/sodiumion.hh>
#include <dumux/material/components/trichloroethene.hh>
#include <dumux/material/components/urea.hh>
#include <dumux/material/components/xylene.hh>
#include <dumux/material/components/componenttraits.hh>
using namespace std;
using namespace Dumux;
namespace Dumux {
//! Helper struct to deactivate static assertions in component's base classes.
struct DisableStaticAssert {};
/*!
* \brief Specialization of Dumux::AlwaysFalse for the struct defined
* above. This is done in order to deactivate the static_assert in
* the base classes of components. If the base class function is compiled
* we do not call it (see below).
*/
template<>
struct AlwaysFalse<DisableStaticAssert> : public std::true_type {};
} // end namespace Dumux
//! Helper structs for detecting if a component has certain functions overloaded
struct checkLiqDen { template<class C> auto operator()(C&& c) -> decltype(C::template liquidDensity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkLiqEnth { template<class C> auto operator()(C&& c) -> decltype(C::template liquidEnthalpy<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkLiqHeatCap { template<class C> auto operator()(C&& c) -> decltype(C::template liquidHeatCapacity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkLiqVisc { template<class C> auto operator()(C&& c) -> decltype(C::template liquidViscosity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkLiqThermCond { template<class C> auto operator()(C&& c) -> decltype(C::template liquidThermalConductivity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkGasDen { template<class C> auto operator()(C&& c) -> decltype(C::template gasDensity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkGasEnth { template<class C> auto operator()(C&& c) -> decltype(C::template gasEnthalpy<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkGasHeatCap { template<class C> auto operator()(C&& c) -> decltype(C::template gasHeatCapacity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkGasVisc { template<class C> auto operator()(C&& c) -> decltype(C::template gasViscosity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkGasThermCond { template<class C> auto operator()(C&& c) -> decltype(C::template gasThermalConductivity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkSolDen { template<class C> auto operator()(C&& c) -> decltype(C::template solidDensity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkSolHeatCap { template<class C> auto operator()(C&& c) -> decltype(C::template solidHeatCapacity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkSolThermCond { template<class C> auto operator()(C&& c) -> decltype(C::template solidThermalConductivity<DisableStaticAssert>(0.0, 0.0)) {} };
struct checkIonCharge { template<class C> auto operator()(C&& c) -> decltype(C::template charge<DisableStaticAssert>(0.0, 0.0)) {} };
//! Plot given values
template<class Functor>
void plot(Functor&& f,
const vector<double>& T,
const double pressure,
const std::string& compName,
const std::string& phaseName,
const std::string& propName,
const std::string& unit,
bool openPlot)
{
vector<double> values(T.size());
for (int i = 0; i < T.size(); ++i)
values[i] = f(T[i], pressure);
const auto minMax = minmax_element(values.begin(), values.end());
Dumux::GnuplotInterface<double> gnuplot(true);
gnuplot.setOpenPlotWindow(openPlot);
gnuplot.setCreateImage(true);
gnuplot.setXRange(T[0], T[T.size()-1]);
gnuplot.setYRange(*(minMax.first)*0.999, *(minMax.second)*1.001);
gnuplot.setXlabel("temperature [K]");
gnuplot.setYlabel(phaseName + " " + propName + " " + unit);
gnuplot.setDatafileSeparator(',');
gnuplot.addDataSetToPlot(T, values, compName + "_" + phaseName + "_" + propName + ".csv");
gnuplot.plot(compName + "_" + phaseName + "_" + propName);
}
//! Plot properties if overloads compile
template<class C, class hasNoDensityOverload = checkLiqDen>
auto plotLiquidDensity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoDensityOverload{})(declval<C>()))::value && ComponentTraits<C>::hasLiquidState, void>
{
auto f = [] (auto T, auto p) { return C::liquidDensity(T, p); };
plot(f, T, p, C::name(), "liquid", "density", "[kg/^3]", openPlot);
}
template<class C, class hasNoEnthalpyOverload = checkLiqEnth>
auto plotLiquidEnthalpy(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoEnthalpyOverload{})(declval<C>()))::value && ComponentTraits<C>::hasLiquidState, void>
{
auto f = [] (auto T, auto p) { return C::liquidEnthalpy(T, p); };
plot(f, T, p, C::name(), "liquid", "enthalpy", "[J/(kg)]", openPlot);}
template<class C, class hasNoHeatCapOverload = checkLiqHeatCap>
auto plotLiquidHeatCapacity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoHeatCapOverload{})(declval<C>()))::value && ComponentTraits<C>::hasLiquidState, void>
{
auto f = [] (auto T, auto p) { return C::liquidHeatCapacity(T, p); };
plot(f, T, p, C::name(), "liquid", "heat capacity", "[J/(kg*K)]", openPlot);
}
template<class C, class hasNoViscOverload = checkLiqVisc>
auto plotLiquidViscosity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoViscOverload{})(declval<C>()))::value && ComponentTraits<C>::hasLiquidState, void>
{
auto f = [] (auto T, auto p) { return C::liquidViscosity(T, p); };
plot(f, T, p, C::name(), "liquid", "viscosity", "[Pa*s]", openPlot);
}
template<class C, class hasNoThermCondOverload = checkLiqThermCond>
auto plotLiquidThermalConductivity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoThermCondOverload{})(declval<C>()))::value && ComponentTraits<C>::hasLiquidState, void>
{
auto f = [] (auto T, auto p) { return C::liquidThermalConductivity(T, p); };
plot(f, T, p, C::name(), "liquid", "thermal conductivity", "[J/(kg*K)]", openPlot);
}
template<class C, class hasNoDensityOverload = checkGasDen>
auto plotGasDensity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoDensityOverload{})(declval<C>()))::value && ComponentTraits<C>::hasGasState, void>
{
auto f = [] (auto T, auto p) { return C::gasDensity(T, p); };
plot(f, T, p, C::name(), "gas", "density", "[kg/^3]", openPlot);
}
template<class C, class hasNoEnthalpyOverload = checkGasEnth>
auto plotGasEnthalpy(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoEnthalpyOverload{})(declval<C>()))::value && ComponentTraits<C>::hasGasState, void>
{
auto f = [] (auto T, auto p) { return C::gasEnthalpy(T, p); };
plot(f, T, p, C::name(), "gas", "enthalpy", "[J/(kg)]", openPlot);
}
template<class C, class hasNoHeatCapOverload = checkGasHeatCap>
auto plotGasHeatCapacity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoHeatCapOverload{})(declval<C>()))::value && ComponentTraits<C>::hasGasState, void>
{
auto f = [] (auto T, auto p) { return C::gasHeatCapacity(T, p); };
plot(f, T, p, C::name(), "gas", "heat capacity", "[J/(kg*K)]", openPlot);
}
template<class C, class hasNoViscOverload = checkGasVisc>
auto plotGasViscosity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoViscOverload{})(declval<C>()))::value && ComponentTraits<C>::hasGasState, void>
{
auto f = [] (auto T, auto p) { return C::gasViscosity(T, p); };
plot(f, T, p, C::name(), "gas", "viscosity", "[Pa*s]", openPlot);
}
template<class C, class hasNoThermCondOverload = checkGasThermCond>
auto plotGasThermalConductivity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoThermCondOverload{})(declval<C>()))::value && ComponentTraits<C>::hasGasState, void>
{
auto f = [] (auto T, auto p) { return C::gasThermalConductivity(T, p); };
plot(f, T, p, C::name(), "gas", "thermal conductivity", "[J/(kg*K)]", openPlot);
}
template<class C, class hasNoDensityOverload = checkSolDen>
auto plotSolidDensity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoDensityOverload{})(declval<C>()))::value && ComponentTraits<C>::hasSolidState, void>
{ auto f = [] (auto T, auto p) { return C::solidDensity(T); };
plot(f, T, p, C::name(), "solid", "density", "[kg/^3]", openPlot);
}
template<class C, class hasNoHeatCapOverload = checkSolHeatCap>
auto plotSolidHeatCapacity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoHeatCapOverload{})(declval<C>()))::value && ComponentTraits<C>::hasSolidState, void>
{
auto f = [] (auto T, auto p) { return C::solidHeatCapacity(T); };
plot(f, T, p, C::name(), "solid", "heat capacity", "[J/(kg*K)]", openPlot);
}
template<class C, class hasNoThermCondOverload = checkSolThermCond>
auto plotSolidThermalConductivity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoThermCondOverload{})(declval<C>()))::value && ComponentTraits<C>::hasSolidState, void>
{
auto f = [] (auto T, auto p) { return C::solidThermalConductivity(T); };
plot(f, T, p, C::name(), "solid", "thermal conductivity", "[J/(kg*K)]", openPlot);
}
template<class C, class hasNoChargeOverload = checkIonCharge>
auto plotIonCharge(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<!decltype(isValid(hasNoChargeOverload{})(declval<C>()))::value && ComponentTraits<C>::isIon, void>
{
auto f = [] (auto T, auto p) { return C::charge(); };
plot(f, T, p, C::name(), "ion", "charge", "[e]", openPlot);
}
//! Do not plot properties if overloads don't compile
template<class C, class hasNoDensityOverload = checkLiqDen>
auto plotLiquidDensity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoDensityOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasLiquidState, void> {}
template<class C, class hasNoEnthalpyOverload = checkLiqEnth>
auto plotLiquidEnthalpy(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoEnthalpyOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasLiquidState, void> {}
template<class C, class hasNoHeatCapOverload = checkLiqHeatCap>
auto plotLiquidHeatCapacity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoHeatCapOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasLiquidState, void> {}
template<class C, class hasNoViscOverload = checkLiqVisc>
auto plotLiquidViscosity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoViscOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasLiquidState, void> {}
template<class C, class hasNoThermCondOverload = checkLiqThermCond>
auto plotLiquidThermalConductivity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoThermCondOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasLiquidState, void> {}
template<class C, class hasNoDensityOverload = checkGasDen>
auto plotGasDensity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoDensityOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasGasState, void> {}
template<class C, class hasNoEnthalpyOverload = checkGasEnth>
auto plotGasEnthalpy(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoEnthalpyOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasGasState, void> {}
template<class C, class hasNoHeatCapOverload = checkGasHeatCap>
auto plotGasHeatCapacity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoHeatCapOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasGasState, void> {}
template<class C, class hasNoViscOverload = checkGasVisc>
auto plotGasViscosity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoViscOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasGasState, void> {}
template<class C, class hasNoThermCondOverload = checkGasThermCond>
auto plotGasThermalConductivity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoThermCondOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasGasState, void> {}
template<class C, class hasNoDensityOverload = checkSolDen>auto plotSolidDensity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoDensityOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasSolidState, void> {}
template<class C, class hasNoHeatCapOverload = checkSolHeatCap>
auto plotSolidHeatCapacity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoHeatCapOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasSolidState, void> {}
template<class C, class hasNoThermCondOverload = checkSolThermCond>
auto plotSolidThermalConductivity(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoThermCondOverload{})(declval<C>()))::value || !ComponentTraits<C>::hasSolidState, void> {}
template<class C, class hasNoChargeOverload = checkIonCharge>
auto plotIonCharge(const vector<double>& T, double p, bool openPlot)
-> typename std::enable_if_t<decltype(isValid(hasNoChargeOverload{})(declval<C>()))::value || !ComponentTraits<C>::isIon, void> {}
//! A number of properties of a component
template<class Component>
void plotStuff(bool openPlotWindow)
{
double pressure = 1e5;
double TMin = 273.15;
double TMax = 323.15;
double TRange = TMax - TMin;
const unsigned int numIntervals = 100;
vector<double> T(numIntervals + 1);
for (int i = 0; i <= numIntervals; i++)
T[i] = TMin + TRange * double(i) /double(numIntervals);
plotLiquidDensity<Component>(T, pressure, openPlotWindow);
plotLiquidEnthalpy<Component>(T, pressure, openPlotWindow);
plotLiquidHeatCapacity<Component>(T, pressure, openPlotWindow);
plotLiquidViscosity<Component>(T, pressure, openPlotWindow);
plotLiquidThermalConductivity<Component>(T, pressure, openPlotWindow);
plotGasDensity<Component>(T, pressure, openPlotWindow);
plotGasEnthalpy<Component>(T, pressure, openPlotWindow);
plotGasHeatCapacity<Component>(T, pressure, openPlotWindow);
plotGasViscosity<Component>(T, pressure, openPlotWindow);
plotGasThermalConductivity<Component>(T, pressure, openPlotWindow);
plotSolidDensity<Component>(T, pressure, openPlotWindow);
plotSolidThermalConductivity<Component>(T, pressure, openPlotWindow);
plotSolidHeatCapacity<Component>(T, pressure, openPlotWindow);
plotIonCharge<Component>(T, pressure, openPlotWindow);
}
////////////////////////
// the main function
////////////////////////
int main(int argc, char *argv[])
{
using namespace Dumux;
bool openPlotWindow = false;
if (argc == 3 && (strcmp(argv[2], "1") || strcmp(argv[2], "true") || strcmp(argv[2], "True")))
openPlotWindow = true;
if (argc < 2)
DUNE_THROW(Dune::InvalidStateException, "At least one argument (the component name) is required!");
const std::string compName = argv[1];
if (compName == "Air")
plotStuff< Components::Air<double> >(openPlotWindow);
else if (compName == "Ammonia")
plotStuff< Components::Ammonia<double> >(openPlotWindow);
else if (compName == "Benzene")
plotStuff< Components::Benzene<double> >(openPlotWindow);
else if (compName == "Brine") {
Parameters::init([](auto& params){ params["Brine.Salinity"] = "0.1"; });
plotStuff< Components::Brine<double> >(openPlotWindow);
}
else if (compName == "Calcite")
plotStuff< Components::Calcite<double> >(openPlotWindow);
else if (compName == "CalciumIon")
plotStuff< Components::CalciumIon<double> >(openPlotWindow);
else if (compName == "CaO")
plotStuff< Components::CaO<double> >(openPlotWindow);
else if (compName == "CaO2H2")
plotStuff< Components::CaO2H2<double> >(openPlotWindow);
else if (compName == "CarbonateIon")
plotStuff< Components::CarbonateIon<double> >(openPlotWindow);
else if (compName == "CH4")
plotStuff< Components::CH4<double> >(openPlotWindow);
else if (compName == "ChlorideIon")
plotStuff< Components::ChlorideIon<double> >(openPlotWindow);
else if (compName == "Constant")
{
Parameters::init([](auto& params){
params["Component.LiquidDensity"] = "1e3";
params["Component.LiquidKinematicViscosity"] = "1e-3";
params["Component.LiquidThermalConductivity"] = "0.679";
params["Component.LiquidHeatCapacity"] = "4.2e3";
params["Component.GasDensity"] = "1";
params["Component.GasKinematicViscosity"] = "1";
params["Component.SolidDensity"] = "1e3";
params["Component.SolidThermalConductivity"] = "0.679";
params["Component.SolidHeatCapacity"] = "4.2e3";
});
plotStuff< Components::Constant<1, double> >(openPlotWindow);
}
else if (compName == "Glucose")
plotStuff< Components::Glucose<double> >(openPlotWindow);
else if (compName == "Granite")
plotStuff< Components::Granite<double> >(openPlotWindow);
else if (compName == "H2")
plotStuff< Components::H2<double> >(openPlotWindow);
else if (compName == "H2O")
plotStuff< Components::H2O<double> >(openPlotWindow);
else if (compName == "HeavyOil")
plotStuff< Components::HeavyOil<double> >(openPlotWindow);
else if (compName == "Mesitylene")
plotStuff< Components::Mesitylene<double> >(openPlotWindow);
else if (compName == "N2")
plotStuff< Components::N2<double> >(openPlotWindow);
else if (compName == "NaCl")
plotStuff< Components::NaCl<double> >(openPlotWindow);
else if (compName == "O2")
plotStuff< Components::O2<double> >(openPlotWindow);
else if (compName == "SimpleH2O")
plotStuff< Components::SimpleH2O<double> >(openPlotWindow);
else if (compName == "SodiumIon")
plotStuff< Components::SodiumIon<double> >(openPlotWindow);
else if (compName == "Trichloroethene")
plotStuff< Components::Trichloroethene<double> >(openPlotWindow);
else if (compName == "Urea")
plotStuff< Components::Urea<double> >(openPlotWindow);
else if (compName == "Xylene")
plotStuff< Components::Xylene<double> >(openPlotWindow);
else
DUNE_THROW(Dune::NotImplemented, "Test for component " << compName);
}