diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp index 9d0c97a38..17a7245cf 100644 --- a/opm/material/constraintsolvers/ChiFlash.hpp +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -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 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 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 + template 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 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_ (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] << " " < 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; - 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] << " " < 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 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 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 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 static void evalJacobian(const ComponentVector& globalComposition,