Removed legacy reorder solver code.

This commit is contained in:
Atgeirr Flø Rasmussen 2012-02-26 20:17:22 +01:00
parent c63d817332
commit 40e2ccd1f0
6 changed files with 4 additions and 534 deletions

View File

@ -7,11 +7,7 @@ SUBDIRS = . examples
lib_LTLIBRARIES = libopmpolymer.la
libopmpolymer_la_SOURCES = \
opm/polymer/TransportModelPolymer.cpp \
opm/polymer/polymermodel.cpp \
opm/polymer/polymertransport.cpp
opm/polymer/TransportModelPolymer.cpp
nobase_include_HEADERS = \
opm/polymer/TransportModelPolymer.hpp \
opm/polymer/polymermodel.hpp \
opm/polymer/polymertransport.hpp
opm/polymer/TransportModelPolymer.hpp

View File

@ -39,8 +39,6 @@
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
// #include <opm/core/linalg/LinearSolverIstl.hpp>
#include <opm/polymer/polymertransport.hpp>
#include <opm/polymer/polymermodel.hpp>
#include <opm/polymer/TransportModelPolymer.hpp>
#include <boost/filesystem/convenience.hpp>
@ -319,8 +317,6 @@ main(int argc, char** argv)
double poly_amount = param.getDefault("poly_amount", 5.0);
PolymerInflow poly_inflow(poly_start, poly_end, poly_amount);
bool new_code = param.getDefault("new_code", false);
// Extra rock init.
std::vector<double> porevol;
compute_porevolume(grid->c_grid(), *props, porevol);
@ -432,23 +428,8 @@ main(int argc, char** argv)
// boundary flows must be accumulated into
// source term following the same convention.
transport_timer.start();
if (new_code) {
tmodel.solve(&state.faceflux()[0], &reorder_src[0], stepsize, inflow_c,
&reorder_sat[0], &state.concentration()[0], &state.cmax()[0]);
} else {
polymertransport(&porevol[0],
props->porosity(),
&reorder_src[0],
stepsize,
inflow_c,
const_cast<UnstructuredGrid*>(grid->c_grid()),
props.get(),
&polydata,
&state.faceflux()[0],
&reorder_sat[0],
&state.concentration()[0],
&state.cmax()[0]);
}
tmodel.solve(&state.faceflux()[0], &reorder_src[0], stepsize, inflow_c,
&reorder_sat[0], &state.concentration()[0], &state.cmax()[0]);
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;

View File

@ -1,344 +0,0 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
/* Copyright 2012 (c) SINTEF */
#include <opm/polymer/polymermodel.hpp>
#include <opm/core/grid.h>
#include <opm/core/transport/reorder/nlsolvers.h>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <cstdlib>
#include <cstdio>
#include <cmath>
// #define EXTRA_DEBUG_OUTPUT
#ifdef EXTRA_DEBUG_OUTPUT
#include <iostream>
#endif
/* Parameters used in solution of single-cell boundary-value problem */
struct ParametersSRes
{
double c;
double s0;
double dtpv; /* dt/pv(i) */
double influx; /* sum_j min(v_ij, 0)*f(s_j) */
double outflux; /* sum_j max(v_ij, 0) */
int cell;
const Opm::IncompPropertiesInterface* props;
const Opm::PolymerData* polydata;
};
struct ParametersCRes
{
double s0;
double c0;
double cmax0;
double dtpv; /* dt/pv(i) */
double influx; /* sum_j min(v_ij, 0)*f(s_j) */
double influx_polymer;
double outflux; /* sum_j max(v_ij, 0) */
double porosity;
PolymerSolverData* psdata;
int cell;
double s;
const Opm::PolymerData* polydata;
};
static struct ParametersSRes get_parameters_s(struct PolymerSolverData *d, int cell);
static struct ParametersCRes get_parameters_c(struct PolymerSolverData *d, int cell);
static double residual_s(double s, void *data);
static double residual_c(double c, void *data);
static double fluxfun_props(double s,
double c,
int cell,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* polydata);
static double compute_mc(double c,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* polydata);
void
destroy_solverdata(struct PolymerSolverData *d)
{
if (d!=NULL)
{
free(d->fractionalflow);
free(d->mc);
}
free(d);
}
struct PolymerSolverData *
init_solverdata(struct UnstructuredGrid *grid,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* polydata,
const double *darcyflux,
const double *porevolume,
const double *porosity,
const double *source,
const double dt,
const double inflow_c,
double *saturation,
double *concentration,
double *cmax)
{
int i;
struct PolymerSolverData *d = (struct PolymerSolverData*) malloc(sizeof *d);
if(d!=NULL)
{
d->grid = grid;
d->props = props;
d->polydata = polydata;
d->darcyflux = darcyflux;
d->porevolume = porevolume;
d->porosity = porosity;
d->source = source;
d->dt = dt;
d->inflow_c = inflow_c;
d->saturation = saturation;
d->concentration = concentration;
d->cmax = cmax;
d->fractionalflow = (double*) malloc(grid->number_of_cells *
sizeof *d->fractionalflow);
d->mc = (double*) malloc(grid->number_of_cells *
sizeof *d->mc);
if (d->fractionalflow == NULL || d->mc == NULL)
{
destroy_solverdata(d);
d = NULL;
}
for(i=0; i<grid->number_of_cells; ++i)
{
d->fractionalflow[i] = 0.0;
}
}
return d;
}
/* Solver for single-cell bvp calls root-finder in nlsolvers.c */
void polymer_solvecell(void *data, int cell)
{
struct PolymerSolverData *d = (struct PolymerSolverData*) data;
struct NonlinearSolverCtrl ctrl;
ctrl.maxiterations = 20;
ctrl.nltolerance = 1e-9;
ctrl.method = NonlinearSolverCtrl::REGULAFALSI;
ctrl.initialguess = 0.3*d->polydata->c_max_limit;
ctrl.min_bracket = 0.0;
ctrl.max_bracket = d->polydata->c_max_limit;
struct ParametersCRes prm = get_parameters_c(d, cell);
d->concentration[cell] = find_zero(residual_c, &prm, &ctrl);
d->cmax[cell] = std::max(d->cmax[cell], d->concentration[cell]);
d->saturation[cell] = prm.s;
d->fractionalflow[cell] = fluxfun_props(d->saturation[cell], d->concentration[cell], cell, d->props, d->polydata);
d->mc[cell] = compute_mc(d->concentration[cell], d->props, d->polydata);
}
/* ====================== Internals =================================*/
/* Residual function r(s) for a single-cell bvp */
/*
* r(s) = s - s0 + dt/pv*(influx - outflux*f(s) )
*/
/* influx is water influx, outflux is total outflux */
static double
residual_s(double s, void *data)
{
struct ParametersSRes *p = (struct ParametersSRes*) data;
double c = p->c;
return s - p->s0 + p->dtpv*(p->outflux*fluxfun_props(s, c, p->cell, p->props, p->polydata) + p->influx);
}
static double
residual_c(double c, void *data)
{
struct ParametersCRes *p = (struct ParametersCRes*) data;
int cell = p->cell;
struct ParametersSRes prm_s = get_parameters_s(p->psdata, cell);
prm_s.c = c;
struct NonlinearSolverCtrl ctrl;
ctrl.maxiterations = 20;
ctrl.nltolerance = 1e-9;
ctrl.method = NonlinearSolverCtrl::REGULAFALSI;
ctrl.initialguess = 0.5;
ctrl.min_bracket = 0.2; // TODO: Make this a proper s_min value.
ctrl.max_bracket = 1.0;
double s = find_zero(residual_s, &prm_s, &ctrl);
p->s = s;
double ff = fluxfun_props(s, c, p->cell, prm_s.props, p->polydata);
double mc = compute_mc(c, prm_s.props, p->polydata);
double dps = p->polydata->dps;
double s0 = p->s0;
double c0 = p->c0;
double rhor = p->polydata->rhor;
double porosity = p->porosity;
double cmax0 = p->cmax0;
double ads0 = p->polydata->adsorbtion(std::max(c0, cmax0));
double ads = p->polydata->adsorbtion(std::max(c, cmax0));
double res = (s - dps)*c - (s0 - dps)*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+ p->dtpv*(p->outflux*ff*mc + p->influx_polymer);
#ifdef EXTRA_DEBUG_OUTPUT
std::cout << "c = " << c << " s = " << s << " c-residual = " << res << std::endl;
#endif
return res;
}
static struct ParametersSRes
get_parameters_s(struct PolymerSolverData *d, int cell)
{
int i;
struct UnstructuredGrid *g = d->grid;
struct ParametersSRes p;
double flux;
int f, other;
p.c = d->concentration[cell];
p.s0 = d->saturation[cell];
p.influx = d->source[cell] > 0 ? -d->source[cell] : 0.0;
p.outflux = d->source[cell] <= 0 ? -d->source[cell] : 0.0;
p.dtpv = d->dt/d->porevolume[cell];
p.cell = cell;
p.props = d->props;
p.polydata = d->polydata;
for (i=g->cell_facepos[cell]; i<g->cell_facepos[cell+1]; ++i) {
f = g->cell_faces[i];
/* Compute cell flux*/
if (cell == g->face_cells[2*f]) {
flux = d->darcyflux[f];
other = g->face_cells[2*f+1];
}
else {
flux =-d->darcyflux[f];
other = g->face_cells[2*f];
}
if (other != -1) {
if (flux < 0.0) {
p.influx += flux*d->fractionalflow[other];
}
else {
p.outflux += flux;
}
}
}
return p;
}
static struct ParametersCRes
get_parameters_c(struct PolymerSolverData *d, int cell)
{
int i;
struct UnstructuredGrid *g = d->grid;
double flux;
int f, other;
ParametersCRes p;
p.c0 = d->concentration[cell];
p.cmax0 = d->cmax[cell];
p.s0 = d->saturation[cell];
p.dtpv = d->dt/d->porevolume[cell];
double src = d->source[cell];
p.influx = src > 0 ? -src : 0.0;
p.influx_polymer = src > 0 ? -src*compute_mc(d->inflow_c, d->props, d->polydata) : 0.0;
p.outflux = src <= 0 ? -src : 0.0;
p.porosity = d->porosity[cell];
p.psdata = d;
p.cell = cell;
p.s = -1e100;
p.polydata = d->polydata;
for (i=g->cell_facepos[cell]; i<g->cell_facepos[cell+1]; ++i) {
f = g->cell_faces[i];
/* Compute cell flux*/
if (cell == g->face_cells[2*f]) {
flux = d->darcyflux[f];
other = g->face_cells[2*f+1];
}
else {
flux =-d->darcyflux[f];
other = g->face_cells[2*f];
}
if (other != -1) {
if (flux < 0.0) {
p.influx += flux*d->fractionalflow[other];
p.influx_polymer += flux*d->fractionalflow[other]*d->mc[other];
}
else {
p.outflux += flux;
}
}
}
#ifdef EXTRA_DEBUG_OUTPUT
std::cout << "Cell: " << cell
<< " in: " << p.influx
<< " in(polymer): " << p.influx_polymer
<< " out: " << p.outflux << '\n'
<< " c0: " << p.c0
<< " s0: " << p.s0 << std::endl;
#endif
return p;
}
static double fluxfun_props(double s, double c, int cell,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* pd)
{
const double* visc = props->viscosity();
double c_max_limit = pd->c_max_limit;
double cbar = c/c_max_limit;
double mu_w = visc[0];
double mu_m = pd->viscMult(c)*mu_w;
double mu_p = pd->viscMult(pd->c_max_limit)*mu_w;
double omega = pd->omega;
double mu_m_omega = std::pow(mu_m, omega);
double mu_w_e = mu_m_omega*std::pow(mu_w, 1.0 - omega);
double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - omega);
double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff;
double inv_visc_eff[2] = { inv_mu_w_eff, 1.0/visc[1] };
double sat[2] = { s, 1.0 - s };
double mob[2];
props->relperm(1, sat, &cell, mob, NULL);
mob[0] *= inv_visc_eff[0];
mob[1] *= inv_visc_eff[1];
return mob[0]/(mob[0] + mob[1]);
}
static double compute_mc(double c,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* pd)
{
const double* visc = props->viscosity();
double c_max_limit = pd->c_max_limit;
double cbar = c/c_max_limit;
double mu_w = visc[0];
double mu_m = pd->viscMult(c)*mu_w;
double mu_p = pd->viscMult(pd->c_max_limit)*mu_w;
double omega = pd->omega;
double mu_m_omega = std::pow(mu_m, omega);
double mu_w_e = mu_m_omega*std::pow(mu_w, 1.0 - omega);
double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - omega);
double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff;
return c/(inv_mu_w_eff*mu_p_eff);
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

View File

@ -1,60 +0,0 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
/* Copyright 2012 (c) SINTEF */
#ifndef POLYMER_HPP_INCLUDED
#define POLYMER_HPP_INCLUDED
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/polymer/TransportModelPolymer.hpp> // For PolymerData.
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
}
struct PolymerSolverData {
struct UnstructuredGrid *grid;
const Opm::IncompPropertiesInterface* props;
const Opm::PolymerData* polydata;
const double *darcyflux; /* one flux per face in cdata::grid*/
const double *porevolume; /* one volume per cell */
const double *porosity;
const double *source; /* one source per cell */
double dt;
double inflow_c;
double *saturation; /* one per cell */
double *concentration; /* one per cell */
double *cmax;
double *fractionalflow; /* one per cell */
double *mc;
};
void
polymer_solvecell(void *data, int cell);
void
destroy_solverdata(struct PolymerSolverData *d);
struct PolymerSolverData *
init_solverdata(struct UnstructuredGrid *grid,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* polydata,
const double *darcyflux,
const double *porevolume,
const double *porosity,
const double *source,
const double dt,
const double inflow_c,
double *saturation,
double *concentration,
double *cmax);
#endif /* POLYMER_H_INCLUDED */
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

View File

@ -1,73 +0,0 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
/* Copyright 2012 (c) SINTEF */
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <opm/polymer/polymermodel.hpp>
#include <opm/polymer/polymertransport.hpp>
#include <opm/core/grid.h>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/transport/reorder/tarjan.h>
#include <opm/core/transport/reorder/nlsolvers.h>
void polymertransport(
const double *porevolume,
const double *porosity,
const double *source,
const double dt,
const double inflow_c,
struct UnstructuredGrid *grid,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* polydata,
const double *darcyflux,
double *saturation,
double *concentration,
double *cmax)
{
int i;
int ncomponents;
int *sequence;
int *components;
PolymerSolverData *data = init_solverdata(grid, props, polydata, darcyflux,
porevolume, porosity, source, dt,
inflow_c, saturation,
concentration, cmax);
/* Compute sequence of single-cell problems */
sequence = (int*) malloc( grid->number_of_cells * sizeof *sequence);
components = (int*) malloc(( grid->number_of_cells+1) * sizeof *components);
compute_sequence(grid, darcyflux, sequence, components, &ncomponents);
assert(ncomponents == grid->number_of_cells);
/* Assume all strong components are single-cell domains. */
for(i=0; i<grid->number_of_cells; ++i)
{
#ifdef MATLAB_MEX_FILE
if (interrupt_signal)
{
mexPrintf("Reorder loop interrupted by user: %d of %d "
"cells finished.\n", i, grid->number_of_cells);
break;
}
#endif
polymer_solvecell(data, sequence[i]);
/*print("iterations:%d residual:%20.16f\n",
* ctrl.iterations, ctrl.residual);*/
}
destroy_solverdata(data);
free(sequence);
free(components);
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

View File

@ -1,30 +0,0 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
/* Copyright 2012 (c) SINTEF */
#ifndef POLYMERTRANSPORT_HPP_INCLUDED
#define POLYMERTRANSPORT_HPP_INCLUDED
namespace Opm
{
class IncompPropertiesInterface;
struct PolymerData;
}
struct UnstructuredGrid;
void polymertransport(
const double *porevolume,
const double *porosity,
const double *source,
const double dt,
const double inflow_c,
struct UnstructuredGrid *grid,
const Opm::IncompPropertiesInterface* props,
const Opm::PolymerData* polydata,
const double *darcyflux,
double *saturation,
double *concentration,
double *cmax);
#endif /* POLYMERTRANSPORT_HPP_INCLUDED */