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

simplify the updateDerivativeSinglephase

Browse Source 
and also finishing the singlephase test.
This commit is contained in:
Kai Bao 2022-07-05 14:20:59 +02:00
parent 6da7cab570
commit eccce0209c
2 changed files with 73 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
opm/material/constraintsolvers/PTFlash.hpp
View File

@ -171,7 +171,6 @@ public:
}
// the flash solution process were performed in scalar form, after the flash calculation finishes,
// we transform the derivatives to the final solution
// ensure that things in fluid_state_scalar is transformed to fluid_state
for (int compIdx=0; compIdx<numComponents; ++compIdx){
const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, compIdx);
@ -185,6 +184,7 @@ public:
fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
}
fluid_state.setLvalue(L_scalar);
// we update the derivatives in fluid_state
updateDerivatives_(fluid_state_scalar, z, fluid_state, is_single_phase);
}//end solve
@ -1109,52 +1109,18 @@ protected:
template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
static void updateDerivativesSinglePhase_(const FlashFluidStateScalar& fluid_state_scalar,
const ComponentVector& z,
FluidState& fluid_state)
const ComponentVector& z,
FluidState& fluid_state)
{
constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
constexpr size_t num_deri = numComponents;
using InputEval = typename FluidState::Scalar;
using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
ComponentVectorMoleFraction x(numComponents), y(numComponents);
// L_eval is converted from a scalar, so all derivatives are zero at this point
InputEval L_eval = fluid_state_scalar.L();;
// for single phase situation, x = y = z;
// and L_eval have all zero derivatives
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);
y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);
}
// set the (trivial) derivatives for x, y and L against P and x.
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
std::vector<double> deriX(num_deri, 0.);
std::vector<double> deriY(num_deri, 0.);
for (unsigned idx = 0; idx < num_deri; ++idx) {
if (idx==0) {
x[compIdx].setDerivative(idx, 0);
y[compIdx].setDerivative(idx, 0);
} else {
if (compIdx==0) {
x[compIdx].setDerivative(1, 1);
y[compIdx].setDerivative(1, 1);
} else {
x[compIdx].setDerivative(1, -1);
y[compIdx].setDerivative(1, -1);
}
}
}
}
std::vector<double> deriL(num_deri, 0.);
for (unsigned idx = 0; idx < num_deri; ++idx) {
L_eval.setDerivative(idx, 0);
}
// update x, y and L in fluid_state
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, z[compIdx]);
fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, z[compIdx]);
}
fluid_state.setLvalue(L_eval);
} //end updateDerivativesSinglePhase

99
tests/test_co2brine_ptflash.cpp
View File

@ -51,7 +51,28 @@ using Evaluation = Opm::DenseAd::Evaluation<double, numComponents>;
typedef Dune::FieldVector<Evaluation, numComponents> ComponentVector;
typedef Opm::CompositionalFluidState<Evaluation, FluidSystem> FluidState;
bool result_okay(const FluidState& fluid_state);
bool result_okay_twophase(const FluidState& fluid_state);
bool result_okay_singlephase(const FluidState& fluid_state, const ComponentVector& z);
bool eval_almost_equal(const Evaluation& val, const Evaluation& ref) {
auto almost_equal = [](const double x, const double y, const double rel_tol = 2.e-3, const double abs_tol = 1.e-3)->bool {
return std::fabs(x - y) <= rel_tol * std::fabs(x + y) * 2 || std::fabs(x - y) < abs_tol;
};
bool equal_okay = true;
if (!almost_equal(val.value(), ref.value())) {
equal_okay = false;
std::cout << " the value are different with " << val.value() << " against the reference " << ref.value() << std::endl;
}
for (int i = 0; i < val.size(); ++i) {
if (!almost_equal(val.derivative(i), ref.derivative(i))) {
equal_okay = false;
std::cout << " the " << i << "th derivative is different with value " << val.derivative(i) << " against the reference " << ref.derivative(i) << std::endl;
}
}
return equal_okay;
};
bool testPTFlash(const std::string& flash_twophase_method)
{
@ -115,7 +136,7 @@ bool testPTFlash(const std::string& flash_twophase_method)
}
const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
const int flash_verbosity = 1;
const int flash_verbosity = 0;
// TODO: should we set these?
// Set initial K and L
@ -130,33 +151,12 @@ bool testPTFlash(const std::string& flash_twophase_method)
using Flash = Opm::PTFlash<double, FluidSystem>;
Flash::solve(fluid_state, z, spatialIdx, flash_twophase_method, flash_tolerance, flash_verbosity);
return result_okay(fluid_state);
return result_okay_twophase(fluid_state);
}
bool result_okay(const FluidState& fluid_state)
bool result_okay_twophase(const FluidState& fluid_state)
{
bool res_okay = true;
auto almost_equal = [](const double x, const double y, const double rel_tol = 2.e-3, const double abs_tol = 1.e-3)->bool {
return std::fabs(x - y) <= rel_tol * std::fabs(x + y) * 2 || std::fabs(x - y) < abs_tol;
};
auto eval_almost_equal = [almost_equal](const Evaluation& val, const Evaluation& ref) -> bool {
bool equal_okay = true;
if (!almost_equal(val.value(), ref.value())) {
equal_okay = false;
std::cout << " the value are different with " << val.value() << " against the reference " << ref.value() << std::endl;
}
for (int i = 0; i < val.size(); ++i) {
if (!almost_equal(val.derivative(i), ref.derivative(i))) {
equal_okay = false;
std::cout << " the " << i << "th derivative is different with value " << val.derivative(i) << " against the reference " << ref.derivative(i) << std::endl;
}
}
return equal_okay;
};
ComponentVector x, y;
const Evaluation L = fluid_state.L();
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
@ -164,7 +164,6 @@ bool result_okay(const FluidState& fluid_state)
y[comp_idx] = fluid_state.moleFraction(FluidSystem::gasPhaseIdx, comp_idx);
}
Evaluation ref_L = 1 - 0.5013878578252918;
ref_L.setDerivative(0, -0.00010420367632860657);
ref_L.setDerivative(1, -1.0043436395393446);
@ -221,9 +220,6 @@ bool testPTFlashSingle(const std::string& flash_twophase_method)
ComponentVector sat;
sat[0] = 1.0; sat[1] = 1.0-sat[0];
// FluidState will be the input for the flash calculation
FluidState fluid_state;
fluid_state.setPressure(FluidSystem::oilPhaseIdx, p_init);
@ -270,8 +266,8 @@ bool testPTFlashSingle(const std::string& flash_twophase_method)
z[numComponents - 1] = z_last;
}
const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
const int flash_verbosity = 1;
const double flash_tolerance = 1.e-12;
const int flash_verbosity = 0;
// TODO: should we set these?
// Set initial K and L
@ -286,9 +282,45 @@ bool testPTFlashSingle(const std::string& flash_twophase_method)
using Flash = Opm::PTFlash<double, FluidSystem>;
Flash::solve(fluid_state, z, spatialIdx, flash_twophase_method, flash_tolerance, flash_verbosity);
return 1;
return result_okay_singlephase(fluid_state, z);
}
//TODO: add when we have something to compare against return result_okay(fluid_state);
bool result_okay_singlephase(const FluidState& fluid_state, const ComponentVector& z)
{
bool res_okay = true;
ComponentVector x, y;
const Evaluation L = fluid_state.L();
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
x[comp_idx] = fluid_state.moleFraction(FluidSystem::oilPhaseIdx, comp_idx);
y[comp_idx] = fluid_state.moleFraction(FluidSystem::gasPhaseIdx, comp_idx);
}
Evaluation ref_L = 1.;
ComponentVector ref_x = z;
ComponentVector ref_y = z;
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
if (!eval_almost_equal(x[comp_idx], ref_x[comp_idx])) {
res_okay = false;
std::cout << " the " << comp_idx << "th x does not match" << std::endl;
}
if (!eval_almost_equal(y[comp_idx], ref_y[comp_idx])) {
res_okay = false;
std::cout << " the " << comp_idx << "th x does not match" << std::endl;
}
}
if (!eval_almost_equal(L, ref_L)) {
res_okay = false;
std::cout << " the L does not match" << std::endl;
}
// TODO: we should also check densities, viscosities, saturations and so on
return res_okay;
}
int main(int argc, char **argv)
@ -307,7 +339,6 @@ int main(int argc, char **argv)
std::cout << method << " solution for PTFlash two-phase case passed " << std::endl;
}
// for the single-phase case (still TODO add refrence and result_okay test)
if (!testPTFlashSingle(method) ) {
std::cout << method << " solution for PTFlash single-phase failed " << std::endl;
test_passed_single = false;