From bae5e96b7ac5f2dab3ed3d53aeac33792d29af84 Mon Sep 17 00:00:00 2001 From: Max Luebke Date: Fri, 21 Jul 2023 12:36:06 +0200 Subject: [PATCH] refactor: remove support for PhreeqcRM support --- README.md | 4 +- app/CMakeLists.txt | 9 +- app/poet_prm.cpp | 285 ------------------------ include/poet/ChemistryModule.hpp | 32 +-- src/ChemistryModule/ChemistryModule.cpp | 47 ---- 5 files changed, 7 insertions(+), 370 deletions(-) delete mode 100644 app/poet_prm.cpp diff --git a/README.md b/README.md index 06a08679b..4e29aea28 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # POET @@ -69,8 +69,6 @@ following available options: DHT usage. Defaults to _OFF_. - **POET_ENABLE_TESTING**=_boolean_ - enables small set of unit tests (more to come). Defaults to _OFF_. -- **POET_USE_PRM_BACKEND**=_bollean_ - use the PhreeqcRM parallelization instead - of POET's one. Intended for debugging purposes for modellers. ### Example: Build from scratch diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 96c88e150..1aa791dd4 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -1,14 +1,7 @@ configure_file(poet.h.in poet.h) -if(POET_USE_PRM_BACKEND) - set(poet_SRC poet_prm.cpp) -else() - set(poet_SRC poet.cpp) -endif() - -add_executable(poet ${poet_SRC}) +add_executable(poet poet.cpp) target_include_directories(poet PUBLIC "${CMAKE_CURRENT_BINARY_DIR}") target_link_libraries(poet PUBLIC poet_lib MPI::MPI_CXX) -#target_compile_definitions(poet PRIVATE OMPI_SKIP_MPICXX) install(TARGETS poet DESTINATION bin) diff --git a/app/poet_prm.cpp b/app/poet_prm.cpp deleted file mode 100644 index 73734dd3d..000000000 --- a/app/poet_prm.cpp +++ /dev/null @@ -1,285 +0,0 @@ -/* -** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of -** Potsdam) -** -** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam) -** -** POET 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 2 of the License, or (at your option) any later -** version. -** -** POET 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 -** this program; if not, write to the Free Software Foundation, Inc., 51 -** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. -*/ - -#include "poet/ChemistryModule.hpp" -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include - -using namespace std; -using namespace poet; -using namespace Rcpp; - -void set_chem_parameters(poet::ChemistryModule &chem, uint32_t wp_size, - const std::string &database_path) { - chem.SetErrorHandlerMode(1); - chem.SetComponentH2O(false); - chem.SetRebalanceFraction(0.5); - chem.SetRebalanceByCell(true); - chem.UseSolutionDensityVolume(false); - chem.SetPartitionUZSolids(false); - - // Set concentration units - // 1, mg/L; 2, mol/L; 3, kg/kgs - chem.SetUnitsSolution(2); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsPPassemblage(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsExchange(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsSurface(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsGasPhase(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsSSassemblage(1); - // 0, mol/L cell; 1, mol/L water; 2 mol/L rock - chem.SetUnitsKinetics(1); - - // Set representative volume - std::vector rv; - rv.resize(wp_size, 1.0); - chem.SetRepresentativeVolume(rv); - - // Set initial porosity - std::vector por; - por.resize(wp_size, 1); - chem.SetPorosity(por); - - // Set initial saturation - std::vector sat; - sat.resize(wp_size, 1.0); - chem.SetSaturation(sat); - - // Load database - chem.LoadDatabase(database_path); -} - -inline double RunMasterLoop(SimParams ¶ms, RInside &R, Grid &grid, - ChemistryParams &chem_params, - const GridParams &g_params, uint32_t nxyz_master) { - - DiffusionModule diffusion(poet::DiffusionParams(R), grid); - /* Iteration Count is dynamic, retrieving value from R (is only needed by - * master for the following loop) */ - uint32_t maxiter = R.parseEval("mysetup$iterations"); - - double sim_time = .0; - - ChemistryModule chem(grid.GetTotalCellCount(), MPI_COMM_WORLD); - set_chem_parameters(chem, nxyz_master, chem_params.database_path); - chem.RunInitFile(chem_params.input_script); - - chem.SetTimeStep(0); - chem.RunCells(); - - StateMemory *chem_state = grid.RegisterState("state_C", chem.GetPropNames()); - auto &chem_field = chem_state->mem; - chem_field = chem.GetField(); - /* SIMULATION LOOP */ - - double dStartTime = MPI_Wtime(); - for (uint32_t iter = 1; iter < maxiter + 1; iter++) { - uint32_t tick = 0; - // cout << "CPP: Evaluating next time step" << endl; - // R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)"); - - double dt = Rcpp::as( - R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]")); - cout << "CPP: Next time step is " << dt << "[s]" << endl; - - /* displaying iteration number, with C++ and R iterator */ - cout << "CPP: Going through iteration " << iter << endl; - cout << "CPP: R's $iter: " << ((uint32_t)(R.parseEval("mysetup$iter"))) - << ". Iteration" << endl; - - /* run transport */ - // TODO: transport to diffusion - diffusion.simulate(dt); - - grid.PreModuleFieldCopy(tick++); - - cout << "CPP: Chemistry" << endl; - - chem.SetTimeStep(dt); - - chem.SetConcentrations(chem_field); - chem.SetTimeStep(dt); - chem.RunCells(); - chem_field = chem.GetField(); - - grid.WriteFieldsToR(R); - grid.PreModuleFieldCopy(tick++); - - R["req_dt"] = dt; - R["simtime"] = (sim_time += dt); - - R.parseEval("mysetup$req_dt <- req_dt"); - R.parseEval("mysetup$simtime <- simtime"); - - // MDL master_iteration_end just writes on disk state_T and - // state_C after every iteration if the cmdline option - // --ignore-results is not given (and thus the R variable - // store_result is TRUE) - R.parseEvalQ("mysetup <- master_iteration_end(setup=mysetup)"); - - cout << endl - << "CPP: End of *coupling* iteration " << iter << "/" << maxiter - << endl - << endl; - - // MPI_Barrier(MPI_COMM_WORLD); - - } // END SIMULATION LOOP - - R.parseEvalQ("profiling <- list()"); - - R["simtime_chemistry"] = chem.GetChemistryTime(); - R.parseEvalQ("profiling$simtime_chemistry <- simtime_chemistry"); - - chem.MpiWorkerBreak(); - diffusion.end(); - - return MPI_Wtime() - dStartTime; -} - -int main(int argc, char *argv[]) { - double dSimTime, sim_end; - - int world_size, world_rank; - - MPI_Init(&argc, &argv); - - MPI_Comm_size(MPI_COMM_WORLD, &world_size); - - MPI_Comm_rank(MPI_COMM_WORLD, &world_rank); - - if (world_rank == 0) { - cout << "Running POET in version " << poet_version << endl << endl; - } - - if (world_rank > 0) { - { - uint32_t nxyz; - MPI_Bcast(&nxyz, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD); - ChemistryModule worker(nxyz, MPI_COMM_WORLD); - worker.MpiWorker(); - } - - MPI_Barrier(MPI_COMM_WORLD); - cout << "CPP: finished, cleanup of process " << world_rank << endl; - MPI_Finalize(); - - return EXIT_SUCCESS; - } - - /* initialize R runtime */ - RInside R(argc, argv); - - /*Loading Dependencies*/ - // TODO: kann raus - std::string r_load_dependencies = "source('../R_lib/kin_r_library.R');"; - R.parseEvalQ(r_load_dependencies); - - SimParams params(world_rank, world_size); - int pret = params.parseFromCmdl(argv, R); - - if (pret == poet::PARSER_ERROR) { - MPI_Finalize(); - return EXIT_FAILURE; - } else if (pret == poet::PARSER_HELP) { - MPI_Finalize(); - return EXIT_SUCCESS; - } - - cout << "CPP: R Init (RInside) on process " << world_rank << endl; - - R.parseEvalQ("mysetup <- setup"); - // if (world_rank == 0) { // get timestep vector from - // grid_init function ... // - std::string master_init_code = "mysetup <- master_init(setup=setup)"; - R.parseEval(master_init_code); - - Grid grid; - GridParams g_params(R); - - grid.InitModuleFromParams(g_params); - grid.PushbackModuleFlow(poet::DIFFUSION_MODULE_NAME, CHEMISTRY_MODULE_NAME); - grid.PushbackModuleFlow(CHEMISTRY_MODULE_NAME, poet::DIFFUSION_MODULE_NAME); - - params.initVectorParams(R, grid.GetSpeciesCount()); - - // MDL: store all parameters - if (world_rank == 0) { - cout << "CPP: Calling R Function to store calling parameters" << endl; - R.parseEvalQ("StoreSetup(setup=mysetup)"); - } - - if (world_rank == 0) { - cout << "CPP: Init done on process with rank " << world_rank << endl; - } - - // MPI_Barrier(MPI_COMM_WORLD); - - poet::ChemistryParams chem_params(R); - - /* THIS IS EXECUTED BY THE MASTER */ - uint32_t nxyz = grid.GetTotalCellCount(); - MPI_Bcast(&nxyz, 1, MPI_UINT32_T, 0, MPI_COMM_WORLD); - uint32_t nxyz_master = grid.GetTotalCellCount(); - - dSimTime = RunMasterLoop(params, R, grid, chem_params, g_params, nxyz_master); - - cout << "CPP: finished simulation loop" << endl; - - cout << "CPP: start timing profiling" << endl; - - R["simtime"] = dSimTime; - R.parseEvalQ("profiling$simtime <- simtime"); - - string r_vis_code; - r_vis_code = "saveRDS(profiling, file=paste0(fileout,'/timings.rds'));"; - R.parseEval(r_vis_code); - - cout << "CPP: Done! Results are stored as R objects into <" - << params.getOutDir() << "/timings.rds>" << endl; - - MPI_Barrier(MPI_COMM_WORLD); - - cout << "CPP: finished, cleanup of process " << world_rank << endl; - MPI_Finalize(); - - if (world_rank == 0) { - cout << "CPP: done, bye!" << endl; - } - - exit(0); -} diff --git a/include/poet/ChemistryModule.hpp b/include/poet/ChemistryModule.hpp index 917300a43..2db3cd17e 100644 --- a/include/poet/ChemistryModule.hpp +++ b/include/poet/ChemistryModule.hpp @@ -1,4 +1,4 @@ -// Time-stamp: "Last modified 2023-07-12 12:50:53 mluebke" +// Time-stamp: "Last modified 2023-07-21 12:35:23 mluebke" #ifndef CHEMISTRYMODULE_H_ #define CHEMISTRYMODULE_H_ @@ -24,20 +24,6 @@ namespace poet { */ class ChemistryModule : public PhreeqcRM { public: -#ifdef POET_USE_PRM - /** - * Creates a new instance of Chemistry module with given grid cell count and - * using MPI communicator. - * - * Constructor is just a wrapper around PhreeqcRM's constructor, so just calls - * the base class constructor. - * - * \param nxyz Count of grid cells. - * \param communicator MPI communicator where work will be distributed. - */ - ChemistryModule(uint32_t nxyz, MPI_Comm communicator) - : PhreeqcRM(nxyz, communicator), n_cells(nxyz){}; -#else /** * Creates a new instance of Chemistry module with given grid cell count, work * package size and communicator. @@ -59,7 +45,6 @@ public: * Deconstructor, which frees DHT data structure if used. */ ~ChemistryModule(); -#endif /** * Parses the input script and extract information needed during runtime. @@ -92,8 +77,6 @@ public: */ auto GetChemistryTime() const { return this->chem_t; } -#ifndef POET_USE_PRM - /** * Create a new worker instance inside given MPI communicator. * @@ -126,9 +109,9 @@ public: * Enumerating DHT file options */ enum { - DHT_FILES_DISABLED = 0, //!< disabled file output - DHT_FILES_SIMEND, //!< only output of snapshot after simulation - DHT_FILES_ITEREND //!< output snapshots after each iteration + DHT_FILES_DISABLED = 0, //!< disabled file output + DHT_FILES_SIMEND, //!< only output of snapshot after simulation + DHT_FILES_ITEREND //!< output snapshots after each iteration }; /** @@ -281,11 +264,8 @@ public: * \param enabled True if print progressbar, false if not. */ void setProgressBarPrintout(bool enabled); -#endif + protected: -#ifdef POET_USE_PRM - void PrmToPoetField(std::vector &field); -#else enum { CHEM_INIT, CHEM_WP_SIZE, @@ -401,8 +381,6 @@ protected: bool print_progessbar{false}; -#endif - double chem_t = 0.; uint32_t n_cells = 0; diff --git a/src/ChemistryModule/ChemistryModule.cpp b/src/ChemistryModule/ChemistryModule.cpp index ad58881ef..a42d8151c 100644 --- a/src/ChemistryModule/ChemistryModule.cpp +++ b/src/ChemistryModule/ChemistryModule.cpp @@ -12,7 +12,6 @@ #include #include -#ifndef POET_USE_PRM poet::ChemistryModule::ChemistryModule(uint32_t nxyz, uint32_t wp_size, MPI_Comm communicator) : PhreeqcRM(nxyz, 1), group_comm(communicator), wp_size(wp_size) { @@ -42,10 +41,7 @@ poet::ChemistryModule::createWorker(MPI_Comm communicator) { return ChemistryModule(wp_size, wp_size, communicator); } -#endif - void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { -#ifndef POET_USE_PRM if (this->is_master) { int f_type = CHEM_INIT; PropagateFunctionType(f_type); @@ -54,7 +50,6 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { ChemBCast(&count, 1, MPI_INT); ChemBCast(const_cast(input_script_path.data()), count, MPI_CHAR); } -#endif this->RunFile(true, true, false, input_script_path); this->RunString(true, false, false, "DELETE; -all; PRINT; -warnings 0;"); @@ -69,7 +64,6 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1); char surface = (speciesPerModule[4] == 0 ? -1 : 1); -#ifdef POET_USE_PRM std::vector ic1; ic1.resize(this->nxyz * 7, -1); // TODO: hardcoded reaction modules @@ -84,25 +78,8 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) { } this->InitialPhreeqc2Module(ic1); -#else - std::vector ic1; - ic1.resize(this->nxyz * 7, -1); - // TODO: hardcoded reaction modules - for (int i = 0; i < nxyz; i++) { - ic1[i] = 1; // Solution 1 - ic1[nxyz + i] = equilibrium; // Equilibrium 1 - ic1[2 * nxyz + i] = exchange; // Exchange none - ic1[3 * nxyz + i] = surface; // Surface none - ic1[4 * nxyz + i] = -1; // Gas phase none - ic1[5 * nxyz + i] = -1; // Solid solutions none - ic1[6 * nxyz + i] = kinetics; // Kinetics 1 - } - - this->InitialPhreeqc2Module(ic1); -#endif } -#ifndef POET_USE_PRM void poet::ChemistryModule::initializeField(const Field &trans_field) { if (is_master) { @@ -368,27 +345,3 @@ void poet::ChemistryModule::unshuffleField(const std::vector &in_buffer, } } } - -#else // POET_USE_PRM - -inline void poet::ChemistryModule::PrmToPoetField(std::vector &field) { - this->getDumpedField(field); -} - -void poet::ChemistryModule::RunCells() { - PhreeqcRM::RunCells(); - - std::vector tmp_field; - - PrmToPoetField(tmp_field); - this->field = tmp_field; - - for (uint32_t i = 0; i < field.size(); i++) { - uint32_t back_i = static_cast(i / this->nxyz); - uint32_t mod_i = i % this->nxyz; - - field[i] = tmp_field[back_i + (mod_i * this->prop_count)]; - } -} - -#endif