Import the remaining code from opm-core

examples/compute_eikonal_from_files.cpp

@@ -0,0 +1,131 @@
+/*
+  Copyright 2014 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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/flowdiagnostics/AnisotropicEikonal.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <boost/filesystem.hpp>
+#include <memory>
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <numeric>
+#include <iterator>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Warning: unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "-------------------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    ParameterGroup param(argc, argv);
+
+    // Read grid.
+    GridManager grid_manager(param.get<std::string>("grid_filename"));
+    const UnstructuredGrid& grid = *grid_manager.c_grid();
+
+    // Read metric tensor.
+    std::vector<double> metric;
+    {
+        std::ifstream metric_stream(param.get<std::string>("metric_filename").c_str());
+        std::istream_iterator<double> beg(metric_stream);
+        std::istream_iterator<double> end;
+        metric.assign(beg, end);
+        if (int(metric.size()) != grid.number_of_cells*grid.dimensions*grid.dimensions) {
+            OPM_THROW(std::runtime_error, "Size of metric field differs from (dim^2 * number of cells).");
+        }
+    }
+
+    // Read starting cells.
+    std::vector<int> startcells;
+    {
+        std::ifstream start_stream(param.get<std::string>("startcells_filename").c_str());
+        std::istream_iterator<int> beg(start_stream);
+        std::istream_iterator<int> end;
+        startcells.assign(beg, end);
+    }
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    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 (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        param.writeParam(output_dir + "/eikonal.param");
+    }
+
+    // Issue a warning if any parameters were unused.
+    warnIfUnusedParams(param);
+
+    // Solve eikonal equation.
+    Opm::time::StopWatch timer;
+    timer.start();
+    std::vector<double> solution;
+    AnisotropicEikonal2d ae(grid);
+    ae.solve(metric.data(), startcells, solution);
+    timer.stop();
+    double tt = timer.secsSinceStart();
+    std::cout << "Eikonal solver took: " << tt << " seconds." << std::endl;
+
+    // Output.
+    if (output) {
+        std::string filename = output_dir + "/solution.txt";
+        std::ofstream stream(filename.c_str());
+        stream.precision(16);
+        std::copy(solution.begin(), solution.end(), std::ostream_iterator<double>(stream, "\n"));
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

examples/compute_initial_state.cpp

@@ -0,0 +1,174 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2017 IRIS
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/simulator/initStateEquil.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/props/phaseUsageFromDeck.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/utility/compressedToCartesian.hpp>
+
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
+
+#include <boost/filesystem.hpp>
+
+#include <fstream>
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+
+    void outputData(const std::string& output_dir,
+                    const std::string& name,
+                    const std::vector<double>& data)
+    {
+        std::ostringstream fname;
+        fname << output_dir << "/" << name;
+        boost::filesystem::path fpath = fname.str();
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        fname << "/" << "initial.txt";
+        std::ofstream file(fname.str().c_str());
+        if (!file) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
+        }
+        std::copy(data.begin(), data.end(), std::ostream_iterator<double>(file, "\n"));
+    }
+
+    /// Convert saturations from a vector of individual phase saturation vectors
+    /// to an interleaved format where all values for a given cell come before all
+    /// values for the next cell, all in a single vector.
+    template <class FluidSystem>
+    void convertSats(std::vector<double>& sat_interleaved, const std::vector< std::vector<double> >& sat, const Opm::PhaseUsage& pu)
+    {
+        assert(sat.size() == 3);
+        const auto nc = sat[0].size();
+        const auto np = sat_interleaved.size() / nc;
+        for (size_t c = 0; c < nc; ++c) {
+            if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                const int opos = pu.phase_pos[Opm::BlackoilPhases::Liquid];
+                const std::vector<double>& sat_p = sat[ FluidSystem::oilPhaseIdx];
+                sat_interleaved[np*c + opos] = sat_p[c];
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                const int wpos = pu.phase_pos[Opm::BlackoilPhases::Aqua];
+                const std::vector<double>& sat_p = sat[ FluidSystem::waterPhaseIdx];
+                sat_interleaved[np*c + wpos] = sat_p[c];
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                const int gpos = pu.phase_pos[Opm::BlackoilPhases::Vapour];
+                const std::vector<double>& sat_p = sat[ FluidSystem::gasPhaseIdx];
+                sat_interleaved[np*c + gpos] = sat_p[c];
+            }
+        }
+    }
+
+} // anon namespace
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    // Setup.
+    ParameterGroup param(argc, argv);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+    const std::string deck_filename = param.get<std::string>("deck_filename");
+    Opm::ParseContext parseContext;
+    Opm::Parser parser;
+    const Opm::Deck& deck = parser.parseFile(deck_filename , parseContext);
+    const Opm::EclipseState eclipseState(deck, parseContext);
+    const double grav = param.getDefault("gravity", unit::gravity);
+    GridManager gm(eclipseState.getInputGrid());
+    const UnstructuredGrid& grid = *gm.c_grid();
+    warnIfUnusedParams(param);
+
+    // Create material law manager.
+    std::vector<int> compressedToCartesianIdx
+        = Opm::compressedToCartesian(grid.number_of_cells, grid.global_cell);
+
+    typedef FluidSystems::BlackOil<double> FluidSystem;
+
+    // Forward declaring the MaterialLawManager template.
+    typedef Opm::ThreePhaseMaterialTraits<double,
+    /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
+    /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
+    /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> MaterialTraits;
+    typedef Opm::EclMaterialLawManager<MaterialTraits> MaterialLawManager;
+
+    MaterialLawManager materialLawManager = MaterialLawManager();
+    materialLawManager.initFromDeck(deck, eclipseState, compressedToCartesianIdx);
+
+    // Initialisation.
+    //initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
+    BlackoilState state( UgGridHelpers::numCells(grid) , UgGridHelpers::numFaces(grid), 3);
+    FluidSystem::initFromDeck(deck, eclipseState);
+    PhaseUsage pu = phaseUsageFromDeck(deck);
+
+    typedef EQUIL::DeckDependent::InitialStateComputer<FluidSystem> ISC;
+
+    ISC isc(materialLawManager, eclipseState, grid, grav);
+
+    const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
+    const int oilpos = FluidSystem::oilPhaseIdx;
+    const int waterpos = FluidSystem::waterPhaseIdx;
+    const int ref_phase = oil ? oilpos : waterpos;
+
+    state.pressure() = isc.press()[ref_phase];
+    convertSats<FluidSystem>(state.saturation(), isc.saturation(), pu);
+    state.gasoilratio() = isc.rs();
+    state.rv() = isc.rv();
+
+    // Output.
+    const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
+    outputData(output_dir, "pressure", state.pressure());
+    outputData(output_dir, "saturation", state.saturation());
+    outputData(output_dir, "rs", state.gasoilratio());
+    outputData(output_dir, "rv", state.rv());
+}
+catch (const std::exception& e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

examples/compute_tof_from_files.cpp

@@ -0,0 +1,213 @@
+/*
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/utility/SparseTable.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/props/IncompPropertiesBasic.hpp>
+#include <opm/core/props/IncompPropertiesFromDeck.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/flowdiagnostics/TofReorder.hpp>
+#include <opm/core/flowdiagnostics/TofDiscGalReorder.hpp>
+
+#include <memory>
+#include <boost/filesystem.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include <iterator>
+#include <fstream>
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Warning: unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "-------------------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    ParameterGroup param(argc, argv);
+
+    // Read grid.
+    GridManager grid_manager(param.get<std::string>("grid_filename"));
+    const UnstructuredGrid& grid = *grid_manager.c_grid();
+
+    // Read porosity, compute pore volume.
+    std::vector<double> porevol;
+    {
+        std::ifstream poro_stream(param.get<std::string>("poro_filename").c_str());
+        std::istream_iterator<double> beg(poro_stream);
+        std::istream_iterator<double> end;
+        porevol.assign(beg, end); // Now contains poro.
+        if (int(porevol.size()) != grid.number_of_cells) {
+            OPM_THROW(std::runtime_error, "Size of porosity field differs from number of cells.");
+        }
+        for (int i = 0; i < grid.number_of_cells; ++i) {
+            porevol[i] *= grid.cell_volumes[i];
+        }
+    }
+
+    // Read flux.
+    std::vector<double> flux;
+    {
+        std::ifstream flux_stream(param.get<std::string>("flux_filename").c_str());
+        std::istream_iterator<double> beg(flux_stream);
+        std::istream_iterator<double> end;
+        flux.assign(beg, end);
+        if (int(flux.size()) != grid.number_of_faces) {
+            OPM_THROW(std::runtime_error, "Size of flux field differs from number of faces.");
+        }
+    }
+
+    // Read source terms.
+    std::vector<double> src;
+    {
+        std::ifstream src_stream(param.get<std::string>("src_filename").c_str());
+        std::istream_iterator<double> beg(src_stream);
+        std::istream_iterator<double> end;
+        src.assign(beg, end);
+        if (int(src.size()) != grid.number_of_cells) {
+            OPM_THROW(std::runtime_error, "Size of source term field differs from number of cells.");
+        }
+    }
+
+    const bool compute_tracer = param.getDefault("compute_tracer", false);
+    Opm::SparseTable<int> tracerheads;
+    if (compute_tracer) {
+        std::ifstream tr_stream(param.get<std::string>("tracerheads_filename").c_str());
+        int num_rows;
+        tr_stream >> num_rows;
+        for (int row = 0; row < num_rows; ++row) {
+            int row_size;
+            tr_stream >> row_size;
+            std::vector<int> rowdata(row_size);
+            for (int elem = 0; elem < row_size; ++elem) {
+                tr_stream >> rowdata[elem];
+            }
+            tracerheads.appendRow(rowdata.begin(), rowdata.end());
+        }
+    }
+
+    // Choice of tof solver.
+    bool use_dg = param.getDefault("use_dg", false);
+    bool use_multidim_upwind = false;
+    // Need to initialize dg solver here, since it uses parameters now.
+    std::unique_ptr<Opm::TofDiscGalReorder> dg_solver;
+    if (use_dg) {
+        dg_solver.reset(new Opm::TofDiscGalReorder(grid, param));
+    } else {
+        use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
+    }
+
+    // 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 (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    // Issue a warning if any parameters were unused.
+    warnIfUnusedParams(param);
+
+    // Solve time-of-flight.
+    Opm::time::StopWatch transport_timer;
+    transport_timer.start();
+    std::vector<double> tof;
+    std::vector<double> tracer;
+    if (use_dg) {
+        if (compute_tracer) {
+            dg_solver->solveTofTracer(&flux[0], &porevol[0], &src[0], tracerheads, tof, tracer);
+        } else {
+            dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof);
+        }
+    } else {
+        Opm::TofReorder tofsolver(grid, use_multidim_upwind);
+        if (compute_tracer) {
+            tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tracerheads, tof, tracer);
+        } else {
+            tofsolver.solveTof(&flux[0], &porevol[0], &src[0], tof);
+        }
+    }
+    transport_timer.stop();
+    double tt = transport_timer.secsSinceStart();
+    std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+
+    // Output.
+    if (output) {
+        std::string tof_filename = output_dir + "/tof.txt";
+        std::ofstream tof_stream(tof_filename.c_str());
+        tof_stream.precision(16);
+        std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
+        if (compute_tracer) {
+            std::string tracer_filename = output_dir + "/tracer.txt";
+            std::ofstream tracer_stream(tracer_filename.c_str());
+            tracer_stream.precision(16);
+            const int nt = tracer.size()/grid.number_of_cells;
+            for (int i = 0; i < nt*grid.number_of_cells; ++i) {
+                tracer_stream << tracer[i] << (((i + 1) % nt == 0) ? '\n' : ' ');
+            }
+        }
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

examples/diagnose_relperm.cpp

@@ -0,0 +1,105 @@
+/*
+  Copyright 2015 Statoil ASA.
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/props/satfunc/RelpermDiagnostics.hpp>
+#include <opm/core/props/satfunc/RelpermDiagnostics_impl.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/OpmLog/EclipsePRTLog.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/algorithm/string/case_conv.hpp>
+#include <memory>
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include <iterator>
+
+
+void usage() {
+    std::cout << std::endl << 
+        "Usage: diagnose_relperm <eclipseFile>" << std::endl;
+}
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+    if (argc <= 1) {
+        usage();
+        exit(1);
+    } 
+    const char* eclipseFilename = argv[1];
+    Parser parser;
+    Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }, 
+                              { ParseContext::PARSE_UNKNOWN_KEYWORD, InputError::IGNORE},
+                              { ParseContext::PARSE_RANDOM_TEXT, InputError::IGNORE},
+                              { ParseContext::UNSUPPORTED_SCHEDULE_GEO_MODIFIER, InputError::IGNORE},
+                              { ParseContext::UNSUPPORTED_COMPORD_TYPE, InputError::IGNORE},
+                              { ParseContext::UNSUPPORTED_INITIAL_THPRES, InputError::IGNORE},
+                              { ParseContext::INTERNAL_ERROR_UNINITIALIZED_THPRES, InputError::IGNORE}
+                             });
+    Opm::Deck deck = parser.parseFile(eclipseFilename, parseContext);
+    Opm::EclipseState eclState( deck, parseContext );
+
+    GridManager gm(eclState.getInputGrid());
+    const UnstructuredGrid& grid = *gm.c_grid();
+    using boost::filesystem::path; 
+    path fpath(eclipseFilename);
+    std::string baseName;
+    if (boost::to_upper_copy(path(fpath.extension()).string())== ".DATA") {
+        baseName = path(fpath.stem()).string();
+    } else {
+        baseName = path(fpath.filename()).string();
+    }
+
+    std::string logFile = baseName + ".SATFUNCLOG";
+    std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile, Log::DefaultMessageTypes);
+    OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
+    prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true, false));
+    std::shared_ptr<StreamLog> streamLog = std::make_shared<EclipsePRTLog>(std::cout, Log::DefaultMessageTypes);
+    OpmLog::addBackend( "STREAMLOG" , streamLog );
+    streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10));
+    streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true, true));
+    RelpermDiagnostics diagnostic;
+    diagnostic.diagnosis(eclState, deck, grid);
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}