From 40e2ccd1f0fc0534b89249806ebad20487b1f301 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sun, 26 Feb 2012 20:17:22 +0100 Subject: [PATCH] Removed legacy reorder solver code. --- Makefile.am | 8 +- examples/polymer_reorder.cpp | 23 +-- opm/polymer/polymermodel.cpp | 344 ------------------------------- opm/polymer/polymermodel.hpp | 60 ------ opm/polymer/polymertransport.cpp | 73 ------- opm/polymer/polymertransport.hpp | 30 --- 6 files changed, 4 insertions(+), 534 deletions(-) delete mode 100644 opm/polymer/polymermodel.cpp delete mode 100644 opm/polymer/polymermodel.hpp delete mode 100644 opm/polymer/polymertransport.cpp delete mode 100644 opm/polymer/polymertransport.hpp diff --git a/Makefile.am b/Makefile.am index 055b0c828..675b02c34 100644 --- a/Makefile.am +++ b/Makefile.am @@ -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 diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index eb546e6d7..c9543f491 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -39,8 +39,6 @@ #include // #include -#include -#include #include #include @@ -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 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(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; diff --git a/opm/polymer/polymermodel.cpp b/opm/polymer/polymermodel.cpp deleted file mode 100644 index c7fadfc94..000000000 --- a/opm/polymer/polymermodel.cpp +++ /dev/null @@ -1,344 +0,0 @@ -/* Copyright 2011 (c) Jostein R. Natvig */ -/* Copyright 2012 (c) SINTEF */ - - -#include -#include -#include -#include - -#include -#include -#include - -// #define EXTRA_DEBUG_OUTPUT -#ifdef EXTRA_DEBUG_OUTPUT -#include -#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; inumber_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]; icell_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]; icell_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: */ diff --git a/opm/polymer/polymermodel.hpp b/opm/polymer/polymermodel.hpp deleted file mode 100644 index a56569e72..000000000 --- a/opm/polymer/polymermodel.hpp +++ /dev/null @@ -1,60 +0,0 @@ -/* Copyright 2011 (c) Jostein R. Natvig */ -/* Copyright 2012 (c) SINTEF */ - -#ifndef POLYMER_HPP_INCLUDED -#define POLYMER_HPP_INCLUDED - -#include -#include // 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: */ diff --git a/opm/polymer/polymertransport.cpp b/opm/polymer/polymertransport.cpp deleted file mode 100644 index c0a8380c2..000000000 --- a/opm/polymer/polymertransport.cpp +++ /dev/null @@ -1,73 +0,0 @@ -/* Copyright 2011 (c) Jostein R. Natvig */ -/* Copyright 2012 (c) SINTEF */ - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - - - -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; inumber_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: */ diff --git a/opm/polymer/polymertransport.hpp b/opm/polymer/polymertransport.hpp deleted file mode 100644 index d370c33b7..000000000 --- a/opm/polymer/polymertransport.hpp +++ /dev/null @@ -1,30 +0,0 @@ -/* Copyright 2011 (c) Jostein R. Natvig */ -/* 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 */