Hi, all.
I am trying to calculate the numerical Jacobian of the reservoir residuals
w.r.t. the initial state. To do this I need to replace the initial solution
of EBOS with a perturbed one. I've been trying for a while, and I am
getting nowhere.
I've attached the code. It is placed at the end of "nonlinearIteration" in
"BlackoilModelEbos.hpp":
if (!report.converged) { ... }
else {
// all the code goes here.
}
I've also made the "nonlinearIteration" take the initial reservoir state as
an input. Thus, we have both the final state and the initial state when the
"else"-case is invoked.
Setting the initial solution as done in the attachment does not work. How
can I pass this information further such that the "linearize()" function
(the one that calculates the residuals and the Jacobian in EBOS) will be
affected?
Let me be more precise. The code runs, and the initial solution is set, but
the "linearize()" function are not aware of the change.
Best regards,
Joakim R. Andersen
else{
ReservoirState final_reservoir_state = reservoir_state;
std::vector<GlobalEqVector> residualsMB(2);
std::vector<Scalar> pert_sizes0 ={-0.000001, -10};
// Calculates the jacobian of the reservoir residuals w.r.t. the intial state.
for(std::size_t cell_block=0; cell_block<9; ++cell_block){ // For each grid block.
for(int state_type = 0; state_type<2; ++state_type){ // For each state_variable (saturation_water, pressure_oil) in that grid block
residualsMB[0] = 0; residualsMB[1] = 0; // reset.
for (int i = 0; i < 2; ++i){ // We are doing central difference
// Copy the states
ReservoirState tmp_initial_reservoir_state = copy_initial_reservoir_state;
ReservoirState tmp_final_reservoir_state = final_reservoir_state ;
WellState tmp_final_well_state = well_state;
dx = 0; // Reset the perturbation vector
// The first iteration calculates the f(x - dx/2)
dx[cell_block][state_type] = (i==0) ? -pert_sizes0[state_type]/2 : pert_sizes0[state_type]/2;
// Perturb the initial state
updateState(dx,tmp_initial_reservoir_state);
// Update the current and the previous solution in ebos with the perturbed initial solution
convertInput( 0, tmp_initial_reservoir_state, ebosSimulator_ );
// Update the current solution in ebos with the final state
convertInput( 1, tmp_final_reservoir_state, ebosSimulator_ );
ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
ebosSimulator_.model().invalidateIntensiveQuantitiesCache(1);
// Calculate the residual and the jacobian
ebosSimulator_.model().linearizer().linearize();
// Need to convert the result to "flow format".
auto& ebosJac = ebosSimulator_.model().linearizer().matrix();
auto& ebosResid = ebosSimulator_.model().linearizer().residual();
convertResults(ebosResid, ebosJac);
//Run the well equations.
double dt = timer.currentStepLength();
wellModel().assemble(ebosSimulator_, iteration, dt, tmp_final_well_state); // This will affect the residuals and the jacobian in ebosSimulator_
// Get a copy of the residuals
auto resMB_perturbed = ebosSimulator_.model().linearizer().residual();
residualsMB[i] = resMB_perturbed;
}
// Calculate the numerical difference
for (std::size_t cell_block_res=0; cell_block_res < 9; ++cell_block_res){
for (std::size_t res_nr=0; res_nr < 2; ++res_nr){
// (f(x+dx/2) - f(x-dx/2)) / (dx)
numA0[cell_block_res][cell_block][res_nr][state_type] = (residualsMB[1][cell_block_res][res_nr] - residualsMB[0][cell_block_res][res_nr])/((-pert_sizes0[state_type]));
}
}
}//end for state_type
} //end for cell_block
} //end else
_______________________________________________
Opm mailing list
[email protected]
http://opm-project.org/cgi-bin/mailman/listinfo/opm
