Implemented computeNew() [temporary name], which uses the new vectorized interfaces.
This commit is contained in:
Atgeirr Flø Rasmussen
parent
995d6542bd
commit
a9add0f757
@ -88,6 +88,7 @@ namespace Opm
|
||||
state.drelperm_.resize(num);
|
||||
state.mobility_.resize(num);
|
||||
state.dmobility_.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
FluidMatrixInteractionBlackoil<double>::kr(state.relperm_[i], fmi_params_,
|
||||
state.saturation_[i], state.temperature_[i]);
|
||||
@ -151,8 +152,7 @@ namespace Opm
|
||||
std::vector<double> phasemobf; // Phase mobilities. Flat storage (numPhases per face).
|
||||
std::vector<PhaseVec> phasemobc; // Phase mobilities per cell.
|
||||
std::vector<double> phasemobf_deriv; // Phase mobility derivatives. Flat storage (numPhases^2 per face).
|
||||
typedef Dune::FieldVector<Dune::FieldVector<Scalar, numPhases>, numPhases> PhaseMat;
|
||||
std::vector<PhaseMat> phasemobc_deriv; // Phase mobilities derivatives per cell.
|
||||
std::vector<PhaseJacobian> phasemobc_deriv; // Phase mobilities derivatives per cell.
|
||||
std::vector<double> gravcapf; // Gravity (\rho g \delta z-ish) contribution per face, flat storage.
|
||||
|
||||
public:
|
||||
@ -205,7 +205,7 @@ namespace Opm
|
||||
viscosity[cell] = state.viscosity_;
|
||||
phasemobc[cell] = state.mobility_;
|
||||
phasemobc_deriv[cell] = state.dmobility_;
|
||||
std::copy(state.phase_to_comp_, state.phase_to_comp_ + nc*np, &cellA[cell*nc*np]);
|
||||
std::copy(&state.phase_to_comp_[0][0], &state.phase_to_comp_[0][0] + nc*np, &cellA[cell*nc*np]);
|
||||
// Fractional flow must be calculated.
|
||||
double total_mobility = 0.0;
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
@ -270,7 +270,7 @@ namespace Opm
|
||||
CompVec face_z(0.0);
|
||||
double face_z_factor = 0.5;
|
||||
PhaseVec phase_mob[2];
|
||||
PhaseMat phasemob_deriv[2];
|
||||
PhaseJacobian phasemob_deriv[2];
|
||||
for (int j = 0; j < 2; ++j) {
|
||||
if (c[j] >= 0) {
|
||||
face_z += cell_z[c[j]];
|
||||
@ -313,9 +313,181 @@ namespace Opm
|
||||
|
||||
// Find faceA.
|
||||
FluidStateBlackoil face_state = fluid.computeState(face_pressure[face], face_z);
|
||||
std::copy(face_state.phase_to_comp_, face_state.phase_to_comp_ + nc*np, &faceA[face*nc*np]);
|
||||
std::copy(&face_state.phase_to_comp_[0][0], &face_state.phase_to_comp_[0][0] + nc*np, &faceA[face*nc*np]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
public:
|
||||
template <class Grid, class Rock>
|
||||
void computeNew(const Grid& grid,
|
||||
const Rock& rock,
|
||||
const BlackoilFluid& fluid,
|
||||
const typename Grid::Vector gravity,
|
||||
const std::vector<PhaseVec>& cell_pressure,
|
||||
const std::vector<PhaseVec>& face_pressure,
|
||||
const std::vector<CompVec>& cell_z,
|
||||
const CompVec& bdy_z,
|
||||
const double dt)
|
||||
{
|
||||
int num_cells = cell_z.size();
|
||||
ASSERT(num_cells == grid.numCells());
|
||||
int num_faces = face_pressure.size();
|
||||
ASSERT(num_faces == grid.numFaces());
|
||||
const int np = numPhases;
|
||||
const int nc = numComponents;
|
||||
bool nonzero_gravity = gravity.two_norm() > 0.0;
|
||||
BOOST_STATIC_ASSERT(np == nc);
|
||||
totcompr.resize(num_cells);
|
||||
totphasevol_density.resize(num_cells);
|
||||
voldiscr.resize(num_cells);
|
||||
relvoldiscr.resize(num_cells);
|
||||
saturation.resize(num_cells);
|
||||
frac_flow.resize(num_cells);
|
||||
rel_perm.resize(num_cells);
|
||||
viscosity.resize(num_cells);
|
||||
expjacterm.resize(num_cells);
|
||||
cellA.resize(num_cells*nc*np);
|
||||
faceA.resize(num_faces*nc*np);
|
||||
phasemobf.resize(np*num_faces);
|
||||
phasemobc.resize(num_cells);
|
||||
phasemobf_deriv.resize(np*np*num_faces);
|
||||
phasemobc_deriv.resize(np*np*num_cells);
|
||||
gravcapf.resize(np*num_faces);
|
||||
BOOST_STATIC_ASSERT(np == 3);
|
||||
|
||||
// Compute and copy properties.
|
||||
ManyFluidStatesBlackoil state;
|
||||
fluid.computeManyStates(cell_pressure, cell_z, state);
|
||||
totcompr = state.total_compressibility_;
|
||||
totphasevol_density = state.total_phase_volume_density_;
|
||||
saturation = state.saturation_;
|
||||
rel_perm = state.relperm_;
|
||||
viscosity = state.viscosity_;
|
||||
phasemobc = state.mobility_;
|
||||
phasemobc_deriv = state.dmobility_;
|
||||
std::copy(&state.phase_to_comp_[0][0][0],
|
||||
&state.phase_to_comp_[0][0][0] + num_cells*nc*np,
|
||||
&cellA[0]);
|
||||
expjacterm = state.experimental_term_;
|
||||
|
||||
#pragma omp parallel for
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
double pv = rock.porosity(cell)*grid.cellVolume(cell);
|
||||
voldiscr[cell] = (totphasevol_density[cell] - 1.0)*pv/dt;
|
||||
relvoldiscr[cell] = std::fabs(totphasevol_density[cell] - 1.0);
|
||||
|
||||
// Fractional flow must be calculated.
|
||||
double total_mobility = 0.0;
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
total_mobility += state.mobility_[cell][phase];
|
||||
}
|
||||
frac_flow[cell] = state.mobility_[cell];
|
||||
frac_flow[cell] /= total_mobility;
|
||||
}
|
||||
|
||||
// Obtain properties from both sides of the face.
|
||||
#pragma omp parallel for
|
||||
for (int face = 0; face < num_faces; ++face) {
|
||||
typedef typename Grid::Vector Vec;
|
||||
Vec fc = grid.faceCentroid(face);
|
||||
int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
|
||||
|
||||
// Get pressures and compute gravity contributions,
|
||||
// to decide upwind directions.
|
||||
PhaseVec phase_p[2];
|
||||
PhaseVec gravcontrib[2];
|
||||
for (int j = 0; j < 2; ++j) {
|
||||
if (c[j] >= 0) {
|
||||
// Pressures
|
||||
phase_p[j] = cell_pressure[c[j]];
|
||||
// Gravity contribution.
|
||||
if (nonzero_gravity) {
|
||||
Vec cdiff = fc;
|
||||
cdiff -= grid.cellCentroid(c[j]);
|
||||
gravcontrib[j] = fluid.phaseDensities(&cellA[np*nc*c[j]]);;
|
||||
gravcontrib[j] *= (cdiff*gravity);
|
||||
} else {
|
||||
gravcontrib[j] = 0.0;
|
||||
}
|
||||
} else {
|
||||
// Pressures
|
||||
phase_p[j] = face_pressure[face];
|
||||
// Gravity contribution.
|
||||
gravcontrib[j] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Gravity contribution:
|
||||
// gravcapf = rho_1*g*(z_12 - z_1) - rho_2*g*(z_12 - z_2)
|
||||
// where _1 and _2 refers to two neigbour cells, z is the
|
||||
// z coordinate of the centroid, and z_12 is the face centroid.
|
||||
// Also compute the potentials.
|
||||
PhaseVec pot[2];
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
gravcapf[np*face + phase] = gravcontrib[0][phase] - gravcontrib[1][phase];
|
||||
pot[0][phase] = phase_p[0][phase] + gravcapf[np*face + phase];
|
||||
pot[1][phase] = phase_p[1][phase];
|
||||
}
|
||||
|
||||
// Now we can easily find the upwind direction for every phase,
|
||||
// we can also tell which boundary faces are inflow bdys.
|
||||
|
||||
// Compute face_z, which is averaged from the cells, unless on outflow or noflow bdy.
|
||||
// Get mobilities and derivatives.
|
||||
CompVec face_z(0.0);
|
||||
double face_z_factor = 0.5;
|
||||
PhaseVec phase_mob[2];
|
||||
PhaseJacobian phasemob_deriv[2];
|
||||
for (int j = 0; j < 2; ++j) {
|
||||
if (c[j] >= 0) {
|
||||
face_z += cell_z[c[j]];
|
||||
phase_mob[j] = phasemobc[c[j]];
|
||||
phasemob_deriv[j] = phasemobc_deriv[c[j]];
|
||||
} else if (pot[j][Liquid] > pot[(j+1)%2][Liquid]) {
|
||||
// Inflow boundary.
|
||||
face_z += bdy_z;
|
||||
FluidStateBlackoil bdy_state = fluid.computeState(face_pressure[face], bdy_z);
|
||||
phase_mob[j] = bdy_state.mobility_;
|
||||
phasemob_deriv[j] = bdy_state.dmobility_;
|
||||
} else {
|
||||
// For outflow or noflow boundaries, only cell z is used.
|
||||
face_z_factor = 1.0;
|
||||
// Also, make sure the boundary data are not used for mobilities.
|
||||
pot[j] = -1e100;
|
||||
}
|
||||
}
|
||||
face_z *= face_z_factor;
|
||||
|
||||
// Computing upwind mobilities and derivatives
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (pot[0][phase] == pot[1][phase]) {
|
||||
// Average.
|
||||
double aver = 0.5*(phase_mob[0][phase] + phase_mob[1][phase]);
|
||||
phasemobf[np*face + phase] = aver;
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
phasemobf_deriv[np*np*face + np*phase + p2] = phasemob_deriv[0][phase][p2]
|
||||
+ phasemob_deriv[1][phase][p2];
|
||||
}
|
||||
} else {
|
||||
// Upwind.
|
||||
int upwind = pot[0][phase] > pot[1][phase] ? 0 : 1;
|
||||
phasemobf[np*face + phase] = phase_mob[upwind][phase];
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
phasemobf_deriv[np*np*face + np*phase + p2] = phasemob_deriv[upwind][phase][p2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Find faceA.
|
||||
FluidStateBlackoil face_state = fluid.computeState(face_pressure[face], face_z);
|
||||
std::copy(&face_state.phase_to_comp_[0][0], &face_state.phase_to_comp_[0][0] + nc*np, &faceA[face*nc*np]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user