Skip to content

Commit

Permalink
added Kais last to commits and took values to fluid_state
Browse files Browse the repository at this point in the history
  • Loading branch information
Trine Mykkeltvedt committed Jun 8, 2022
1 parent e84435d commit 60564ef
Show file tree
Hide file tree
Showing 2 changed files with 119 additions and 63 deletions.
1 change: 0 additions & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,4 +112,3 @@ opm_add_test(test_fluidsystems)
opm_add_test(test_immiscibleflash)
opm_add_test(test_chiflash)
opm_add_test(test_chiflash_scalar)
opm_add_test(test_chiflashspe5)
181 changes: 119 additions & 62 deletions opm/material/constraintsolvers/ChiFlash.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ class ChiFlash
// TODO: should be able the handle the derivatives directly, which will affect the final organization
// of the code
// TODO: Does fluid_state_scalar contain z with derivatives?
updateDerivatives_(fluid_state_scalar, z_scalar, fluid_state);
updateDerivatives_(fluid_state_scalar, z, fluid_state);

// Update phases
/* typename FluidSystem::template ParameterCache<InputEval> paramCache;
Expand Down Expand Up @@ -869,6 +869,13 @@ class ChiFlash
}
}
}
// for (unsigned i = 0; i < num_equations; ++i) {
// for (unsigned j = 0; j < num_primary_variables; ++j) {
// std::cout << " " << jac[i][j] ;
// }
// std::cout << std::endl;
// }
std::cout << std::endl;
if (!converged) {
throw std::runtime_error(" Newton composition update did not converge within maxIterations");
}
Expand All @@ -881,9 +888,11 @@ class ChiFlash
fluidState.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
const auto K_i = Opm::getValue(flash_fluid_state.K(idx));
fluidState.setKvalue(idx, K_i);
// TODO: not sure we need K and L here, because they are in the flash_fluid_state anyway.
K[idx] = K_i;
}
fluidState.setLvalue(Opm::getValue(flash_fluid_state.L()));
}
L = Opm::getValue(l);
fluidState.setLvalue(L); }

template <class DefectVector>
static void updateCurrentSol_(DefectVector& x, DefectVector& d)
Expand Down Expand Up @@ -955,9 +964,6 @@ class ChiFlash
auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
fluid_state.fugacity(gasPhaseIdx, compIdx));
res[compIdx + numComponents] = Opm::getValue(local_res);
//std::cout << "fugacity oil = " << local_res.derivative(i) << " gas = " << fluid_state.fugacity(gasPhaseIdx, compIdx) << " comp " << compIdx << std::endl; trine


for (unsigned i = 0; i < num_primary; ++i) {
jac[compIdx + numComponents][i] = local_res.derivative(i);
}
Expand All @@ -977,10 +983,10 @@ class ChiFlash
}
}

template <typename FlashFluidStateScalar, typename FlashFluidState, typename ComponentVector>
template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
static void updateDerivatives_(const FlashFluidStateScalar& fluid_state_scalar,
const ComponentVector& z,
FlashFluidState& fluid_state)
FluidState& fluid_state)
{
// getting the secondary Jocobian matrix
constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
Expand All @@ -1001,7 +1007,7 @@ class ChiFlash
secondary_fluid_state.setTemperature(Opm::getValue(fluid_state_scalar.temperature(0)));

for (unsigned idx = 0; idx < numComponents; ++idx) {
secondary_z[idx] = SecondaryEval(z[idx], idx + 1);
secondary_z[idx] = SecondaryEval(Opm::getValue(z[idx]), idx + 1);
}
// set up the mole fractions
for (unsigned idx = 0; idx < num_equations; ++idx) {
Expand Down Expand Up @@ -1051,6 +1057,7 @@ class ChiFlash
primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
}
// TODO: x and y are not needed here
{
std::vector<PrimaryEval> x(numComponents), y(numComponents);
PrimaryEval l;
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
Expand All @@ -1062,8 +1069,8 @@ class ChiFlash
primary_fluid_state.setKvalue(comp_idx, y[comp_idx] / x[comp_idx]);
}
l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
// TODO: do we need to setL?
primary_fluid_state.setLvalue(l);
}
primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));

Expand All @@ -1087,66 +1094,116 @@ class ChiFlash
assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
(primary_fluid_state, primary_z, pri_jac, pri_res);

//corresponds to julias J_s
std::cout << "sec_jac:" << std::endl;
std::cout << "[" << sec_jac[0][0] << " " << sec_jac[0][1] << " " << sec_jac[0][2] << " " << sec_jac[0][3] << "] " << std::endl;
std::cout << "[" << sec_jac[1][0] << " " << sec_jac[1][1] << " " << sec_jac[1][2] << " " << sec_jac[1][3] << "] " << std::endl;
std::cout << "[" << sec_jac[2][0] << " " << sec_jac[2][1] << " " << sec_jac[2][2] << " " << sec_jac[2][3] << "] " << std::endl;
std::cout << "[" << sec_jac[3][0] << " " << sec_jac[3][1] << " " << sec_jac[3][2] << " " << sec_jac[3][3] << "] " << std::endl;
std::cout << "[" << sec_jac[4][0] << " " << sec_jac[4][1] << " " << sec_jac[4][2] << " " << sec_jac[4][3] << "] " << std::endl;
std::cout << "[" << sec_jac[5][0] << " " << sec_jac[5][1] << " " << sec_jac[5][2] << " " << sec_jac[5][3] << "] " << std::endl;
std::cout << "[" << sec_jac[6][0] << " " << sec_jac[6][1] << " " << sec_jac[6][2] << " " << sec_jac[6][3] << "] " << std::endl;

//corresponds to julias J_p (we miss d/dt, and have d/dL instead of d/dV)
std::cout << "pri_jac:" << std::endl;
std::cout << "[" << pri_jac[0][0] << " " << pri_jac[0][1] << " " << pri_jac[0][2] << " " << pri_jac[0][3] << " " << pri_jac[0][4]<< " " << pri_jac[0][5] << " " << pri_jac[0][6]<< "] " << std::endl;
std::cout << "[" << pri_jac[1][0] << " " << pri_jac[1][1] << " " << pri_jac[1][2] << " " << pri_jac[1][3] << " " << pri_jac[1][4]<< " " << pri_jac[1][5] << " " << pri_jac[1][6]<< "] " << std::endl;
std::cout << "[" << pri_jac[2][0] << " " << pri_jac[2][1] << " " << pri_jac[2][2] << " " << pri_jac[2][3] << " " << pri_jac[2][4]<< " " << pri_jac[2][5] << " " << pri_jac[2][6]<< "] " << std::endl;
std::cout << "[" << pri_jac[3][0] << " " << pri_jac[3][1] << " " << pri_jac[3][2] << " " << pri_jac[3][3] << " " << pri_jac[3][4]<< " " << pri_jac[3][5] << " " << pri_jac[3][6]<< "] " << std::endl;
std::cout << "[" << pri_jac[4][0] << " " << pri_jac[4][1] << " " << pri_jac[4][2] << " " << pri_jac[4][3] << " " << pri_jac[4][4]<< " " << pri_jac[4][5] << " " << pri_jac[4][6]<< "] " << std::endl;
std::cout << "[" << pri_jac[5][0] << " " << pri_jac[5][1] << " " << pri_jac[5][2] << " " << pri_jac[5][3] << " " << pri_jac[5][4]<< " " << pri_jac[5][5] << " " << pri_jac[5][6]<< "] " << std::endl;
std::cout << "[" << pri_jac[6][0] << " " << pri_jac[6][1] << " " << pri_jac[6][2] << " " << pri_jac[6][3] << " " << pri_jac[6][4]<< " " << pri_jac[6][5] << " " << pri_jac[6][6]<< "] " << std::endl;
// for (unsigned i =0; i < num_equations; ++i) {
// for (unsigned j = 0; j < primary_num_pv; ++j) {
// std::cout << " " << pri_jac[i][j];
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;

//corresponds to julias J_s
// for (unsigned i = 0; i < num_equations; ++i) {
// for (unsigned j = 0; j < secondary_num_pv; ++j) {
// std::cout << " " << sec_jac[i][j] ;
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;

SecondaryNewtonMatrix xx;
pri_jac.solve(xx,sec_jac);
std::cout << " corresponding to julia-code value and updated J_s " << std::endl;
std::cout << "x1 = [" << x[0] << " " << xx[0][0] << " " << xx[0][1] << " " << xx[0][2] << " " << xx[0][3] << "] " << std::endl;
std::cout << "x2 = [" << x[1] << " " << xx[1][0] << " " <<xx[1][1] << " " << xx[1][2] << " " << xx[1][3] << "] " << std::endl;
std::cout << "x3 = [" << x[2] << " " << xx[2][0] << " " << xx[2][1] << " " << xx[2][2] << " " << xx[2][3] << "] " << std::endl;
std::cout << "y1 = [" << y[0] << " " << xx[3][0] << " " << xx[3][1] << " " << xx[3][2] << " " << xx[3][3] << "] " << std::endl;
std::cout << "y2 = [" << y[1] << " " << xx[4][0] << " " << xx[4][1] << " " << xx[4][2] << " " << xx[4][3] << "] " << std::endl;
std::cout << "y3 = [" << y[2] << " " << xx[5][0] << " " << xx[5][1] << " " << xx[5][2] << " " << xx[5][3] << "] " << std::endl;
std::cout << "L = [" << L << " " << xx[6][0] << " " << xx[6][1] << " " << xx[6][2] << " " << xx[6][3] << "] " << std::endl;

// rewrite like Olav xx --> xxx (do this properly to clean up =)
// z3 = 1 -z2 -z1;
// dx1/dp dx1/dt (dx1/dz1-dx1/dz3) (dx1/dz2 - dx1/dz3)
using TertiaryMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv-1>;
TertiaryMatrix xxx;
xxx[0][0] = xx[0][0];
xxx[1][0] = xx[1][0];
xxx[2][0] = xx[2][0];
xxx[3][0] = xx[0][0];
xxx[4][0] = xx[1][0];
xxx[5][0] = xx[2][0];
xxx[6][0] = xx[3][0];
for (unsigned i = 0; i < primary_num_pv; ++i) { // 7 rekker
xxx[i][1] = Opm::getValue(xx[i][1])-Opm::getValue(xx[i][3]);
xxx[i][2] = Opm::getValue(xx[i][2])-Opm::getValue(xx[i][3]);
}
std::cout << " corresponding to julia-code value and derivatives listed in test_setup NB CHANGE SIGN, and we dont have d/dT " << std::endl;
std::cout << "x1 = [" << x[0] << " " << xxx[0][0] << " " << xxx[0][1] << " " << xxx[0][2] << "] " << std::endl;
std::cout << "x2 = [" << x[1] << " " << xxx[1][0] << " " <<xxx[1][1] << " " << xxx[1][2] << "] " << std::endl;
std::cout << "x3 = [" << x[2] << " " << xxx[2][0] << " " << xxx[2][1] << " " << xxx[2][2] << "] " << std::endl;
std::cout << "y1 = [" << y[0] << " " << xxx[3][0] << " " << xxx[3][1] << " " << xxx[3][2] << "] " << std::endl;
std::cout << "y2 = [" << y[1] << " " << xxx[4][0] << " " << xxx[4][1] << " " << xxx[4][2] << "] " << std::endl;
std::cout << "y3 = [" << y[2] << " " << xxx[5][0] << " " << xxx[5][1] << " " << xxx[5][2] << "] " << std::endl;
std::cout << "L = [" << L << " " << xxx[6][0] << " " << xxx[6][1] << " " << xxx[6][2] << " " << xx[6][3] << "] " << std::endl;

// for (unsigned i = 0; i < num_equations; ++i) {
// for (unsigned j = 0; j < secondary_num_pv; ++j) {
// std::cout << " " << xx[i][j] ;
// }
// std::cout << std::endl;
// }

using InputEval = typename FluidState::Scalar;
std::vector<InputEval> x(numComponents), y(numComponents);
InputEval L_eval = L;
// TODO: then beginning from that point
{
const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
std::vector<double> K(numComponents);

// const double L = fluid_state_scalar.L();
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
K[compIdx] = fluid_state_scalar.K(compIdx);
x[compIdx] = z[compIdx] * 1. / (L + (1 - L) * K[compIdx]);
y[compIdx] = x[compIdx] * K[compIdx];
}
// then we try to set the derivatives for x, y and K against P and x.
// p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary
// pressure joins
{
constexpr size_t num_deri = numComponents;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
std::vector<double> deri(num_deri, 0.);
// derivatives from P
// for (unsigned idx = 0; idx < num_deri; ++idx) {
// probably can use some DUNE operator for vectors or matrics
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx][0] * p_l.derivative(idx);
}
// }

for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
//std::cout << "HEI x[" << compIdx << "] |" << idx << "| " << deri[idx] << " from: " << xx[compIdx][0] << ", " << p_l.derivative(idx) << ", " << pz << ", " << zi << std::endl;
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
x[compIdx].setDerivative(idx, deri[idx]);
}
// handling y
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx + numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
y[compIdx].setDerivative(idx, deri[idx]);
}
}
// handling derivatives of L
std::vector<double> deri(num_deri, 0.);
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[2*numComponents][0] * p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[2*numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
L_eval.setDerivative(idx, deri[idx]);
}
}
}
// x, y og L_eval

}//end updateDerivative
// set up the mole fractions
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.setLvalue(L_eval);
}

/* template <class Vector, class Matrix, class Eval, class ComponentVector>
static void evalJacobian(const ComponentVector& globalComposition,
Expand Down

0 comments on commit 60564ef

Please sign in to comment.