From 6c34ddd59d4c7a5f2a697535c368b9306e447e28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 24 Aug 2012 16:08:26 +0200 Subject: [PATCH 01/15] Initial commit of tof computation by reordering. --- Makefile.am | 2 + .../reorder/TransportModelTracerTof.cpp | 100 ++++++++++++++++++ .../reorder/TransportModelTracerTof.hpp | 63 +++++++++++ 3 files changed, 165 insertions(+) create mode 100644 opm/core/transport/reorder/TransportModelTracerTof.cpp create mode 100644 opm/core/transport/reorder/TransportModelTracerTof.hpp diff --git a/Makefile.am b/Makefile.am index 712d7c65..547528f6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -98,6 +98,7 @@ opm/core/simulator/SimulatorReport.cpp \ opm/core/simulator/SimulatorTimer.cpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \ opm/core/transport/reorder/TransportModelInterface.cpp \ +opm/core/transport/reorder/TransportModelTracerTof.cpp \ opm/core/transport/reorder/TransportModelTwophase.cpp \ opm/core/transport/reorder/nlsolvers.c \ opm/core/transport/reorder/reordersequence.cpp \ @@ -214,6 +215,7 @@ opm/core/transport/SimpleFluid2pWrapper.hpp \ opm/core/transport/SinglePointUpwindTwoPhase.hpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \ opm/core/transport/reorder/TransportModelInterface.hpp \ +opm/core/transport/reorder/TransportModelTracerTof.hpp \ opm/core/transport/reorder/TransportModelTwophase.hpp \ opm/core/transport/reorder/nlsolvers.h \ opm/core/transport/reorder/reordersequence.h \ diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp new file mode 100644 index 00000000..683b732c --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -0,0 +1,100 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include + +namespace Opm +{ + + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid) + : grid_(grid) + { + } + + + + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in, out] tof Array of time-of-flight values. + void TransportModelTracerTof::solveTof(const double* darcyflux, + const double* porevolume, + std::vector& tof) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + tof.resize(grid_.number_of_cells); + std::fill(tof.begin(), tof.end(), 0.0); + tof_ = &tof[0]; + reorderAndTransport(grid_, darcyflux); + } + + + + + void TransportModelTracerTof::solveSingleCell(const int cell) + { + // Compute flux terms. + double upwind_term = 0.0; + double downwind_flux = 0.0; + for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { + int f = grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == grid_.face_cells[2*f]) { + flux = darcyflux_[f]; + other = grid_.face_cells[2*f+1]; + } else { + flux =-darcyflux_[f]; + other = grid_.face_cells[2*f]; + } + // Add flux to upwind_term or downwind_flux, if interior. + if (other != -1) { + if (flux < 0.0) { + upwind_term += flux*tof_[other]; + } else { + downwind_flux += flux; + } + } + } + + // Compute tof. + tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux; + } + + + + + void TransportModelTracerTof::solveMultiCell(const int /*num_cells*/, const int* /*cells*/) + { + THROW("Not yet implemented."); + } + + + + +} // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp new file mode 100644 index 00000000..1ab7ef3b --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -0,0 +1,63 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED + +#include +#include +#include +#include +struct UnstructuredGrid; + +namespace Opm +{ + + class IncompPropertiesInterface; + + /// Implements a reordering transport solver for incompressible two-phase flow. + class TransportModelTracerTof : public TransportModelInterface + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTof(const UnstructuredGrid& grid); + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in, out] tof Array of time-of-flight values. + void solveTof(const double* darcyflux, + const double* porevolume, + std::vector& tof); + + private: + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + + private: + const UnstructuredGrid& grid_; + const double* darcyflux_; // one flux per grid face + const double* porevolume_; // one volume per cell + double* tof_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED From dfef531fc861448a86df60c65b1a1eb2debf0c63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 09:34:03 +0200 Subject: [PATCH 02/15] Program compute_tof added. --- examples/Makefile.am | 2 + examples/compute_tof.cpp | 241 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 243 insertions(+) create mode 100644 examples/compute_tof.cpp diff --git a/examples/Makefile.am b/examples/Makefile.am index b01f7674..dbc56ae6 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -15,6 +15,7 @@ $(BOOST_SYSTEM_LIB) # Please keep the list sorted. noinst_PROGRAMS = \ +compute_tof \ refine_wells \ scaneclipsedeck \ sim_2p_comp_reorder \ @@ -29,6 +30,7 @@ wells_example # # Please maintain sort order from "noinst_PROGRAMS". +compute_tof_SOURCES = compute_tof.cpp refine_wells_SOURCES = refine_wells.cpp sim_2p_comp_reorder_SOURCES = sim_2p_comp_reorder.cpp sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp new file mode 100644 index 00000000..6ee6d3aa --- /dev/null +++ b/examples/compute_tof.cpp @@ -0,0 +1,241 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + + +namespace +{ + void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) + { + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + } +} // anon namespace + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; + parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr wells; + TwophaseState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + deck.reset(new EclipseGridParser(deck_filename)); + // Grid init + grid.reset(new GridManager(*deck)); + // Rock and fluid init + props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid())); + // Wells init. + wells.reset(new Opm::WellsManager(*deck, *grid->c_grid(), props->permeability())); + // Gravity. + gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state); + } + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Wells init. + wells.reset(new Opm::WellsManager()); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + if (use_gravity) { + if (props->density()[0] == props->density()[1]) { + std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; + } + } + const double *grav = use_gravity ? &gravity[0] : 0; + + // Initialising src + std::vector porevol; + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + int num_cells = grid->c_grid()->number_of_cells; + std::vector src(num_cells, 0.0); + if (use_deck) { + // Do nothing, wells will be the driving force, not source terms. + } else { + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + FlowBCManager bcs; + if (param.getDefault("use_pside", false)) { + int pside = param.get("pside"); + double pside_pressure = param.get("pside_pressure"); + bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + LinearSolverFactory linsolver(param); + + // Pressure solver. + Opm::IncompTpfa psolver(*grid->c_grid(), *props, 0, linsolver, + 0.0, 0.0, 0, + grav, wells->c_wells(), src, bcs.c_bcs()); + + // Tof solver. + Opm::TransportModelTracerTof tofsolver(*grid->c_grid()); + + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::ofstream epoch_os; + std::string output_dir; + if (output) { + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + std::string filename = output_dir + "/epoch_timing.param"; + epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out); + // open file to clean it. The file is appended to in SimulatorTwophase + filename = output_dir + "/step_timing.param"; + std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out); + step_os.close(); + param.writeParam(output_dir + "/simulation.param"); + } + + // Init wells. + Opm::WellState well_state; + well_state.init(wells->c_wells(), state); + + // Main solvers. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::cout << "\n\n================ Starting main solvers ===============" << std::endl; + + // Solve pressure. + pressure_timer.start(); + psolver.solve(1.0, state, well_state); + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + + // Process transport sources (to include bdy terms and well flows). + std::vector transport_src; + Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, + wells->c_wells(), well_state.perfRates(), transport_src); + + // Solve time-of-flight. + std::vector tof(num_cells, 0.0); + transport_timer.start(); + tofsolver.solveTof(&state.faceflux()[0], &porevol[0], tof); + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + total_timer.stop(); + + // Output. + if (output) { + std::string tof_filename = output_dir + "/tof.txt"; + std::ofstream tof_stream(tof_filename.c_str()); + std::copy(tof.begin(), tof.end(), std::ostream_iterator(tof_stream, "\n")); + } + + std::cout << "\n\n================ End of simulation ===============\n" + << "Total time taken: " << total_timer.secsSinceStart() + << "\n Pressure time: " << ptime + << "\n Transport time: " << ttime << std::endl; +} From eeca5abed1dc1a407966a46cbcc6a13ab0614d42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 27 Aug 2012 11:26:51 +0200 Subject: [PATCH 03/15] Implemented rudimentary solveMultiCell(). Simply calls solveSingleCell() once for each cell in block. --- opm/core/transport/reorder/TransportModelTracerTof.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 683b732c..0825706e 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -89,9 +89,12 @@ namespace Opm - void TransportModelTracerTof::solveMultiCell(const int /*num_cells*/, const int* /*cells*/) + void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells) { - THROW("Not yet implemented."); + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + for (int i = 0; i < num_cells; ++i) { + solveSingleCell(cells[i]); + } } From a473422a82df73a8a5ec15b37fd180bf941ddc90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 25 Sep 2012 14:00:17 +0200 Subject: [PATCH 04/15] Add proper support for source terms. This fixes the problem with infinite tofs at sinks. --- examples/compute_tof.cpp | 2 +- .../reorder/TransportModelTracerTof.cpp | 25 ++++++++++++++++--- .../reorder/TransportModelTracerTof.hpp | 9 ++++--- 3 files changed, 28 insertions(+), 8 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 6ee6d3aa..3b10b8c2 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -220,7 +220,7 @@ main(int argc, char** argv) // Solve time-of-flight. std::vector tof(num_cells, 0.0); transport_timer.start(); - tofsolver.solveTof(&state.faceflux()[0], &porevol[0], tof); + tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); transport_timer.stop(); double tt = transport_timer.secsSinceStart(); std::cout << "Transport solver took: " << tt << " seconds." << std::endl; diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 0825706e..70c2fc40 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -21,6 +21,8 @@ #include #include #include +#include +#include namespace Opm { @@ -37,15 +39,25 @@ namespace Opm /// Solve for time-of-flight at next timestep. - /// \param[in] darcyflux Array of signed face fluxes. - /// \param[in] porevolume Array of pore volumes. - /// \param[in, out] tof Array of time-of-flight values. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[out] tof Array of time-of-flight values. void TransportModelTracerTof::solveTof(const double* darcyflux, const double* porevolume, + const double* source, std::vector& tof) { darcyflux_ = darcyflux; porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + THROW("Sources do not sum to zero: " << cum_src); + } +#endif tof.resize(grid_.number_of_cells); std::fill(tof.begin(), tof.end(), 0.0); tof_ = &tof[0]; @@ -58,8 +70,13 @@ namespace Opm void TransportModelTracerTof::solveSingleCell(const int cell) { // Compute flux terms. + // Sources have zero tof, and therefore do not contribute + // to upwind_term. Sinks on the other hand, must be added + // to the downwind_flux (note sign change resulting from + // different sign conventions: pos. source is injection, + // pos. flux is outflow). double upwind_term = 0.0; - double downwind_flux = 0.0; + double downwind_flux = std::max(-source_[cell], 0.0); for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) { int f = grid_.cell_faces[i]; double flux; diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 1ab7ef3b..21375d89 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -40,11 +40,13 @@ namespace Opm TransportModelTracerTof(const UnstructuredGrid& grid); /// Solve for time-of-flight at next timestep. - /// \param[in] darcyflux Array of signed face fluxes. - /// \param[in] porevolume Array of pore volumes. - /// \param[in, out] tof Array of time-of-flight values. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[out] tof Array of time-of-flight values. void solveTof(const double* darcyflux, const double* porevolume, + const double* source, std::vector& tof); private: @@ -55,6 +57,7 @@ namespace Opm const UnstructuredGrid& grid_; const double* darcyflux_; // one flux per grid face const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell double* tof_; }; From 44a4d68816af9b6625f322cacc48772b5b893664 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 08:58:03 +0200 Subject: [PATCH 05/15] Added skeleton of general order DG tof solver. --- Makefile.am | 2 + .../TransportModelTracerTofDiscGal.cpp | 188 ++++++++++++++++++ .../TransportModelTracerTofDiscGal.hpp | 74 +++++++ 3 files changed, 264 insertions(+) create mode 100644 opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp create mode 100644 opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp diff --git a/Makefile.am b/Makefile.am index 4b7657ce..9d13aa1e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -102,6 +102,7 @@ opm/core/simulator/SimulatorTimer.cpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \ opm/core/transport/reorder/TransportModelInterface.cpp \ opm/core/transport/reorder/TransportModelTracerTof.cpp \ +opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp \ opm/core/transport/reorder/TransportModelTwophase.cpp \ opm/core/transport/reorder/nlsolvers.c \ opm/core/transport/reorder/reordersequence.cpp \ @@ -224,6 +225,7 @@ opm/core/transport/SinglePointUpwindTwoPhase.hpp \ opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \ opm/core/transport/reorder/TransportModelInterface.hpp \ opm/core/transport/reorder/TransportModelTracerTof.hpp \ +opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp \ opm/core/transport/reorder/TransportModelTwophase.hpp \ opm/core/transport/reorder/nlsolvers.h \ opm/core/transport/reorder/reordersequence.h \ diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp new file mode 100644 index 00000000..05ba59b0 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -0,0 +1,188 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + + /// A class providing discontinuous Galerkin basis functions. + struct DGBasis + { + static int numBasisFunc(const int dimensions, + const int degree) + { + switch (dimensions) { + case 1: + return degree + 1; + case 2: + return (degree + 2)*(degree + 1)/2; + case 3: + return (degree + 3)*(degree + 2)*(degree + 1)/6; + default: + THROW("Dimensions must be 1, 2 or 3."); + } + } + + /// Evaluate all nonzero basis functions at x, + /// writing to f_x. The array f_x must have + /// size numBasisFunc(grid.dimensions, degree). + /// + /// The basis functions are the following + /// Degree 0: 1. + /// Degree 1: x - xc, y - yc, z - zc etc. + /// Further degrees await development. + static void eval(const UnstructuredGrid& grid, + const int cell, + const int degree, + const double* x, + double* f_x) + { + const int dim = grid.dimensions; + const double* cc = grid.cell_centroids + dim*cell; + // Note intentional fallthrough in this switch statement! + switch (degree) { + case 1: + for (int ix = 0; ix < dim; ++ix) { + f_x[1 + ix] = x[ix] - cc[ix]; + } + case 0: + f_x[0] = 1; + break; + default: + THROW("Maximum degree is 1 for now."); + } + } + + /// Evaluate gradients of all nonzero basis functions at x, + /// writing to grad_f_x. The array grad_f_x must have size + /// numBasisFunc(grid.dimensions, degree) * grid.dimensions. + /// The components of the first basis function + /// gradient come before the components of the second etc. + static void evalGrad(const UnstructuredGrid& grid, + const int /*cell*/, + const int degree, + const double* /*x*/, + double* grad_f_x) + { + const int dim = grid.dimensions; + const int num_basis = numBasisFunc(dim, degree); + std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0); + if (degree > 1) { + THROW("Maximum degree is 1 for now."); + } else if (degree == 1) { + for (int ix = 0; ix < dim; ++ix) { + grad_f_x[dim*(ix + 1) + ix] = 1.0; + } + } + } + }; + + + + + /// A class providing numerical quadrature functions. + struct Quadrature + { + static void x(){} + }; + + + + + + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid) + : grid_(grid) + { + } + + + + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[in] degree Polynomial degree of DG basis functions used. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + const int degree, + std::vector& tof_coeff) + { + darcyflux_ = darcyflux; + porevolume_ = porevolume; + source_ = source; +#ifndef NDEBUG + // Sanity check for sources. + const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); + if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { + THROW("Sources do not sum to zero: " << cum_src); + } +#endif + degree_ = degree; + tof_coeff.resize(grid_.number_of_cells); + std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); + tof_coeff_ = &tof_coeff[0]; + reorderAndTransport(grid_, darcyflux); + } + + + + + void TransportModelTracerTofDiscGal::solveSingleCell(const int cell) + { + // Residual: + // For each cell K, basis function b_j (spanning V_h), + // writing the solution u_h|K = \sum_i c_i b_i + // Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx + // + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds + // - \int_K \phi b_j + THROW("Not implemented yet!"); + } + + + + + void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells) + { + std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl; + for (int i = 0; i < num_cells; ++i) { + solveSingleCell(cells[i]); + } + } + + + + +} // namespace Opm diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp new file mode 100644 index 00000000..94d5edf1 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -0,0 +1,74 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED + +#include +#include +#include +#include +struct UnstructuredGrid; + +namespace Opm +{ + + class IncompPropertiesInterface; + + /// Implements a reordering transport solver for incompressible two-phase flow. + class TransportModelTracerTofDiscGal : public TransportModelInterface + { + public: + /// Construct solver. + /// \param[in] grid A 2d or 3d grid. + TransportModelTracerTofDiscGal(const UnstructuredGrid& grid); + + /// Solve for time-of-flight at next timestep. + /// \param[in] darcyflux Array of signed face fluxes. + /// \param[in] porevolume Array of pore volumes. + /// \param[in] source Transport source term. + /// \param[in] degree Polynomial degree of DG basis functions used. + /// \param[out] tof_coeff Array of time-of-flight solution coefficients. + /// The values are ordered by cell, meaning that + /// the K coefficients corresponding to the first + /// cell comes before the K coefficients corresponding + /// to the second cell etc. + /// K depends on degree and grid dimension. + void solveTof(const double* darcyflux, + const double* porevolume, + const double* source, + const int degree, + std::vector& tof_coeff); + + private: + virtual void solveSingleCell(const int cell); + virtual void solveMultiCell(const int num_cells, const int* cells); + + private: + const UnstructuredGrid& grid_; + const double* darcyflux_; // one flux per grid face + const double* porevolume_; // one volume per cell + const double* source_; // one volumetric source term per cell + int degree_; + double* tof_coeff_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED From 3ba50eb2d100dd7de58f091d7056dc39f21f4b46 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 13:30:54 +0200 Subject: [PATCH 06/15] Initial version of DG(1) for tof implemented. Basis functions, quadratures and velocity interpolation are basic versions, not handling any higher than DG(1) for now. These are currently in helper classes and functions. The code in the main solver class is written with the aim of supporting DG(n) generally. --- .../TransportModelTracerTofDiscGal.cpp | 228 +++++++++++++++++- .../TransportModelTracerTofDiscGal.hpp | 8 + 2 files changed, 230 insertions(+), 6 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 05ba59b0..e626a2b3 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -103,20 +104,115 @@ namespace Opm - /// A class providing numerical quadrature functions. - struct Quadrature + /// A class providing numerical quadrature for cells. + class CellQuadrature { - static void x(){} + public: + CellQuadrature(const UnstructuredGrid& grid, + const int cell, + const int degree) + : grid_(grid), cell_(cell), degree_(degree) + { + if (degree > 1) { + THROW("Only quadrature degree up to 1 for now."); + } + } + + int numQuadPts() const + { + return 1; + } + + void quadPtCoord(int /*index*/, double* coord) const + { + const double* cc = grid_.cell_centroids + grid_.dimensions*cell_; + std::copy(cc, cc + grid_.dimensions, coord); + } + + double quadPtWeight(int /*index*/) const + { + return 1.0; + } + + private: + const UnstructuredGrid& grid_; + const int cell_; + const int degree_; }; + /// A class providing numerical quadrature for faces. + class FaceQuadrature + { + public: + FaceQuadrature(const UnstructuredGrid& grid, + const int face, + const int degree) + : grid_(grid), face_(face), degree_(degree) + { + if (degree > 1) { + THROW("Only quadrature degree up to 1 for now."); + } + } + + int numQuadPts() const + { + return 1; + } + + void quadPtCoord(int /*index*/, double* coord) const + { + const double* fc = grid_.face_centroids + grid_.dimensions*face_; + std::copy(fc, fc + grid_.dimensions, coord); + } + + double quadPtWeight(int /*index*/) const + { + return 1.0; + } + + private: + const UnstructuredGrid& grid_; + const int face_; + const int degree_; + }; + + + // Initial version: only a constant interpolation. + static void interpolateVelocity(const UnstructuredGrid& grid, + const int cell, + const double* darcyflux, + const double* /*x*/, + double* v) + { + const int dim = grid.dimensions; + std::fill(v, v + dim, 0.0); + const double* cc = grid.cell_centroids + cell*dim; + for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) { + const int face = grid.cell_faces[hface]; + const double* fc = grid.face_centroids + face*dim; + double flux = 0.0; + if (cell == grid.face_cells[2*face]) { + flux = darcyflux[face]; + } else { + flux = -darcyflux[face]; + } + for (int dd = 0; dd < dim; ++dd) { + v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell]; + } + } + } + + /// Construct solver. /// \param[in] grid A 2d or 3d grid. TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid) - : grid_(grid) + : grid_(grid), + coord_(grid.dimensions), + velocity_(grid.dimensions) { } @@ -151,9 +247,15 @@ namespace Opm } #endif degree_ = degree; - tof_coeff.resize(grid_.number_of_cells); + const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_); + tof_coeff.resize(num_basis*grid_.number_of_cells); std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0); tof_coeff_ = &tof_coeff[0]; + rhs_.resize(num_basis); + jac_.resize(num_basis*num_basis); + basis_.resize(num_basis); + basis_nb_.resize(num_basis); + grad_basis_.resize(num_basis*grid_.dimensions); reorderAndTransport(grid_, darcyflux); } @@ -168,7 +270,121 @@ namespace Opm // Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx // + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds // - \int_K \phi b_j - THROW("Not implemented yet!"); + // This is linear in c_i, so we do not need any nonlinear iterations. + // We assemble the jacobian and the right-hand side. The residual is + // equal to Res = Jac*c - rhs, and we compute rhs directly. + + const int dim = grid_.dimensions; + const int num_basis = DGBasis::numBasisFunc(degree_, dim); + + std::fill(rhs_.begin(), rhs_.end(), 0.0); + std::fill(jac_.begin(), jac_.end(), 0.0); + + // Compute cell residual contribution. + // Note: Assumes that \int_K b_j = 0 for all j > 0 + rhs_[0] += porevolume_[cell]; + + // Compute upstream residual contribution. + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + int upstream_cell = -1; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face+1]; + } else { + flux = -darcyflux_[face]; + upstream_cell = grid_.face_cells[2*face]; + } + if (upstream_cell < 0) { + // This is an outer boundary. Assumed tof = 0 on inflow, so no contribution. + continue; + } + if (flux >= 0.0) { + // This is an outflow boundary. + continue; + } + // Do quadrature over the face to compute + // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds + // (where u_h^{ext} is the upstream unknown (tof)). + FaceQuadrature quad(grid_, face, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // u^ext flux B (B = {b_j}) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]); + const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(), + tof_coeff_ + num_basis*upstream_cell, 0.0); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + rhs_[j] -= w * tof_upstream * flux * basis_[j]; + } + } + } + + // Compute cell jacobian contribution. We use Fortran ordering + // for jac_, i.e. rows cycling fastest. + { + CellQuadrature quad(grid_, cell, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // b_i (v \cdot \grad b_j) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]); + interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + for (int dd = 0; dd < dim; ++dd) { + jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd]; + } + } + } + } + } + + // Compute downstream jacobian contribution. + for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { + const int face = grid_.cell_faces[hface]; + double flux = 0.0; + if (cell == grid_.face_cells[2*face]) { + flux = darcyflux_[face]; + } else { + flux = -darcyflux_[face]; + } + if (flux <= 0.0) { + // This is an inflow boundary. + continue; + } + // Do quadrature over the face to compute + // \int_{\partial K} b_i (v(x) \cdot n) b_j ds + FaceQuadrature quad(grid_, face, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + // u^ext flux B (B = {b_j}) + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + } + } + } + } + + // Solve linear equation. + MAT_SIZE_T n = num_basis; + MAT_SIZE_T nrhs = 1; + MAT_SIZE_T lda = num_basis; + std::vector piv(num_basis); + MAT_SIZE_T ldb = num_basis; + MAT_SIZE_T info = 0; + dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); + if (info != 0) { + THROW("Lapack error: " << info); + } + // The solution ends up in rhs_, so we must copy it. + std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); } diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index 94d5edf1..b9163f8e 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -67,6 +67,14 @@ namespace Opm const double* source_; // one volumetric source term per cell int degree_; double* tof_coeff_; + std::vector rhs_; // single-cell right-hand-side + std::vector jac_; // single-cell jacobian + // Below: storage for quantities needed by solveSingleCell(). + std::vector coord_; + std::vector basis_; + std::vector basis_nb_; + std::vector grad_basis_; + std::vector velocity_; }; } // namespace Opm From 713bb630451d61662b3154c00018dd14ef02df8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 14:21:00 +0200 Subject: [PATCH 07/15] Fix argument order in call. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index e626a2b3..c042614a 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -275,7 +275,7 @@ namespace Opm // equal to Res = Jac*c - rhs, and we compute rhs directly. const int dim = grid_.dimensions; - const int num_basis = DGBasis::numBasisFunc(degree_, dim); + const int num_basis = DGBasis::numBasisFunc(dim, degree_); std::fill(rhs_.begin(), rhs_.end(), 0.0); std::fill(jac_.begin(), jac_.end(), 0.0); From 74b74bca99c7574f09253cfb7ee07d73bcf9689d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 14:21:16 +0200 Subject: [PATCH 08/15] Make DG(n) solver an option for test program. For now, only DG(0,1) is expected to work. --- examples/compute_tof.cpp | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index 3b10b8c2..ea55b172 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include @@ -165,8 +166,12 @@ main(int argc, char** argv) 0.0, 0.0, 0, grav, wells->c_wells(), src, bcs.c_bcs()); - // Tof solver. - Opm::TransportModelTracerTof tofsolver(*grid->c_grid()); + // Choice of tof solver. + bool use_dg = param.getDefault("use_dg", false); + int dg_degree = -1; + if (use_dg) { + dg_degree = param.getDefault("dg_degree", 0); + } // Write parameters used for later reference. bool output = param.getDefault("output", true); @@ -218,10 +223,18 @@ main(int argc, char** argv) wells->c_wells(), well_state.perfRates(), transport_src); // Solve time-of-flight. - std::vector tof(num_cells, 0.0); - transport_timer.start(); - tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); - transport_timer.stop(); + std::vector tof; + if (use_dg) { + Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid()); + transport_timer.start(); + tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof); + transport_timer.stop(); + } else { + Opm::TransportModelTracerTof tofsolver(*grid->c_grid()); + transport_timer.start(); + tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof); + transport_timer.stop(); + } double tt = transport_timer.secsSinceStart(); std::cout << "Transport solver took: " << tt << " seconds." << std::endl; ttime += tt; From 8cfb8f005e3643d88b1d14028de40e6c7797bc2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 26 Sep 2012 15:15:04 +0200 Subject: [PATCH 09/15] Add sink term contribution. --- .../TransportModelTracerTofDiscGal.cpp | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index c042614a..22e60a36 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -343,7 +343,7 @@ namespace Opm } } - // Compute downstream jacobian contribution. + // Compute downstream jacobian contribution from faces. for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) { const int face = grid_.cell_faces[hface]; double flux = 0.0; @@ -372,6 +372,25 @@ namespace Opm } } + // Compute downstream jacobian contribution from sink terms. + if (source_[cell] < 0.0) { + // A sink. + const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. + // Do quadrature over the cell to compute + // \int_{K} b_i flux b_j dx + CellQuadrature quad(grid_, cell, degree_); + for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { + quad.quadPtCoord(quad_pt, &coord_[0]); + DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); + const double w = quad.quadPtWeight(quad_pt); + for (int j = 0; j < num_basis; ++j) { + for (int i = 0; i < num_basis; ++i) { + jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + } + } + } + } + // Solve linear equation. MAT_SIZE_T n = num_basis; MAT_SIZE_T nrhs = 1; From cf2b975d459c5af3c9cf5a85171fca77d9cb90c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 27 Sep 2012 09:49:36 +0200 Subject: [PATCH 10/15] Fix: forgotten multiply by cell volume in a quadrature. --- opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 22e60a36..cdf909b3 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -325,6 +325,7 @@ namespace Opm // Compute cell jacobian contribution. We use Fortran ordering // for jac_, i.e. rows cycling fastest. { + const double cvol = grid_.cell_volumes[cell]; CellQuadrature quad(grid_, cell, degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // b_i (v \cdot \grad b_j) @@ -336,7 +337,7 @@ namespace Opm for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { for (int dd = 0; dd < dim; ++dd) { - jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd]; + jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd] * cvol; } } } From 756ab674efd5980fd4ce218d9f879326456b5716 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 28 Sep 2012 14:44:04 +0200 Subject: [PATCH 11/15] Work in progress: degree 2 quadratures. Also, changed quadrature degrees used to get exact quadratures for all terms. --- .../TransportModelTracerTofDiscGal.cpp | 282 ++++++++++++++++-- 1 file changed, 257 insertions(+), 25 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index cdf909b3..e9e812ef 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -29,6 +29,10 @@ namespace Opm { + // --------------- Helpers for TransportModelTracerTofDiscGal --------------- + + + /// A class providing discontinuous Galerkin basis functions. struct DGBasis { @@ -103,8 +107,83 @@ namespace Opm + static void cross(const double* a, const double* b, double* res) + { + res[0] = a[1]*b[2] - a[2]*b[1]; + res[1] = a[2]*b[0] - a[0]*b[2]; + res[2] = a[0]*b[1] - a[1]*b[0]; + } + + + + + static double triangleArea3d(const double* p0, + const double* p1, + const double* p2) + { + double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; + double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; + double cr[3]; + cross(a, b, cr); + return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]); + } + + + + + /// Calculates the determinant of a 3 x 3 matrix, represented as + /// three three-dimensional arrays. + static double determinantOf(const double* a0, + const double* a1, + const double* a2) + { + return + a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - + a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + + a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); + } + + + + + /// Computes the volume of a tetrahedron consisting of 4 vertices + /// with 3-dimensional coordinates + static double tetVolume(const double* p0, + const double* p1, + const double* p2, + const double* p3) + { + double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; + double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; + double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] }; + return std::fabs(determinantOf(a, b, c) / 6.0); + } + + + /// A class providing numerical quadrature for cells. + /// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). + /// Note that this class does multiply weights by cell volume, + /// so weights always sum to cell volume. + /// Degree 1 method: + /// Midpoint (centroid) method. + /// n = 1, w_0 = cell volume, x_0 = cell centroid + /// Degree 2 method: + /// Based on subdivision of each cell face into triangles + /// with the face centroid as a common vertex, and then + /// subdividing the cell into tetrahedra with the cell + /// centroid as a common vertex. Then apply the tetrahedron + /// rule with the following 4 nodes (uniform weights): + /// a = 0.138196601125010515179541316563436 + /// x_i has all barycentric coordinates = a, except for + /// the i'th coordinate which is = 1 - 3a. + /// This rule is from http://nines.cs.kuleuven.be/ecf, + /// it is the second degree 2 4-point rule for tets, + /// referenced to Stroud(1971). + /// The tetrahedra are numbered T_{i,j}, and are given by the + /// cell centroid C, the face centroid FC_i, and two nodes + /// of face i: FN_{i,j}, FN_{i,j+1}. class CellQuadrature { public: @@ -113,25 +192,93 @@ namespace Opm const int degree) : grid_(grid), cell_(cell), degree_(degree) { - if (degree > 1) { - THROW("Only quadrature degree up to 1 for now."); + if (degree > 2) { + THROW("CellQuadrature exact for polynomial degrees > 1 not implemented."); + } + if (degree == 2) { + // Prepare subdivision. } } int numQuadPts() const { - return 1; + if (degree_ < 2) { + return 1; + } + // Degree 2 case. + int sumnodes = 0; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + } + return 4*sumnodes; } - void quadPtCoord(int /*index*/, double* coord) const + void quadPtCoord(const int index, double* coord) const { - const double* cc = grid_.cell_centroids + grid_.dimensions*cell_; - std::copy(cc, cc + grid_.dimensions, coord); + const int dim = grid_.dimensions; + const double* cc = grid_.cell_centroids + dim*cell_; + if (degree_ < 2) { + std::copy(cc, cc + dim, coord); + return; + } + // Degree 2 case. + int tetindex = index / 4; + const int subindex = index % 4; + const double* nc = grid_.node_coordinates; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + if (nfn <= tetindex) { + // Our tet is not associated with this face. + tetindex -= nfn; + continue; + } + const double* fc = grid_.face_centroids + dim*face; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; + const int node0 = fnodes[tetindex]; + const int node1 = fnodes[(tetindex + 1) % nfn]; + const double* n0c = nc + dim*node0; + const double* n1c = nc + dim*node1; + const double a = 0.138196601125010515179541316563436; + // Barycentric coordinates of our point in the tet. + double baryc[4] = { a, a, a, a }; + baryc[subindex] = 1.0 - 3.0*a; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd]; + } + return; + } + THROW("Should never reach this point."); } - double quadPtWeight(int /*index*/) const + double quadPtWeight(const int index) const { - return 1.0; + if (degree_ < 2) { + return grid_.cell_volumes[cell_]; + } + // Degree 2 case. + const int dim = grid_.dimensions; + const double* cc = grid_.cell_centroids + dim*cell_; + int tetindex = index / 4; + const double* nc = grid_.node_coordinates; + for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) { + const int face = grid_.cell_faces[hf]; + const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face]; + if (nfn <= tetindex) { + // Our tet is not associated with this face. + tetindex -= nfn; + continue; + } + const double* fc = grid_.face_centroids + dim*face; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face]; + const int node0 = fnodes[tetindex]; + const int node1 = fnodes[(tetindex + 1) % nfn]; + const double* n0c = nc + dim*node0; + const double* n1c = nc + dim*node1; + return 0.25*tetVolume(cc, fc, n0c, n1c); + } + THROW("Should never reach this point."); } private: @@ -142,7 +289,28 @@ namespace Opm + + /// A class providing numerical quadrature for faces. + /// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i). + /// Note that this class does multiply weights by face area, + /// so weights always sum to face area. + /// Degree 1 method: + /// Midpoint (centroid) method. + /// n = 1, w_0 = face area, x_0 = face centroid + /// Degree 2 method: + /// Based on subdivision of the face into triangles, + /// with the centroid as a common vertex, and the triangle + /// edge midpoint rule. + /// Triangle i consists of the centroid C, nodes N_i and N_{i+1}. + /// Its area is A_i. + /// n = 2 * nn (nn = num nodes in face) + /// For i = 0..(nn-1): + /// w_i = 1/3 A_i. + /// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i + /// x_i = (N_i + N_{i+1})/2 + /// x_{nn+i} = (C + N_i)/2 + /// All N and A indices are interpreted cyclic, modulus nn. class FaceQuadrature { public: @@ -151,25 +319,79 @@ namespace Opm const int degree) : grid_(grid), face_(face), degree_(degree) { - if (degree > 1) { - THROW("Only quadrature degree up to 1 for now."); + if (grid_.dimensions != 3) { + THROW("FaceQuadrature only implemented for 3D case."); + } + if (degree_ > 2) { + THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented."); } } int numQuadPts() const { - return 1; + if (degree_ < 2) { + return 1; + } + // Degree 2 case. + return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]); } - void quadPtCoord(int /*index*/, double* coord) const + void quadPtCoord(const int index, double* coord) const { - const double* fc = grid_.face_centroids + grid_.dimensions*face_; - std::copy(fc, fc + grid_.dimensions, coord); + const int dim = grid_.dimensions; + const double* fc = grid_.face_centroids + dim*face_; + if (degree_ < 2) { + std::copy(fc, fc + dim, coord); + return; + } + // Degree 2 case. + const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; + const double* nc = grid_.node_coordinates; + if (index < nn) { + // Boundary edge midpoint. + const int node0 = fnodes[index]; + const int node1 = fnodes[(index + 1)%nn]; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]); + } + } else { + // Interiour edge midpoint. + // Recall that index is now in [nn, 2*nn). + const int node = fnodes[index - nn]; + for (int dd = 0; dd < dim; ++dd) { + coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]); + } + } } - double quadPtWeight(int /*index*/) const + double quadPtWeight(const int index) const { - return 1.0; + if (degree_ < 2) { + return grid_.face_areas[face_]; + } + // Degree 2 case. + const int dim = grid_.dimensions; + const double* fc = grid_.face_centroids + dim*face_; + const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]; + const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_]; + const double* nc = grid_.node_coordinates; + if (index < nn) { + // Boundary edge midpoint. + const int node0 = fnodes[index]; + const int node1 = fnodes[(index + 1)%nn]; + const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); + return area / 3.0; + } else { + // Interiour edge midpoint. + // Recall that index is now in [nn, 2*nn). + const int node0 = fnodes[(index - 1) % nn]; + const int node1 = fnodes[index - nn]; + const int node2 = fnodes[(index + 1) % nn]; + const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc); + const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc); + return (area0 + area1) / 3.0; + } } private: @@ -179,6 +401,8 @@ namespace Opm }; + + // Initial version: only a constant interpolation. static void interpolateVelocity(const UnstructuredGrid& grid, const int cell, @@ -206,6 +430,9 @@ namespace Opm + // --------------- Methods of TransportModelTracerTofDiscGal --------------- + + /// Construct solver. /// \param[in] grid A 2d or 3d grid. @@ -307,9 +534,9 @@ namespace Opm // Do quadrature over the face to compute // \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds // (where u_h^{ext} is the upstream unknown (tof)). + const double normal_velocity = flux / grid_.face_areas[face]; FaceQuadrature quad(grid_, face, degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { - // u^ext flux B (B = {b_j}) quad.quadPtCoord(quad_pt, &coord_[0]); DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]); @@ -317,7 +544,7 @@ namespace Opm tof_coeff_ + num_basis*upstream_cell, 0.0); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { - rhs_[j] -= w * tof_upstream * flux * basis_[j]; + rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j]; } } } @@ -325,8 +552,7 @@ namespace Opm // Compute cell jacobian contribution. We use Fortran ordering // for jac_, i.e. rows cycling fastest. { - const double cvol = grid_.cell_volumes[cell]; - CellQuadrature quad(grid_, cell, degree_); + CellQuadrature quad(grid_, cell, 2*degree_ - 1); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // b_i (v \cdot \grad b_j) quad.quadPtCoord(quad_pt, &coord_[0]); @@ -337,7 +563,7 @@ namespace Opm for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { for (int dd = 0; dd < dim; ++dd) { - jac_[j*num_basis + i] += w * basis_[i] * grad_basis_[dim*j + dd] * velocity_[dd] * cvol; + jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd]; } } } @@ -359,7 +585,8 @@ namespace Opm } // Do quadrature over the face to compute // \int_{\partial K} b_i (v(x) \cdot n) b_j ds - FaceQuadrature quad(grid_, face, degree_); + const double normal_velocity = flux / grid_.face_areas[face]; + FaceQuadrature quad(grid_, face, 2*degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { // u^ext flux B (B = {b_j}) quad.quadPtCoord(quad_pt, &coord_[0]); @@ -367,26 +594,31 @@ namespace Opm const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { - jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j]; } } } } // Compute downstream jacobian contribution from sink terms. + // Contribution from inflow sources would be + // similar to the contribution from upstream faces, but + // it is zero since we let all external inflow be associated + // with a zero tof. if (source_[cell] < 0.0) { // A sink. const double flux = -source_[cell]; // Sign convention for flux: outflux > 0. + const double flux_density = flux / grid_.cell_volumes[cell]; // Do quadrature over the cell to compute // \int_{K} b_i flux b_j dx - CellQuadrature quad(grid_, cell, degree_); + CellQuadrature quad(grid_, cell, 2*degree_); for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) { quad.quadPtCoord(quad_pt, &coord_[0]); DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]); const double w = quad.quadPtWeight(quad_pt); for (int j = 0; j < num_basis; ++j) { for (int i = 0; i < num_basis; ++i) { - jac_[j*num_basis + i] += w * basis_[i] * flux * basis_[j]; + jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j]; } } } From d1be859828e403be9a5443a02ef839ebc1f49b81 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 09:54:26 +0200 Subject: [PATCH 12/15] Add timing of topological sort. --- opm/core/transport/reorder/TransportModelInterface.cpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index f8e138a3..8a90ae34 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g std::vector sequence(grid.number_of_cells); std::vector components(grid.number_of_cells + 1); int ncomponents; + time::StopWatch clock; + clock.start(); compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents); + clock.stop(); + std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl; // Invoke appropriate solve method for each interdependent component. for (int comp = 0; comp < ncomponents; ++comp) { From f761e73eaef927b84bffc3126c5ed419e43a66ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 09:54:54 +0200 Subject: [PATCH 13/15] Improve docs and give more info on error. --- .../reorder/TransportModelTracerTof.cpp | 6 ++++-- .../reorder/TransportModelTracerTof.hpp | 14 ++++++++++--- .../TransportModelTracerTofDiscGal.cpp | 21 ++++++++++++++++--- .../TransportModelTracerTofDiscGal.hpp | 17 ++++++++++++--- 4 files changed, 47 insertions(+), 11 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp index 70c2fc40..8d93f248 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -38,10 +38,12 @@ namespace Opm - /// Solve for time-of-flight at next timestep. + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[out] tof Array of time-of-flight values. void TransportModelTracerTof::solveTof(const double* darcyflux, const double* porevolume, diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp index 21375d89..270ac328 100644 --- a/opm/core/transport/reorder/TransportModelTracerTof.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -31,7 +31,13 @@ namespace Opm class IncompPropertiesInterface; - /// Implements a reordering transport solver for incompressible two-phase flow. + /// Implements a first-order finite volume solver for + /// (single-phase) time-of-flight using reordering. + /// The equation solved is: + /// v \cdot \grad\tau = \phi + /// where v is the fluid velocity, \tau is time-of-flight and + /// \phi is the porosity. This is a boundary value problem, where + /// \tau is specified to be zero on all inflow boundaries. class TransportModelTracerTof : public TransportModelInterface { public: @@ -39,10 +45,12 @@ namespace Opm /// \param[in] grid A 2d or 3d grid. TransportModelTracerTof(const UnstructuredGrid& grid); - /// Solve for time-of-flight at next timestep. + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[out] tof Array of time-of-flight values. void solveTof(const double* darcyflux, const double* porevolume, diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index e9e812ef..4bc2f1ef 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -446,10 +446,12 @@ namespace Opm - /// Solve for time-of-flight at next timestep. + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[in] degree Polynomial degree of DG basis functions used. /// \param[out] tof_coeff Array of time-of-flight solution coefficients. /// The values are ordered by cell, meaning that @@ -633,7 +635,20 @@ namespace Opm MAT_SIZE_T info = 0; dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info); if (info != 0) { - THROW("Lapack error: " << info); + // Print the local matrix and rhs. + std::cerr << "Failed solving single-cell system Ax = b in cell " << cell + << " with A = \n"; + for (int row = 0; row < n; ++row) { + for (int col = 0; col < n; ++col) { + std::cerr << " " << jac_[row + n*col]; + } + std::cerr << '\n'; + } + std::cerr << "and b = \n"; + for (int row = 0; row < n; ++row) { + std::cerr << " " << rhs_[row] << '\n'; + } + THROW("Lapack error: " << info << " encountered in cell " << cell); } // The solution ends up in rhs_, so we must copy it. std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell); diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp index b9163f8e..0285e1f0 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -31,7 +31,15 @@ namespace Opm class IncompPropertiesInterface; - /// Implements a reordering transport solver for incompressible two-phase flow. + /// Implements a discontinuous Galerkin solver for + /// (single-phase) time-of-flight using reordering. + /// The equation solved is: + /// v \cdot \grad\tau = \phi + /// where v is the fluid velocity, \tau is time-of-flight and + /// \phi is the porosity. This is a boundary value problem, where + /// \tau is specified to be zero on all inflow boundaries. + /// The user may specify the polynomial degree of the basis function space + /// used, but only degrees 0 and 1 are supported so far. class TransportModelTracerTofDiscGal : public TransportModelInterface { public: @@ -39,10 +47,13 @@ namespace Opm /// \param[in] grid A 2d or 3d grid. TransportModelTracerTofDiscGal(const UnstructuredGrid& grid); - /// Solve for time-of-flight at next timestep. + + /// Solve for time-of-flight. /// \param[in] darcyflux Array of signed face fluxes. /// \param[in] porevolume Array of pore volumes. - /// \param[in] source Transport source term. + /// \param[in] source Source term. Sign convention is: + /// (+) inflow flux, + /// (-) outflow flux. /// \param[in] degree Polynomial degree of DG basis functions used. /// \param[out] tof_coeff Array of time-of-flight solution coefficients. /// The values are ordered by cell, meaning that From 22cb94415c7fec78b56d3063c40e1f0f3183c721 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 10:02:47 +0200 Subject: [PATCH 14/15] Fix comment. --- examples/compute_tof.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp index ea55b172..046dae38 100644 --- a/examples/compute_tof.cpp +++ b/examples/compute_tof.cpp @@ -75,7 +75,7 @@ main(int argc, char** argv) { using namespace Opm; - std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; + std::cout << "\n================ Test program for incompressible tof computations ===============\n\n"; parameter::ParameterGroup param(argc, argv, false); std::cout << "--------------- Reading parameters ---------------" << std::endl; From ae04319fa36bfd08f4c60c2958e8282e033f1c16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 9 Oct 2012 12:21:17 +0200 Subject: [PATCH 15/15] Whitespace cleanup. --- .../reorder/TransportModelTracerTofDiscGal.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp index 4bc2f1ef..2d58ec91 100644 --- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -109,9 +109,9 @@ namespace Opm static void cross(const double* a, const double* b, double* res) { - res[0] = a[1]*b[2] - a[2]*b[1]; - res[1] = a[2]*b[0] - a[0]*b[2]; - res[2] = a[0]*b[1] - a[1]*b[0]; + res[0] = a[1]*b[2] - a[2]*b[1]; + res[1] = a[2]*b[0] - a[0]*b[2]; + res[2] = a[0]*b[1] - a[1]*b[0]; } @@ -135,12 +135,12 @@ namespace Opm /// three three-dimensional arrays. static double determinantOf(const double* a0, const double* a1, - const double* a2) + const double* a2) { - return - a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - - a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + - a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); + return + a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) - + a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) + + a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]); } @@ -156,7 +156,7 @@ namespace Opm double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] }; double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] }; double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] }; - return std::fabs(determinantOf(a, b, c) / 6.0); + return std::fabs(determinantOf(a, b, c) / 6.0); }