added Kais last to commits and took values to fluid_state
This commit is contained in:
Trine Mykkeltvedt
parent
83d4194876
commit
783b5011db
@ -183,7 +183,7 @@ public:
|
||||
// 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;
|
||||
@ -869,6 +869,13 @@ protected:
|
||||
}
|
||||
}
|
||||
}
|
||||
// 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");
|
||||
}
|
||||
@ -881,9 +888,11 @@ protected:
|
||||
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)
|
||||
@ -955,9 +964,6 @@ protected:
|
||||
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);
|
||||
}
|
||||
@ -977,10 +983,10 @@ protected:
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
@ -1001,7 +1007,7 @@ protected:
|
||||
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) {
|
||||
@ -1051,6 +1057,7 @@ protected:
|
||||
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) {
|
||||
@ -1062,8 +1069,8 @@ protected:
|
||||
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));
|
||||
|
||||
@ -1087,66 +1094,116 @@ protected:
|
||||
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);
|
||||
|
||||
}//end updateDerivative
|
||||
// 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
|
||||
|
||||
// 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,
|
||||
|
||||
Loading…
Reference in New Issue
Block a user