OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Blackoil fluid test programs now compile.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2011-12-22 12:59:42 +01:00
parent 7d705bb763
commit 9b2bfe40d3
10 changed files with 178 additions and 238 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
opm/core/fluid/BlackoilFluid.hpp
View File

@ -21,11 +21,10 @@
#define OPM_BLACKOILFLUID_HEADER_INCLUDED
#include <dune/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp>
#include <dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp>
#include <dune/porsol/blackoil/fluid/BlackoilPVT.hpp>
#include <dune/common/EclipseGridParser.hpp>
#include <dune/common/fvector.hh>
#include <opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp>
#include <opm/core/fluid/blackoil/FluidStateBlackoil.hpp>
#include <opm/core/fluid/blackoil/BlackoilPVT.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <vector>
@ -174,26 +173,34 @@ namespace Opm
states.phase_volume_density.resize(num);
states.total_phase_volume_density.resize(num);
states.saturation.resize(num);
#ifdef COMPUTE_OLD_TERMS
states.phase_compressibility.resize(num);
states.total_compressibility.resize(num);
states.experimental_term.resize(num);
#endif
#pragma omp parallel for
for (int i = 0; i < num; ++i) {
const CompVec& z = states.surface_volume_density[i];
const PhaseVec& B = states.formation_volume_factor[i];
const PhaseVec& dB = states.formation_volume_factor_deriv[i];
const PhaseVec& R = states.solution_factor[i];
const PhaseVec& dR = states.solution_factor_deriv[i];
PhaseToCompMatrix& At = states.state_matrix[i];
PhaseVec& u = states.phase_volume_density[i];
double& tot_phase_vol_dens = states.total_phase_volume_density[i];
PhaseVec& s = states.saturation[i];
#ifdef COMPUTE_OLD_TERMS
const PhaseVec& dB = states.formation_volume_factor_deriv[i];
const PhaseVec& dR = states.solution_factor_deriv[i];
PhaseVec& cp = states.phase_compressibility[i];
double& tot_comp = states.total_compressibility[i];
double& exp_term = states.experimental_term[i];
computeSingleEquilibrium(B, dB, R, dR, z,
At, u, tot_phase_vol_dens,
s, cp, tot_comp, exp_term);
#else
computeSingleEquilibrium(B, R, z,
At, u, tot_phase_vol_dens,
s);
#endif
}
}
@ -276,6 +283,13 @@ namespace Opm
R[Aqua] = 0.0;
R[Vapour] = pvt_.R(p[Vapour], z, Vapour);
R[Liquid] = pvt_.R(p[Liquid], z, Liquid);
// Convenience vars.
PhaseToCompMatrix& At = fluid_state.phase_to_comp_;
PhaseVec& u = fluid_state.phase_volume_density_;
double& tot_phase_vol_dens = fluid_state.total_phase_volume_density_;
PhaseVec& s = fluid_state.saturation_;
#ifdef COMPUTE_OLD_TERMS
PhaseVec dB;
dB[Aqua] = pvt_.dBdp(p[Aqua], z, Aqua);
dB[Vapour] = pvt_.dBdp(p[Vapour], z, Vapour);
@ -284,19 +298,17 @@ namespace Opm
dR[Aqua] = 0.0;
dR[Vapour] = pvt_.dRdp(p[Vapour], z, Vapour);
dR[Liquid] = pvt_.dRdp(p[Liquid], z, Liquid);
// Convenience vars.
PhaseToCompMatrix& At = fluid_state.phase_to_comp_;
PhaseVec& u = fluid_state.phase_volume_density_;
double& tot_phase_vol_dens = fluid_state.total_phase_volume_density_;
PhaseVec& s = fluid_state.saturation_;
PhaseVec& cp = fluid_state.phase_compressibility_;
double& tot_comp = fluid_state.total_compressibility_;
double& exp_term = fluid_state.experimental_term_;
computeSingleEquilibrium(B, dB, R, dR, z,
At, u, tot_phase_vol_dens,
s, cp, tot_comp, exp_term);
#else
computeSingleEquilibrium(B, R, z,
At, u, tot_phase_vol_dens,
s);
#endif
// Compute viscosities.
PhaseVec& mu = fluid_state.viscosity_;
@ -307,6 +319,7 @@ namespace Opm
#ifdef COMPUTE_OLD_TERMS
static void computeSingleEquilibrium(const PhaseVec& B,
const PhaseVec& dB,
const PhaseVec& R,
@ -319,6 +332,15 @@ namespace Opm
PhaseVec& cp,
double& tot_comp,
double& exp_term)
#else
static void computeSingleEquilibrium(const PhaseVec& B,
const PhaseVec& R,
const CompVec& z,
PhaseToCompMatrix& At,
PhaseVec& u,
double& tot_phase_vol_dens,
PhaseVec& s)
#endif
{
// Set the A matrix (A = RB^{-1})
// Using At since we really want Fortran ordering
@ -345,6 +367,7 @@ namespace Opm
s[phase] = u[phase]/tot_phase_vol_dens;
}
#ifdef COMPUTE_OLD_TERMS
// Phase compressibilities.
// PhaseVec& cp = fluid_state.phase_compressibility_[i];
// Set the derivative of the A matrix (A = RB^{-1})
@ -372,6 +395,7 @@ namespace Opm
dAt.mtv(tmp1, tmp2);
Ait.mtv(tmp2, tmp3);
exp_term = tmp3[Aqua] + tmp3[Liquid] + tmp3[Gas];
#endif
}

165
opm/core/fluid/blackoil/BlackoilComponent.hpp
View File

@ -1,165 +0,0 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
#define OPM_BLACKOILCOMPONENT_HEADER_INCLUDED
#include <dune/common/stdstreams.hh>
namespace Dumux
{
/*!
* \ingroup Components
* \brief
* A component class for the black oil model, intended to be used
* for all three components.
*
* \tparam Scalar The type used for scalar values
*/
template <class Scalar>
class BlackoilComponent
{
public:
/*!
* \brief A human readable name for the component.
*/
static const char *name()
{
return "BlackoilComponent";
}
/*!
* \brief The molar mass in [kg] of the component.
*/
static Scalar molarMass()
{ DUNE_THROW(Dune::NotImplemented, "Component::molarMass()"); }
/*!
* \brief Returns the critical temperature in [K] of the component.
*/
static Scalar criticalTemperature()
{ DUNE_THROW(Dune::NotImplemented, "Component::criticalTemperature()"); }
/*!
* \brief Returns the critical pressure in [Pa] of the component.
*/
static Scalar criticalPressure()
{ DUNE_THROW(Dune::NotImplemented, "Component::criticalPressure()"); }
/*!
* \brief Returns the temperature in [K] at the component's triple point.
*/
static Scalar tripleTemperature()
{ DUNE_THROW(Dune::NotImplemented, "Component::tripleTemperature()"); }
/*!
* \brief Returns the pressure in [Pa] at the component's triple point.
*/
static Scalar triplePressure()
{ DUNE_THROW(Dune::NotImplemented, "Component::triplePressure()"); }
/*!
* \brief The vapor pressure in [Pa] of the component at a given
* temperature in [K].
*
* \param T temperature of the component in [K]
*/
static Scalar vaporPressure(Scalar T)
{ DUNE_THROW(Dune::NotImplemented, "Component::vaporPressure()"); }
/*!
* \brief The density in [kg/m^3] of the component at a given pressure in [Pa] and temperature in [K].
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar gasDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::density()"); }
/*!
* \brief The density [kg/m^3] of the liquid component at a given pressure in [Pa] and temperature in [K].
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar liquidDensity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::density()"); }
/*!
* \brief Specific enthalpy [J/kg] of the pure component in gas.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::gasEnthalpy()"); }
/*!
* \brief Specific enthalpy [J/kg] of the pure component in liquid.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::liquidEnthalpy()"); }
/*!
* \brief Specific internal energy [J/kg] of the pure component in gas.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::gasInternalEnergy()"); }
/*!
* \brief Specific internal energy [J/kg] of pure the pure component in liquid.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::liquidInternalEnergy()"); }
/*!
* \brief The dynamic viscosity [Pa*s] of the pure component at a given pressure in [Pa] and temperature in [K].
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar gasViscosity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::gasViscosity()"); }
/*!
* \brief The dynamic liquid viscosity [Pa*s] of the pure component.
*
* \param temperature temperature of component in [K]
* \param pressure pressure of component in [Pa]
*/
static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
{ DUNE_THROW(Dune::NotImplemented, "Component::liquidViscosity()"); }
};
} // end namepace
#endif // OPM_BLACKOILCOMPONENT_HEADER_INCLUDED

64
opm/core/fluid/blackoil/BlackoilDefs.hpp
View File

@ -22,6 +22,7 @@
#include <tr1/array>
#include <iostream>
#include <boost/static_assert.hpp>
namespace Opm
@ -36,16 +37,63 @@ namespace Opm
enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
// We need types with operator= and constructor taking scalars
// for the small vectors and matrices, to save us from having to
// rewrite a large amount of code.
template <typename T, int N>
class SmallVec
{
public:
SmallVec()
{}
SmallVec(const T& elem)
{ data_.assign(elem); } // In C++11, assign -> fill
SmallVec& operator=(const T& elem)
{ data_.assign(elem); return *this; }
const T& operator[](int i) const
{ return data_[i]; }
T& operator[](int i)
{ return data_[i]; }
template <typename U>
void assign(const U& elem)
{
for (int i = 0; i < N; ++i) {
data_[i] = elem;
}
}
private:
std::tr1::array<T, N> data_;
};
template <typename T, int Rows, int Cols>
class SmallMat
{
public:
SmallMat()
{}
SmallMat(const T& elem)
{ data_.assign(elem); } // In C++11, assign -> fill
SmallMat& operator=(const T& elem)
{ data_.assign(elem); return *this; }
typedef SmallVec<T, Cols> RowType;
const RowType& operator[](int i) const
{ return data_[i]; }
RowType& operator[](int i)
{ return data_[i]; }
private:
SmallVec<RowType, Rows> data_;
};
typedef double Scalar;
// typedef Dune::FieldVector<Scalar, numComponents> CompVec;
// typedef Dune::FieldVector<Scalar, numPhases> PhaseVec;
typedef std::tr1::array<Scalar, numComponents> CompVec;
typedef std::tr1::array<Scalar, numPhases> PhaseVec;
typedef SmallVec<Scalar, numComponents> CompVec;
typedef SmallVec<Scalar, numPhases> PhaseVec;
BOOST_STATIC_ASSERT(int(numComponents) == int(numPhases));
// typedef Dune::FieldMatrix<Scalar, numComponents, numPhases> PhaseToCompMatrix;
// typedef Dune::FieldMatrix<Scalar, numPhases, numPhases> PhaseJacobian;
typedef std::tr1::array<PhaseVec, numComponents> PhaseToCompMatrix;
typedef std::tr1::array<PhaseVec, numPhases> PhaseJacobian;
typedef SmallMat<Scalar, numComponents, numPhases> PhaseToCompMatrix;
typedef SmallMat<Scalar, numPhases, numPhases> PhaseJacobian;
// Attempting to guard against alignment issues.
BOOST_STATIC_ASSERT(sizeof(CompVec) == numComponents*sizeof(Scalar));
BOOST_STATIC_ASSERT(sizeof(PhaseVec) == numPhases*sizeof(Scalar));
BOOST_STATIC_ASSERT(sizeof(PhaseToCompMatrix) == numComponents*numPhases*sizeof(Scalar));
BOOST_STATIC_ASSERT(sizeof(PhaseJacobian) == numPhases*numPhases*sizeof(Scalar));
};
} // namespace Opm

8
opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp
View File

@ -20,10 +20,10 @@
#ifndef OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
#define OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
#include <dune/common/EclipseGridParser.hpp>
#include <dune/porsol/common/UniformTableLinear.hpp>
#include <dune/porsol/common/buildUniformMonotoneTable.hpp>
#include "BlackoilDefs.hpp"
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/fluid/blackoil/BlackoilDefs.hpp>
#include <iostream>
#include <stdexcept>

8
opm/core/fluid/blackoil/FluidStateBlackoil.hpp
View File

@ -40,9 +40,11 @@ namespace Opm
PhaseVec solution_factor_;
PhaseToCompMatrix phase_to_comp_; // RB^{-1} in Fortran ordering
PhaseVec saturation_;
#ifdef COMPUTE_OLD_TERMS
PhaseVec phase_compressibility_;
Scalar total_compressibility_;
Scalar experimental_term_;
#endif
PhaseVec viscosity_;
PhaseVec relperm_;
PhaseJacobian drelperm_;
@ -62,9 +64,7 @@ namespace Opm
// Variables from PVT functions.
std::vector<PhaseVec> formation_volume_factor; // B
std::vector<PhaseVec> formation_volume_factor_deriv; // dB/dp
std::vector<PhaseVec> solution_factor; // R
std::vector<PhaseVec> solution_factor_deriv; // dR/dp
std::vector<PhaseVec> viscosity; // mu
// Variables computed from PVT data.
@ -74,9 +74,13 @@ namespace Opm
std::vector<PhaseVec> phase_volume_density; // u
std::vector<Scalar> total_phase_volume_density; // sum(u)
std::vector<PhaseVec> saturation; // s = u/sum(u)
#ifdef COMPUTE_OLD_TERMS
std::vector<PhaseVec> formation_volume_factor_deriv; // dB/dp
std::vector<PhaseVec> solution_factor_deriv; // dR/dp
std::vector<PhaseVec> phase_compressibility; // c
std::vector<Scalar> total_compressibility; // cT
std::vector<Scalar> experimental_term; // ex = sum(Ai*dA*Ai*z)
#endif
// Variables computed from saturation.
std::vector<PhaseVec> relperm; // kr

4
opm/core/fluid/blackoil/MiscibilityDead.cpp
View File

@ -82,7 +82,7 @@ namespace Opm
{
}
double MiscibilityDead::getViscosity(int region, double press, const surfvol_t& /*surfvol*/) const
double MiscibilityDead::getViscosity(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
{
return viscosity_(press);
}
@ -100,7 +100,7 @@ namespace Opm
}
}
double MiscibilityDead::B(int region, double press, const surfvol_t& /*surfvol*/) const
double MiscibilityDead::B(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
{
// Interpolate 1/B
return 1.0/one_over_B_(press);

6
tests/Makefile.am
View File

@ -15,6 +15,9 @@ $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
noinst_PROGRAMS = \
bo_fluid_p_and_z_deps \
bo_fluid_pressuredeps \
bo_fluid_test \
monotcubicinterpolator_test \
param_test \
sparsetable_test \
@ -23,6 +26,9 @@ unit_test \
test_readvector \
test_sf2p
bo_fluid_p_and_z_deps_SOURCES = bo_fluid_p_and_z_deps.cpp
bo_fluid_pressuredeps_SOURCES = bo_fluid_pressuredeps.cpp
bo_fluid_test_SOURCES = bo_fluid_test.cpp
monotcubicinterpolator_test_SOURCES = monotcubicinterpolator_test.cpp
param_test_SOURCES = param_test.cpp
sparsetable_test_SOURCES = sparsetable_test.cpp

12
tests/bo_fluid_p_and_z_deps.cpp
View File

@ -20,9 +20,9 @@
#include "config.h"
#include <dune/common/param/ParameterGroup.hpp>
#include <dune/common/EclipseGridParser.hpp>
#include <dune/porsol/blackoil/BlackoilFluid.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/fluid/BlackoilFluid.hpp>
int main(int argc, char** argv)
@ -49,7 +49,11 @@ int main(int argc, char** argv)
int changing_component = param.getDefault("changing_component", int(Opm::BlackoilFluid::Gas));
double min_z = param.getDefault("min_z", 0.0);
double max_z = param.getDefault("max_z", 500.0);
#ifdef COMPUTE_OLD_TERMS
int variable = param.getDefault("variable", 0);
#else
int variable = param.getDefault("variable", 2);
#endif
Opm::BlackoilFluid::CompVec z = z0;
std::cout << "%}\n"
<< "data = [\n";
@ -67,12 +71,14 @@ int main(int argc, char** argv)
std::cout.fill(' ');
double var = 0.0;
switch (variable) {
#ifdef COMPUTE_OLD_TERMS
case 0:
var = state.total_compressibility_;
break;
case 1:
var = state.experimental_term_;
break;
#endif
case 2:
var = state.saturation_[0];
break;

87
tests/bo_fluid_pressuredeps.cpp
View File

@ -19,10 +19,51 @@
#include "config.h"
#include <dune/common/param/ParameterGroup.hpp>
#include <dune/common/EclipseGridParser.hpp>
#include <dune/porsol/blackoil/BlackoilFluid.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/fluid/BlackoilFluid.hpp>
#ifdef COMPUTE_OLD_TERMS
void printCompressibilityTerms(double p, const Opm::BlackoilFluid::FluidState& state)
{
std::cout.precision(6);
std::cout.width(15);
std::cout.fill(' ');
std::cout << p << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.total_compressibility_ << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << totmob << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.total_phase_volume_density_ << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.experimental_term_ << '\n';
}
#endif
void printSatsEtc(double p, const Opm::BlackoilFluid::FluidState& state)
{
std::cout.precision(6);
std::cout.width(15);
std::cout.fill(' ');
std::cout << p << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.saturation_[0] << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.saturation_[1] << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.saturation_[2] << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.total_phase_volume_density_ << '\n';
}
int main(int argc, char** argv)
{
@ -43,7 +84,9 @@ int main(int argc, char** argv)
int num_pts = param.getDefault("num_pts", 41);
double min_press = param.getDefault("min_press", 1e7);
double max_press = param.getDefault("max_press", 3e7);
#ifdef COMPUTE_OLD_TERMS
bool print_compr = param.getDefault("print_compr", true);
#endif
for (int i = 0; i < num_pts; ++i) {
double factor = double(i)/double(num_pts - 1);
double p = (1.0 - factor)*min_press + factor*max_press;
@ -52,40 +95,14 @@ int main(int argc, char** argv)
for (int phase = 0; phase < Opm::BlackoilFluid::numPhases; ++phase) {
totmob += state.mobility_[phase];
}
#ifdef COMPUTE_OLD_TERMS
if (print_compr) {
std::cout.precision(6);
std::cout.width(15);
std::cout.fill(' ');
std::cout << p << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.total_compressibility_ << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << totmob << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.total_phase_volume_density_ << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.experimental_term_ << '\n';
printCompressibilityTerms(p, state);
} else {
std::cout.precision(6);
std::cout.width(15);
std::cout.fill(' ');
std::cout << p << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.saturation_[0] << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.saturation_[1] << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.saturation_[2] << " ";
std::cout.width(15);
std::cout.fill(' ');
std::cout << state.total_phase_volume_density_ << '\n';
printSatsEtc(p, state);
}
#else
printSatsEtc(p, state);
#endif
}
}

10
tests/bo_fluid_test.cpp
View File

@ -19,9 +19,9 @@
#include "config.h"
#include <dune/porsol/blackoil/fluid/FluidMatrixInteractionBlackoil.hpp>
#include <dune/common/param/ParameterGroup.hpp>
#include <dune/common/fvector.hh>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp>
int main(int argc, char** argv)
@ -37,7 +37,7 @@ int main(int argc, char** argv)
Opm::FluidMatrixInteractionBlackoilParams<double> fluid_params;
fluid_params.init(parser);
typedef Opm::FluidMatrixInteractionBlackoil<double> Law;
Dune::FieldVector<double, 3> s, kr;
Law::PhaseVec s, kr;
const double temp = 300; // [K]
int num = 41;
for (int i = 0; i < num; ++i) {
@ -47,7 +47,7 @@ int main(int argc, char** argv)
Law::kr(kr, fluid_params, s, temp);
std::cout.width(6);
std::cout.fill(' ');
std::cout << s[Law::Liquid] << " " << kr << '\n';
std::cout << s[Law::Liquid] << " " << kr[0] << ' ' << kr[1] << ' ' << kr[2] << '\n';
}
}