diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index 7ad97fad9..bcf136218 100644 --- a/src/executables/CMakeLists.txt +++ b/src/executables/CMakeLists.txt @@ -42,6 +42,7 @@ include_directories(${ATS_SOURCE_DIR}/src/pks/energy) include_directories(${ATS_SOURCE_DIR}/src/pks/flow) include_directories(${ATS_SOURCE_DIR}/src/pks/deform) include_directories(${ATS_SOURCE_DIR}/src/pks/transport) +include_directories(${ATS_SOURCE_DIR}/src/pks/ecosim) include_directories(${ATS_SOURCE_DIR}/src/operators/upwinding) include_directories(${ATS_SOURCE_DIR}/src/operators/advection) include_directories(${ATS_SOURCE_DIR}/src/operators/deformation) @@ -63,6 +64,12 @@ include_evaluators_directories(LISTNAME ATS_BGC_REG_INCLUDES) include_evaluators_directories(LISTNAME ATS_MPC_REG_INCLUDES) include_evaluators_directories(LISTNAME SED_TRANSPORT_REG_INCLUDES) + +include_evaluators_directories(LISTNAME ATS_ECOSIM_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_ECOSIM_RELATIONS_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_ECOSIM_DATA_REG_INCLUDES) + + set(ats_src_files ats_mesh_factory.cc coordinator.cc @@ -122,6 +129,19 @@ set(ats_link_libs ats_transport_relations ) +if (ENABLE_ECOSIM) + list(APPEND ats_link_libs + ats_ecosim + ats_ecosim_data + # ats_ecosim_relations + ) +endif() + +#In theory covered by ECOSIM_LIBRARIES: +#If not place this in an if statement? +# ats_ecosim +# ats_ecosim_data +# ats_ecosim_relations # note, we can be inclusive here, because if they aren't enabled, # these won't be defined and will result in empty strings. @@ -137,7 +157,8 @@ set(tpl_link_libs ${HYPRE_LIBRARIES} ${HDF5_LIBRARIES} ${CLM_LIBRARIES} - ) + ${ECOSIM_LIBRARIES} + ) add_amanzi_library(ats_executable SOURCE ${ats_src_files} diff --git a/src/executables/ats_registration_files.hh b/src/executables/ats_registration_files.hh index 9dcae1462..6d348689f 100644 --- a/src/executables/ats_registration_files.hh +++ b/src/executables/ats_registration_files.hh @@ -25,3 +25,8 @@ #ifdef ALQUIMIA_ENABLED #include "pks_chemistry_reg.hh" #endif +#ifdef ECOSIM_ENABLED +# include "ats_ecosim_registration.hh" +# include "ats_ecosim_relations_registration.hh" +# include "ats_ecosim_data_registration.hh" +#endif diff --git a/src/pks/CMakeLists.txt b/src/pks/CMakeLists.txt index f6c5fc364..99421a746 100644 --- a/src/pks/CMakeLists.txt +++ b/src/pks/CMakeLists.txt @@ -54,3 +54,7 @@ add_subdirectory(deform) add_subdirectory(surface_balance) add_subdirectory(biogeochemistry) add_subdirectory(mpc) + +if (ENABLE_ECOSIM) + add_subdirectory(ecosim) +endif() diff --git a/src/pks/ecosim/BGCEngine.cc b/src/pks/ecosim/BGCEngine.cc new file mode 100644 index 000000000..2900bcca5 --- /dev/null +++ b/src/pks/ecosim/BGCEngine.cc @@ -0,0 +1,119 @@ +/* + Basic architecture based on the Alquima interfece adapted for + use in the ATSEcoSIM PK + + Copyright 2010-202x held jointly by LANS/LANL, LBNL, and PNNL. + Amanzi is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Authors: Jeffrey Johnson + Sergi Molins + + This implements the Alquimia chemistry engine. +*/ + +#include +#include +#include +#include +#include "BGCEngine.hh" +#include "errors.hh" +#include "exceptions.hh" + +// Support for manipulating floating point exception handling. +#ifdef _GNU_SOURCE +#define AMANZI_USE_FENV +#include +#endif + +namespace Amanzi { +namespace EcoSIM { + +BGCEngine::BGCEngine(const std::string& engineName, + const std::string& inputFile) : + bgc_engine_name_(engineName), + bgc_engine_inputfile_(inputFile) +{ + Errors::Message msg; + + CreateBGCInterface(bgc_engine_name_.c_str(), + &bgc_); + +} + +BGCEngine::~BGCEngine() +{ + bgc_.Shutdown(); + + //Did I forget to implement this? + //FreeBGCProperties(&props); + //FreeBGCState(&state); + //FreeBGCAuxiliaryData(&aux_data); + //FreeAlquimiaEngineStatus(&chem_status_); +} + +const BGCSizes& +BGCEngine::Sizes() const +{ + return sizes_; +} + +void BGCEngine::InitState(BGCProperties& properties, + BGCState& state, + BGCAuxiliaryData& aux_data, + int ncells_per_col_, + int num_components, + int num_columns) +{ + AllocateBGCProperties(&sizes_, &properties, ncells_per_col_, num_columns); + AllocateBGCState(&sizes_, &state, ncells_per_col_, num_components, num_columns); +} + +void BGCEngine::FreeState(BGCProperties& properties, + BGCState& state, + BGCAuxiliaryData& aux_data) +{ + FreeBGCProperties(&properties); + FreeBGCState(&state); +} + +void BGCEngine::DataTest() { + + bgc_.DataTest(); +} + +bool BGCEngine::Setup(BGCProperties& properties, + BGCState& state, + BGCSizes& sizes_, + int num_iterations, + int num_columns, + int ncells_per_col_) +{ + bgc_.Setup(&properties, + &state, + &sizes_, + num_iterations, + num_columns, + ncells_per_col_); + +} + +bool BGCEngine::Advance(const double delta_time, + BGCProperties& properties, + BGCState& state, + BGCSizes& sizes_, + int num_iterations, + int num_columns) +{ + bgc_.Advance(delta_time, + &properties, + &state, + &sizes_, + num_iterations, + num_columns); + +} + +} // namespace +} // namespace diff --git a/src/pks/ecosim/BGCEngine.hh b/src/pks/ecosim/BGCEngine.hh new file mode 100644 index 000000000..75d0dc527 --- /dev/null +++ b/src/pks/ecosim/BGCEngine.hh @@ -0,0 +1,107 @@ +/* + ATS-EcoSIM, Code Adapted for use from Alquimia + + Copyright 2010-202x held jointly by LANS/LANL, LBNL, and PNNL. + Amanzi is released under the three-clause BSD License. + The terms of use and "as is" disclaimer for this license are + provided in the top-level COPYRIGHT file. + + Author: Jeffrey Johnson + + This is a point of contact for the chemistry engine exposed by Alquimia + to the rest of Amanzi--it provides the ability to enforce geochemical + conditions and to integrate reactions given a chemical configuration. +*/ + +#ifndef BGC_ENGINE_HH_ +#define BGC_ENGINE_HH_ + +#include +#include +#include + +#include "BGC_memory.hh" +#include "BGC_containers.hh" + +#include "VerboseObject.hh" + +namespace Amanzi { +namespace EcoSIM { + +class BGCEngine { + public: + + // Constructs a chemistry engine using the given engine (backend) name and input file. + BGCEngine(const std::string& engineName, const std::string& inputFile); + + // Destructor. + ~BGCEngine(); + + // Returns the name of the backend that does the chemistry. + const std::string& Name() const; + + // Returns true if the chemistry engine is thread-safe, false if not. + bool IsThreadSafe() const; + + // Returns a reference to a "sizes" object that can be queried to find the sizes of the various + // arrays representing the geochemical state within the engine. + const BGCSizes& Sizes() const; + + // Initializes the data structures that hold the chemical state information. + void InitState(BGCProperties& properties, + BGCState& state, + BGCAuxiliaryData& aux_data, + int ncells_per_col_, + int num_components, + int num_columns); + + // Frees the data structures that hold the chemical state information. + void FreeState(BGCProperties& properties, + BGCState& state, + BGCAuxiliaryData& aux_data); + + void DataTest(); + + bool Setup(BGCProperties& properties, + BGCState& state, + BGCSizes& sizes, + int num_iterations, + int num_columns, + int ncells_per_col_); + + bool Advance(const double delta_time, + BGCProperties& properties, + BGCState& state, + BGCSizes& sizes, + int num_iterations, + int num_columns); + + void CopyBGCState(const BGCState* const source, + BGCState* destination); + void CopyBGCProperties(const BGCProperties* const source, + BGCProperties* destination); + + private: + + // bgc data structures. + bool bgc_initialized_; + void* engine_state_; + BGCSizes sizes_; + BGCInterface bgc_; + + Teuchos::RCP vo_; + // Back-end engine name and input file. + std::string bgc_engine_name_; + std::string bgc_engine_inputfile_; + + // forbidden. + BGCEngine(); + BGCEngine(const BGCEngine&); + BGCEngine& operator=(const BGCEngine&); + +}; + +} // namespace +} // namespace + +#endif diff --git a/src/pks/ecosim/CMakeLists.txt b/src/pks/ecosim/CMakeLists.txt new file mode 100644 index 000000000..d6ccc750c --- /dev/null +++ b/src/pks/ecosim/CMakeLists.txt @@ -0,0 +1,140 @@ +# -*- mode: cmake -*- +#Everything here depends on EcoSIM +#so we put everything in an if statement + +if(ENABLE_ECOSIM) + add_subdirectory(constitutive_relations) + add_subdirectory(data) + + get_property(AMANZI_TPLS_DIR GLOBAL PROPERTY AMANZI_TPLS_DIR) + + set(ECOSIM_INSTALL_PREFIX ${ECOSIM_DIR}/ecosim) + set(ECOSIM_LIB_LOCATION ${ECOSIM_DIR}/ecosim/local/lib) + set(ECOSIM_BUILD_PREFIX ${ECOSIM_DIR}/ecosim/build) + #set(ECOSIM_CMAKE_BINARY_DIR ${ECOSIM_DIR}/ecosim/build/Linux-x86_64-static-not-set-mpicc-Release) + set(NETCDF_LIB ${ECOSIM_DIR}/lib) + + message("In ATS-EcoSIM CMakeLists:") + message("ECOSIM_DIR:" ${ECOSIM_DIR}) + message("ECOSIM_INSTALL_PREFIX:" ${ECOSIM_INSTALL_PREFIX}) + message("ECOSIM_LIB_LOCATION:" ${ECOSIM_LIB_LOCATION}) + message("ECOSIM_BUILD_PREFIX:" ${ECOSIM_BUILD_PREFIX}) + + include_directories(${ATS_SOURCE_DIR}/src/pks) + include_directories(${ATS_SOURCE_DIR}/src/pks/ecosim) + include_directories(${ATS_SOURCE_DIR}/src/pks/ecosim/data) + #include_directories(${ATS_SOURCE_DIR}/src/pks/ecosim/constitutive_relations/bulk_density) + include_directories(${ATS_SOURCE_DIR}/src/pks/ecosim/constitutive_relations/hydraulic_conductivity) + #include_directories(${ATS_SOURCE_DIR}/src/pks/ecosim/constitutive_relations/matric_pressure) + + message("At include_directories") + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Utils/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Minimath/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Modelconfig/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Modelforc/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Mesh/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Modelpars/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Balances/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Ecosim_datatype/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/SoilPhys/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/SurfPhys/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/PhysData/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/SnowPhys/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Ecosim_mods/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Prescribed_pheno/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Plant_bgc/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/APIs/) + include_directories(${ECOSIM_BUILD_PREFIX}/f90src/APIData/) + + include_directories(${ECOSIM_INCLUDE_DIRS}) + + message("EcoSIM inc dirs: ") + message(eco_inc_file="${ECOSIM_INCLUDE_DIRS}") + + # ATS EcoSIM pk + # For adding F90 files add the F90 file to the ecosim source files and + # inc files. The inc file also needs a header + + #set(ats_ecosim_src_files + # EcoSIM_ATS_interface.cc + # BGCEngine.cc + # ecosim_wrappers.F90 + # data/bgc_fortran_memory_mod.F90 + #) + + #testing using link libs: + #file(GLOB ECOLIBS ${ECOSIM_LIB_LOCATION}/*.a) + + set(ats_ecosim_src_files + EcoSIM_ATS_interface.cc + BGCEngine.cc + data/bgc_fortran_memory_mod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/BGC_containers.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSCPLMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSEcoSIMInitMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSEcoSIMAdvanceMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSUtilsMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/SharedDataMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/c_f_interface_module.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ecosim_wrappers.F90 + ) + + set(ats_ecosim_inc_files + EcoSIM_ATS_interface.hh + BGCEngine.hh + ecosim_interface.h + ) + + file(GLOB ECOSIM_LIBRARIES + ${ECOSIM_LIB_LOCATION}/*.a + ) + + + find_package(NetCDF REQUIRED) + + set(ats_ecosim_link_libs + ${Teuchos_LIBRARIES} + ${Epetra_LIBRARIES} + ${ECOSIM_LIBRARIES} + error_handling + atk + mesh + data_structures + whetstone + operators + solvers + time_integration + state + pks + chemistry_pk + ats_pks + ats_eos + ats_operators + ats_ecosim_data + ats_ecosim_relations + gfortran + ) + + message(STATUS "ats_ecosim_link_libs: ${ats_ecosim_link_libs}") + + message(inc_files="${ats_ecosim_inc_files}") + + add_amanzi_library(ats_ecosim + SOURCE ${ats_ecosim_src_files} + HEADERS ${ats_ecosim_inc_files} + LINK_LIBS ${ats_ecosim_link_libs}) + + #================================================ + # register evaluators/factories/pks + + register_evaluator_with_factory( + HEADERFILE EcoSIM_ATS_interface_reg.hh + LISTNAME ATS_ECOSIM_REG + ) + + generate_evaluators_registration_header( + HEADERFILE ats_ecosim_registration.hh + LISTNAME ATS_ECOSIM_REG + INSTALL True + ) +endif() diff --git a/src/pks/ecosim/EcoSIM_ATS_interface.cc b/src/pks/ecosim/EcoSIM_ATS_interface.cc new file mode 100644 index 000000000..5b784e9cd --- /dev/null +++ b/src/pks/ecosim/EcoSIM_ATS_interface.cc @@ -0,0 +1,1378 @@ +/*-------------------------------------------------------------------------- + ATS + + License: see $ATS_DIR/COPYRIGHT + Author: Andrew Graus + + This is the main PK for the EcoSIM-ATS interface. This code is written + following the example of Alquimia with some additional code from the + SimpleBGC code for walking the columns. + + The idea is to take the basic code used by alquimia and repurpose it so + that it works on a column by column basis instead of a cell by cell basis + + --------------------------------------------------------------------------*/ + +#include +#include +#include + +// TPLs +#include "Epetra_MultiVector.h" +#include "Epetra_Vector.h" +#include "Epetra_SerialDenseVector.h" +#include "Epetra_SerialDenseMatrix.h" +#include "Epetra_MpiComm.h" +#include "Epetra_Map.h" +#include "Teuchos_RCPDecl.hpp" +#include "Teuchos_ParameterList.hpp" + +// Amanzi +#include "errors.hh" +#include "exceptions.hh" +#include "Mesh.hh" + +// include custom evaluators here +// #include "hydraulic_conductivity_evaluator.hh" + +#include "PK_Helpers.hh" +#include "EcoSIM_ATS_interface.hh" + +namespace Amanzi { +namespace EcoSIM { + +EcoSIM::EcoSIM(Teuchos::ParameterList& pk_tree, + const Teuchos::RCP& global_list, + const Teuchos::RCP& S, + const Teuchos::RCP& solution): + PK_Physical_Default(pk_tree, global_list, S, solution), + PK(pk_tree, global_list, S, solution), + ncells_per_col_(-1), + saved_time_(0.0) + { + //grab the surface and subsurface domains + domain_ = Keys::readDomain(*plist_, "domain", "domain"); + domain_surface_ = Keys::readDomainHint(*plist_, domain_, "subsurface", "surface"); + + // transport + mole_fraction_key_ = Keys::readKey(*plist_, domain_, "mole fraction", "mole_fraction"); + //mole_fraction components are accessed by mole_fraction[i][c] where i is the component and c is the cell + + //Flow + porosity_key_ = Keys::readKey(*plist_, domain_, "porosity", "porosity"); + saturation_liquid_key_ = Keys::readKey(*plist_, domain_, "saturation liquid", "saturation_liquid"); + saturation_gas_key_ = Keys::readKey(*plist_,domain_,"saturation gas", "saturation_gas"); + saturation_ice_key_ = Keys::readKey(*plist_,domain_,"saturation ice", "saturation_ice"); + water_content_key_ = Keys::readKey(*plist_,domain_,"water content","water_content"); + //relative_permeability_key_ = Keys::readKey(*plist_,domain_,"relative permeability","relative_permeability"); + //matric_pressure_key_ = Keys::readKey(*plist_,domain_,"matric pressure","matric_pressure"); + cap_pres_key_ = Keys::readKey(*plist_, domain_, "capillary pressure key", "capillary_pressure_gas_liq"); + //liquid_density_key_ = Keys::readKey(*plist_, domain_, "mass density liquid", "mass_density_liquid"); + liquid_density_key_ = Keys::readKey(*plist_, domain_, "molar density liquid", "molar_density_liquid"); + ice_density_key_ = Keys::readKey(*plist_, domain_, "mass density ice", "mass_density_ice"); + gas_density_key_ = Keys::readKey(*plist_, domain_,"mass density gas", "mass_density_gas"); + gas_density_key_test_ = Keys::readKey(*plist_, domain_, "mass density gas", "mass_density_gas"); + rock_density_key_ = Keys::readKey(*plist_, domain_, "density rock", "density_rock"); + + //energy + T_key_ = Keys::readKey(*plist_, domain_, "temperature", "temperature"); + thermal_conductivity_key_ = Keys::readKey(*plist_, domain_, "thermal conductivity", "thermal_conductivity"); + + //Sources + surface_water_source_key_ = Keys::readKey(*plist_, domain_surface_, "surface water source", "water_source"); + surface_energy_source_key_ = + Keys::readKey(*plist_, domain_surface_, "surface energy source", "total_energy_source"); + subsurface_water_source_key_ = + Keys::readKey(*plist_, domain_, "subsurface water source", "water_source"); + subsurface_energy_source_key_ = + Keys::readKey(*plist_, domain_, "subsurface energy source", "total_energy_source"); + surface_energy_source_ecosim_key_ = + Keys::readKey(*plist_, domain_surface_, "surface energy source ecosim", "ecosim_source"); + surface_water_source_ecosim_key_ = + Keys::readKey(*plist_, domain_surface_, "surface water source ecosim", "ecosim_water_source"); + + subsurface_energy_source_ecosim_key_ = + Keys::readKey(*plist_, domain_, "subsurface energy source ecosim", "subsurface_ecosim_source"); + subsurface_water_source_ecosim_key_ = + Keys::readKey(*plist_, domain_, "subsurface water source ecosim", "subsurface_ecosim_water_source"); + + //Other + cell_volume_key_ = Keys::readKey(*plist_, domain_, "cell volume", "cell_volume"); + //ecosim_aux_data_key_ = Keys::readKey(*plist_, domain_, "ecosim aux data", "ecosim_aux_data"); + f_wp_key_ = Keys::readKey(*plist_, domain_, "porosity", "porosity"); + f_root_key_ = Keys::readKey(*plist_, domain_, "porosity", "porosity"); + + //Custom Evaluator keys + hydraulic_conductivity_key_ = Keys::readKey(*plist_, domain_, "hydraulic conductivity", "hydraulic_conductivity"); + //bulk_density_key_ = Keys::readKey(*plist_, domain_, "bulk density", "bulk_density"); + + //Surface balance items + sw_key_ = + Keys::readKey(*plist_, domain_surface_, "incoming shortwave radiation", "incoming_shortwave_radiation"); + lw_key_ = + Keys::readKey(*plist_,domain_surface_, "incoming longwave radiation", "incoming_longwave_radiation"); + air_temp_key_ = Keys::readKey(*plist_, domain_surface_, "air temperature", "air_temperature"); + vp_air_key_ = Keys::readKey(*plist_, domain_surface_, "vapor pressure air", "vapor_pressure_air"); + wind_speed_key_ = Keys::readKey(*plist_, domain_surface_, "wind speed", "wind_speed"); + p_rain_key_ = Keys::readKey(*plist_, domain_surface_, "precipitation rain", "precipitation_rain"); + p_snow_key_ = Keys::readKey(*plist_, domain_surface_, "precipitation snow", "precipitation_snow"); + p_total_key_ = Keys::readKey(*plist_, domain_surface_, "precipitation total", "precipitation_total"); + elev_key_ = Keys::readKey(*plist_, domain_surface_, "elevation", "elevation"); + aspect_key_ = Keys::readKey(*plist_, domain_surface_, "aspect", "aspect"); + slope_key_ = Keys::readKey(*plist_, domain_surface_, "slope", "slope_magnitude"); + snow_depth_key_ = Keys::readKey(*plist_, domain_surface_, "snow depth", "snow_depth"); + snow_albedo_key_ = Keys::readKey(*plist_, domain_surface_, "snow_albedo", "snow_albedo"); + snow_temperature_key_ = Keys::readKey(*plist_, domain_surface_, "snow temperature", "snow_temperature"); + + //Canopy hold over vars for EcoSIM + canopy_lw_key_ = Keys::readKey(*plist_, domain_surface_, "canopy longwave radiation", "canopy_longwave_radiation"); + canopy_latent_heat_key_ = Keys::readKey(*plist_, domain_surface_, "canopy latent heat", "canopy_latent_heat"); + canopy_sensible_heat_key_ = Keys::readKey(*plist_, domain_surface_, "canopy sensible heat", "canopy_sensible_heat"); + canopy_surface_water_key_ = Keys::readKey(*plist_, domain_surface_, "canopy surface water", "canopy_surface_water"); + transpiration_key_ = Keys::readKey(*plist_, domain_surface_, "transpiration", "transpiration"); + evaporation_canopy_key_ = Keys::readKey(*plist_, domain_surface_, "evaporation canopy", "evaporation_canopy"); + evaporation_ground_key_ = Keys::readKey(*plist_, domain_surface_, "evaporation ground", "evaporation_ground"); + evaporation_litter_key_ = Keys::readKey(*plist_, domain_surface_, "evaporation litter", "evaporation_litter"); + evaporation_snow_key_ = Keys::readKey(*plist_, domain_surface_, "evaporation snow", "evaporation_snow"); + sublimation_snow_key_ = Keys::readKey(*plist_, domain_surface_, "sublimation snow", "sublimation_snow"); + + //Plant Phenology Datasets + lai_key_ = Keys::readKey(*plist_, domain_surface_, "LAI", "LAI"); + sai_key_ = Keys::readKey(*plist_, domain_surface_, "SAI", "SAI"); + v_type_key_ = Keys::readKey(*plist_, domain_surface_, "vegetation type", "vegetation_type"); + + //Atmospheric abundance keys + /*atm_n2_ = plist_->get("atmospheric N2"); + atm_o2_ = plist_->get("atmospheric O2"); + atm_co2_ = plist_->get("atmospheric CO2"); + atm_ch4_ = plist_->get("atmospheric CH4"); + atm_n2o_ = plist_->get("atmospheric N2O"); + atm_h2_ = plist_->get("atmospheric H2"); + atm_nh3_ = plist_->get("atmospheric NH3");*/ + + //Starting values and parameters for precribed phenology / albedo + + pressure_at_field_capacity = plist_->get("field capacity [Mpa]"); + pressure_at_wilting_point = plist_->get("wilting point [Mpa]"); + p_bool = plist_->get("EcoSIM precipitation"); + a_bool = plist_->get("prescribe snow albedo"); + pheno_bool = plist_->get("prescribe phenology"); + + //Parameters for times and time of year + dt_ = plist_->get("initial time step"); + c_m_ = plist_->get("heat capacity [MJ mol^-1 K^-1]"); + day0_ = plist_->get("starting day of year [0-364]"); + year0_ = plist_->get("starting year"); + + curr_day_ = day0_; + curr_year_ = year0_; + + //This initializes the engine (found in BGCEngine.cc) This is the code that + //actually points to the driver + if (!plist_->isParameter("engine")) { + Errors::Message msg; + msg << "No 'engine' parameter found in the parameter list for 'BGC'.\n"; + Exceptions::amanzi_throw(msg); + } + if (!plist_->isParameter("engine input file")) { + Errors::Message msg; + msg << "No 'engine input file' parameter found in the parameter list for 'BGC'.\n"; + Exceptions::amanzi_throw(msg); + } + std::string engine_name = plist_->get("engine"); + std::string engine_inputfile = plist_->get("engine input file"); + bgc_engine_ = Teuchos::rcp(new BGCEngine(engine_name, engine_inputfile)); + } + + +// -- Destroy ansilary data structures. +EcoSIM::~EcoSIM() + { + if (bgc_initialized_) + bgc_engine_->FreeState(bgc_props_, bgc_state_, bgc_aux_data_); + } + +// -- Setup step +void EcoSIM::Setup() { + PK_Physical_Default::Setup(); + //Need to do some basic setup of the columns: + mesh_surf_ = S_->GetMesh(domain_surface_); + mesh_ = S_->GetMesh(domain_); + int num_columns_ = + mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + + for (unsigned int column = 0; column != num_columns_; ++column) { + int f = mesh_surf_->getEntityParent(AmanziMesh::Entity_kind::CELL, column); + auto col_iter = mesh_->columns.getCells(column); + std::size_t ncol_cells = col_iter.size(); + + double column_area = mesh_->getFaceArea(f); + + if (ncells_per_col_ < 0) { + ncells_per_col_ = ncol_cells; + } else { + AMANZI_ASSERT(ncol_cells == ncells_per_col_); + } + } + + //Setting records for variables ONLY used by the EcoSIM PK, this includes, weather forcings, + // ecosim surface and subsurface forces, and surface variables from EcoSIM that are saved + // over to ATS (evaporation fluxes, snow related variables) + if (!S_->HasRecord(snow_depth_key_,tag_next_)) { + S_->Require(snow_depth_key_, tag_next_, snow_depth_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + } + + S_->Require(canopy_lw_key_ , tag_next_, canopy_lw_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(canopy_latent_heat_key_ , tag_next_, canopy_latent_heat_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(canopy_sensible_heat_key_, tag_next_, canopy_sensible_heat_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(canopy_surface_water_key_ , tag_next_, canopy_surface_water_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(transpiration_key_ , tag_next_, transpiration_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(evaporation_canopy_key_ , tag_next_, evaporation_canopy_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(evaporation_ground_key_ , tag_next_, evaporation_ground_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(evaporation_litter_key_ , tag_next_, evaporation_litter_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(evaporation_snow_key_ , tag_next_, evaporation_snow_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(sublimation_snow_key_ , tag_next_, sublimation_snow_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(surface_energy_source_ecosim_key_ , tag_next_, name_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(surface_water_source_ecosim_key_ , tag_next_, surface_water_source_ecosim_key_) + .SetMesh(mesh_surf_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(subsurface_energy_source_ecosim_key_ , tag_next_, subsurface_energy_source_ecosim_key_) + .SetMesh(mesh_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(subsurface_water_source_ecosim_key_ , tag_next_, subsurface_water_source_ecosim_key_) + .SetMesh(mesh_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->Require(snow_temperature_key_ , tag_next_, snow_temperature_key_) + .SetMesh(mesh_) + ->SetGhosted(false) + ->SetComponent("cell", AmanziMesh::CELL, 1); + + S_->RequireEvaluator(snow_albedo_key_, tag_next_); + S_->Require(snow_albedo_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + + S_->RequireEvaluator(sw_key_, tag_next_); + S_->Require(sw_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + + S_->RequireEvaluator(lai_key_, tag_next_); + S_->Require(lai_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + + S_->RequireEvaluator(sai_key_, tag_next_); + S_->Require(sai_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + + S_->RequireEvaluator(v_type_key_, tag_next_); + S_->Require(v_type_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + + Teuchos::OSTab tab = vo_->getOSTab(); + + //EcoSIM can do its own precipitation partitioning so you can put in total precipitation if + // you want EcoSIM to do it, or snow/rain if the forcing is already split. + if (p_bool) { + S_->RequireEvaluator(p_total_key_, tag_next_); + S_->Require(p_total_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + } else { + S_->RequireEvaluator(p_snow_key_, tag_next_); + S_->Require(p_snow_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + S_->RequireEvaluator(p_rain_key_, tag_next_); + S_->Require(p_rain_key_, tag_next_).SetMesh(mesh_surf_) + ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1); + } + + //Setup custom evaluators for EcoSIM, found in constitutive relations + requireEvaluatorAtNext(hydraulic_conductivity_key_, tag_next_, *S_) + .SetMesh(mesh_) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + + requireEvaluatorAtCurrent(hydraulic_conductivity_key_, tag_current_, *S_, name_); + +//Setup variables that were owned by ATS SEB, but now are controlled by EcoSIM +// Can remove SEB from the cycle_driver + requireEvaluatorAtNext(lw_key_, tag_next_, *S_) + .SetMesh(mesh_surf_) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + + requireEvaluatorAtCurrent(lw_key_, tag_current_, *S_, name_); + + requireEvaluatorAtNext(air_temp_key_, tag_next_, *S_) + .SetMesh(mesh_surf_) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + + requireEvaluatorAtCurrent(air_temp_key_, tag_current_, *S_, name_); + + requireEvaluatorAtNext(vp_air_key_, tag_next_, *S_) + .SetMesh(mesh_surf_) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + + requireEvaluatorAtCurrent(vp_air_key_, tag_current_, *S_, name_); + + requireEvaluatorAtNext(wind_speed_key_, tag_next_, *S_) + .SetMesh(mesh_surf_) + ->SetGhosted() + ->AddComponent("cell", AmanziMesh::CELL, 1); + + requireEvaluatorAtCurrent(wind_speed_key_, tag_current_, *S_, name_); + + if (vo_->os_OK(Teuchos::VERB_MEDIUM)) { + Teuchos::OSTab tab = vo_->getOSTab(); + *vo_->os() << vo_->color("green") << "Setup of PK was successful" + << vo_->reset() << std::endl << std::endl; + } +} + +// -- Initialize owned (dependent) variables. +void EcoSIM::Initialize() { + PK_Physical_Default::Initialize(); + //Need to know the number of components to initialize data structures + + //Transport removal: + /*const Epetra_MultiVector& mole_fraction= *(S_->GetPtr(mole_fraction_key_, Tags::DEFAULT)->ViewComponent("cell")); + int mole_fraction_num = mole_fraction.NumVectors();*/ + int mole_fraction_num = 1; + Teuchos::OSTab tab = vo_->getOSTab(); + + num_columns_ = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + + //Now we call the engine's init state function which allocates the data + bgc_engine_->InitState(bgc_props_, bgc_state_, bgc_aux_data_, ncells_per_col_, mole_fraction_num, num_columns_); + + int ierr = 0; + + if (S_->HasRecord(ice_density_key_, Tags::DEFAULT)) { + S_->GetEvaluator(saturation_ice_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(ice_density_key_, Tags::DEFAULT).Update(*S_, name_); + has_ice = true; + } else { + Teuchos::OSTab tab = vo_->getOSTab(); + //*vo_->os() << "Did not find ice key" << std::endl; + has_ice = false; + } + + //Check for total precipitation and set the record if it's there + + if (p_bool) { + S_->GetW(p_total_key_, Tags::DEFAULT, "surface-precipitation_total").PutScalar(0.0); + S_->GetRecordW(p_total_key_, Tags::DEFAULT, "surface-precipitation_total").set_initialized(); + } else { + S_->GetW(p_snow_key_, Tags::DEFAULT, "surface-precipitation_snow").PutScalar(0.0); + S_->GetRecordW(p_snow_key_, Tags::DEFAULT, "surface-precipitation_snow").set_initialized(); + S_->GetW(p_rain_key_, Tags::DEFAULT, "surface-precipitation_rain").PutScalar(0.0); + S_->GetRecordW(p_rain_key_, Tags::DEFAULT, "surface-precipitation_rain").set_initialized(); + } + + S_->GetW(snow_depth_key_, Tags::DEFAULT, "surface-snow_depth").PutScalar(0.0); + S_->GetRecordW(snow_depth_key_, Tags::DEFAULT, "surface-snow_depth").set_initialized(); + + S_->GetW(canopy_lw_key_, Tags::DEFAULT, "surface-canopy_longwave_radiation").PutScalar(0.0); + S_->GetRecordW(canopy_lw_key_, Tags::DEFAULT, "surface-canopy_longwave_radiation").set_initialized(); + + S_->GetW(canopy_latent_heat_key_, Tags::DEFAULT, "surface-canopy_latent_heat").PutScalar(0.0); + S_->GetRecordW(canopy_latent_heat_key_, Tags::DEFAULT, "surface-canopy_latent_heat").set_initialized(); + + S_->GetW(canopy_sensible_heat_key_, Tags::DEFAULT, "surface-canopy_sensible_heat").PutScalar(0.0); + S_->GetRecordW(canopy_sensible_heat_key_, Tags::DEFAULT, "surface-canopy_sensible_heat").set_initialized(); + + S_->GetW(canopy_surface_water_key_, Tags::DEFAULT, "surface-canopy_surface_water").PutScalar(0.0); + S_->GetRecordW(canopy_surface_water_key_, Tags::DEFAULT, "surface-canopy_surface_water").set_initialized(); + + S_->GetW(transpiration_key_, Tags::DEFAULT, "surface-transpiration").PutScalar(0.0); + S_->GetRecordW(transpiration_key_, Tags::DEFAULT, "surface-transpiration").set_initialized(); + + S_->GetW(evaporation_canopy_key_, Tags::DEFAULT, "surface-evaporation_canopy").PutScalar(0.0); + S_->GetRecordW(evaporation_canopy_key_, Tags::DEFAULT, "surface-evaporation_canopy").set_initialized(); + + S_->GetW(evaporation_ground_key_, Tags::DEFAULT, "surface-evaporation_ground").PutScalar(0.0); + S_->GetRecordW(evaporation_ground_key_, Tags::DEFAULT, "surface-evaporation_ground").set_initialized(); + + S_->GetW(evaporation_litter_key_, Tags::DEFAULT, "surface-evaporation_litter").PutScalar(0.0); + S_->GetRecordW(evaporation_litter_key_, Tags::DEFAULT, "surface-evaporation_litter").set_initialized(); + + S_->GetW(evaporation_snow_key_, Tags::DEFAULT, "surface-evaporation_snow").PutScalar(0.0); + S_->GetRecordW(evaporation_snow_key_, Tags::DEFAULT, "surface-evaporation_snow").set_initialized(); + + S_->GetW(sublimation_snow_key_, Tags::DEFAULT, "surface-sublimation_snow").PutScalar(0.0); + S_->GetRecordW(sublimation_snow_key_, Tags::DEFAULT, "surface-sublimation_snow").set_initialized(); + + S_->GetW(surface_water_source_ecosim_key_, Tags::DEFAULT, surface_water_source_ecosim_key_).PutScalar(0.0); + S_->GetRecordW(surface_water_source_ecosim_key_, Tags::DEFAULT, surface_water_source_ecosim_key_).set_initialized(); + + S_->GetW(subsurface_water_source_ecosim_key_, Tags::DEFAULT, subsurface_water_source_ecosim_key_).PutScalar(0.0); + S_->GetRecordW(subsurface_water_source_ecosim_key_, Tags::DEFAULT, subsurface_water_source_ecosim_key_).set_initialized(); + + S_->GetW(subsurface_energy_source_ecosim_key_, Tags::DEFAULT, subsurface_energy_source_ecosim_key_).PutScalar(0.0); + S_->GetRecordW(subsurface_energy_source_ecosim_key_, Tags::DEFAULT, subsurface_energy_source_ecosim_key_).set_initialized(); + + S_->GetW(snow_temperature_key_, Tags::DEFAULT, snow_temperature_key_).PutScalar(0.0); + S_->GetRecordW(snow_temperature_key_, Tags::DEFAULT, snow_temperature_key_).set_initialized(); + + //Initialize owned evaluators + S_->GetW(hydraulic_conductivity_key_, Tags::DEFAULT, "hydraulic_conductivity").PutScalar(1.0); + S_->GetRecordW(hydraulic_conductivity_key_, Tags::DEFAULT, "hydraulic_conductivity").set_initialized(); + + //S_->GetW(bulk_density_key_, Tags::DEFAULT, "bulk_density").PutScalar(1.0); + //S_->GetRecordW(bulk_density_key_, Tags::DEFAULT, "bulk_density").set_initialized(); + + //S_->GetW(matric_pressure_key_, Tags::DEFAULT, "matric_pressure").PutScalar(1.0); + //S_->GetRecordW(matric_pressure_key_, Tags::DEFAULT, "matric_pressure").set_initialized(); + + int num_columns_ = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + + //loop over processes instead: + num_columns_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL, false).NumGlobalElements(); + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + num_columns_global_ptype = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::ALL); + + //Loop over processes and Initalize EcoSIM on that process + int numProcesses, p_rank; + MPI_Comm_size(MPI_COMM_WORLD, &numProcesses); + MPI_Comm_rank(MPI_COMM_WORLD, &p_rank); + for (int k = 0; k < numProcesses; ++k) { + MPI_Barrier(MPI_COMM_WORLD); + if (p_rank==k) { + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + + InitializeSingleProcess(p_rank); + } + } + + // verbose message + if (vo_->os_OK(Teuchos::VERB_MEDIUM)) { + Teuchos::OSTab tab = vo_->getOSTab(); + *vo_->os() << vo_->color("green") << "Initialization of PK was successful, T=" + << S_->get_time() << vo_->reset() << std::endl << std::endl; + } +} + +void EcoSIM::CommitStep(double t_old, double t_new, const Tag& tag) { + + // I don't know that we will have much to do here. In SimpleBGC they just copy + // Data to the pfts, which we won't be doing. In Alquimia they just save the time + // As below. + + saved_time_ = t_new; + +} + +bool EcoSIM::AdvanceStep(double t_old, double t_new, bool reinit) { + double dt = t_new - t_old; + current_time_ = saved_time_ + dt; + + + Teuchos::OSTab out = vo_->getOSTab(); + if (vo_->os_OK(Teuchos::VERB_HIGH)) + *vo_->os() << "----------------------------------------------------------------" << std::endl + << "Advancing: t0 = " << S_->get_time(tag_current_) + << " t1 = " << S_->get_time(tag_next_) << " h = " << dt << std::endl + << "Current day: " << curr_day_ << "Current year: " << curr_year_ << std::endl + << "----------------------------------------------------------------" << std::endl; + + + // Ensure dependencies are filled + // Transport removal + + //S_->GetEvaluator(mole_fraction_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(porosity_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(saturation_liquid_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(water_content_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(relative_permeability_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(liquid_density_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(rock_density_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(T_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(cell_volume_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(f_wp_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(f_root_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(cap_pres_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(subsurface_energy_source_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(subsurface_water_source_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(matric_pressure_key_, Tags::DEFAULT).Update(*S_, name_); + + + //Surface data + S_->GetEvaluator(sw_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(lw_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(air_temp_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(vp_air_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(wind_speed_key_, Tags::DEFAULT).Update(*S_, name_); + + S_->GetEvaluator(elev_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(aspect_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(slope_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(snow_albedo_key_, Tags::DEFAULT).Update(*S_, name_); + + S_->GetEvaluator(lai_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(sai_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(v_type_key_, Tags::DEFAULT).Update(*S_, name_); + + if (p_bool){ + S_->GetEvaluator(p_total_key_, Tags::DEFAULT).Update(*S_, name_); + } else { + S_->GetEvaluator(p_rain_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(p_snow_key_, Tags::DEFAULT).Update(*S_, name_); + } + + //S_->GetEvaluator(surface_energy_source_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(surface_water_source_key_, Tags::DEFAULT).Update(*S_, name_); + + if (has_gas) { + S_->GetEvaluator(saturation_gas_key_, Tags::DEFAULT).Update(*S_, name_); + //S_->GetEvaluator(gas_density_key_, Tags::DEFAULT).Update(*S_, name_); + } + + if (has_ice) { + S_->GetEvaluator(saturation_ice_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(ice_density_key_, Tags::DEFAULT).Update(*S_, name_); + } + + S_->GetEvaluator(T_key_, Tags::DEFAULT).Update(*S_, name_); + S_->GetEvaluator(thermal_conductivity_key_, Tags::DEFAULT).Update(*S_, name_); + + //Update owned evaluators + /*Teuchos::RCP hydra_cond = S_->GetPtr(hydraulic_conductivity_key_, Tags::DEFAULT); + S_->GetEvaluator(hydraulic_conductivity_key_, Tags::DEFAULT).Update(*S_, name_); + const Epetra_MultiVector& hydraulic_conductivity = *(*S_->Get("hydraulic_conductivity", tag_next_) + .ViewComponent("cell",false))(0);*/ + + AmanziMesh::Entity_ID num_columns_ = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + + // grab the required fields + + S_->GetEvaluator("porosity", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& porosity = *(*S_->Get("porosity", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("saturation_liquid", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& liquid_saturation = *(*S_->Get("saturation_liquid", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("capillary_pressure_gas_liq", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& capillary_pressure = *(*S_->Get("capillary_pressure_gas_liq", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("water_content", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& water_content = *(*S_->Get("water_content", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("mass_density_liquid", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& liquid_density = *(*S_->Get("mass_density_liquid", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("density_rock", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& rock_density = *(*S_->Get("density_rock", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("cell_volume", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& cell_volume = *(*S_->Get("cell_volume", tag_next_) + .ViewComponent("cell",false))(0); + + if (has_gas) { + S_->GetEvaluator("saturation_gas", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& gas_saturation = *(*S_->Get("saturation_gas", tag_next_) + .ViewComponent("cell",false))(0); + } + + //Atm abundances + S_->GetEvaluator("surface-incoming_shortwave_radiation", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& sw_rad = *(*S_->Get("surface-incoming_shortwave_radiation", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("surface-incoming_longwave_radiation", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& lw_rad = *(*S_->Get("surface-incoming_longwave_radiation", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("surface-air_temperature", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& t_air = *(*S_->Get("surface-air_temperature", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("surface-vapor_pressure_air", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& p_vap = *(*S_->Get("surface-vapor_pressure_air", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("surface-wind_speed", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& v_wind = *(*S_->Get("surface-wind_speed", tag_next_) + .ViewComponent("cell",false))(0); + + //Define before loop to prevent scope issues: + const Epetra_MultiVector* p_tot = nullptr; + const Epetra_MultiVector* p_rain = nullptr; + const Epetra_MultiVector* p_snow = nullptr; + + if(p_bool){ + S_->GetEvaluator("surface-precipitation_total", tag_next_).Update(*S_, name_); + p_tot = &(*(*S_->Get("surface-precipitation_total", tag_next_) + .ViewComponent("cell",false))(0)); + } else { + S_->GetEvaluator("surface-precipitation_rain", tag_next_).Update(*S_, name_); + p_rain = &(*(*S_->Get("surface-precipitation_rain", tag_next_) + .ViewComponent("cell",false))(0)); + S_->GetEvaluator("surface-precipitation_snow", tag_next_).Update(*S_, name_); + p_snow = &(*(*S_->Get("surface-precipitation_snow", tag_next_) + .ViewComponent("cell",false))(0)); + } + + S_->GetEvaluator("surface-elevation", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& elevation = *S_->Get("surface-elevation", tag_next_) + .ViewComponent("cell",false); + + S_->GetEvaluator("surface-aspect", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& aspect = *S_->Get("surface-aspect", tag_next_) + .ViewComponent("cell",false); + + S_->GetEvaluator("surface-slope_magnitude", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& slope = *S_->Get("surface-slope_magnitude", tag_next_) + .ViewComponent("cell",false); + + S_->GetEvaluator("surface-LAI", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& LAI = *(*S_->Get("surface-LAI", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("surface-SAI", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& SAI = *(*S_->Get("surface-SAI", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("surface-vegetation_type", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& vegetation_type = *(*S_->Get("surface-vegetation_type", tag_next_) + .ViewComponent("cell",false))(0); + + if (has_ice) { + S_->GetEvaluator("mass_density_ice", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& ice_density = *(*S_->Get("mass_density_ice", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("saturation_ice", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& ice_saturation = *(*S_->Get("saturation_ice", tag_next_) + .ViewComponent("cell",false))(0); + } + + S_->GetEvaluator("temperature", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& temp = *(*S_->Get("temperature", tag_next_) + .ViewComponent("cell",false))(0); + + S_->GetEvaluator("thermal_conductivity", tag_next_).Update(*S_, name_); + const Epetra_MultiVector& thermal_conductivity = *(*S_->Get("thermal_conductivity", tag_next_) + .ViewComponent("cell",false))(0); + + //loop over processes instead: + num_columns_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL,false).NumGlobalElements(); + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + num_columns_global_ptype = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::ALL); + + //Trying to loop over processors now: + int numProcesses, p_rank; + MPI_Comm_size(MPI_COMM_WORLD, &numProcesses); + MPI_Comm_rank(MPI_COMM_WORLD, &p_rank); + for (int k = 0; k < numProcesses; ++k) { + MPI_Barrier(MPI_COMM_WORLD); + if (p_rank==k) { + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + + AdvanceSingleProcess(dt, p_rank); + } + } + // PLACE TIME AND YEAR ITERATOR HERE + if (curr_day_ != 364) { + curr_day_ = curr_day_ + 1; + } else { + curr_day_ = 0; + curr_year_ = curr_year_ + 1; + } + +} + +// helper function for pushing field to column +void EcoSIM::FieldToColumn_(AmanziMesh::Entity_ID column, const Epetra_Vector& vec, + Teuchos::Ptr col_vec) +{ + auto col_iter = mesh_->columns.getCells(column); + + for (std::size_t i=0; i!=col_iter.size(); ++i) { + std::size_t vec_index = col_iter[i]; + + (*col_vec)[i] = vec[vec_index]; + } +} + +void EcoSIM::FieldToColumn_(AmanziMesh::Entity_ID column, const Teuchos::Ptr vec, + Teuchos::Ptr col_vec) +{ + auto col_iter = mesh_->columns.getCells(column); + + for (std::size_t i=0; i!=col_iter.size(); ++i) { + std::size_t vec_index = col_iter[i]; + + (*col_vec)[i] = (*vec)[vec_index]; + } +} + +//Helper function but for datasets that are multivalued in every cell (concentrations) +void EcoSIM::MatrixFieldToColumn_(AmanziMesh::Entity_ID column, const Epetra_MultiVector& m_arr, + Teuchos::Ptr col_arr) + { + int n_comp = m_arr.NumVectors(); + auto col_iter = mesh_->columns.getCells(column); + + for (int j=0; j!=n_comp; ++j){ + for (std::size_t i=0; i!=col_iter.size(); ++i) { + (*col_arr)(i,j) = m_arr[j][col_iter[i]]; + } + } + } + +// helper function for pushing column back to field +void EcoSIM::ColumnToField_(AmanziMesh::Entity_ID column, Epetra_Vector& vec, + Teuchos::Ptr col_vec) +{ + auto col_iter = mesh_->columns.getCells(column); + for (std::size_t i=0; i!=col_iter.size(); ++i) { + vec[col_iter[i]] = (*col_vec)[i]; + } +} + +void EcoSIM::ColumnToField_(AmanziMesh::Entity_ID column, Teuchos::Ptr vec, + Teuchos::Ptr col_vec) +{ + auto col_iter = mesh_->columns.getCells(column); + for (std::size_t i=0; i!=col_iter.size(); ++i) { + (*vec)[col_iter[i]] = (*col_vec)[i]; + } +} + +void EcoSIM::MatrixColumnToField_(AmanziMesh::Entity_ID column, Epetra_MultiVector& m_arr, + Teuchos::Ptr col_arr) { + + int n_comp = m_arr.NumVectors(); + auto col_iter = mesh_->columns.getCells(column); + + for (int j=0; j!=n_comp; ++j){ + for (std::size_t i=0; i!=col_iter.size(); ++i) { + m_arr[j][col_iter[i]] = (*col_arr)(i,j); + } + } + + } + +// helper function for collecting column dz and depth +void EcoSIM::ColDepthDz_(AmanziMesh::Entity_ID column, + Teuchos::Ptr depth, + Teuchos::Ptr dz) { + AmanziMesh::Entity_ID f_above = mesh_surf_->getEntityParent(AmanziMesh::Entity_kind::CELL, column); + auto col_iter = mesh_->columns.getCells(column); + ncells_per_col_ = col_iter.size(); + + AmanziGeometry::Point surf_centroid = mesh_->getFaceCentroid(f_above); + AmanziGeometry::Point neg_z(3); + neg_z.set(0.,0.,-1); + + Teuchos::OSTab tab = vo_->getOSTab(); + + for (std::size_t i=0; i!=col_iter.size(); ++i) { + // depth centroid + (*depth)[i] = surf_centroid[2] - mesh_->getCellCentroid(col_iter[i])[2]; + + const auto& [faces, dirs] = mesh_->getCellFacesAndDirections(col_iter[i]); + + // -- mimics implementation of build_columns() in Mesh + double mindp = 999.0; + AmanziMesh::Entity_ID f_below = -1; + for (std::size_t j=0; j!=faces.size(); ++j) { + AmanziGeometry::Point normal = mesh_->getFaceNormal(faces[j]); + if (dirs[j] == -1) normal *= -1; + normal /= AmanziGeometry::norm(normal); + + double dp = -normal * neg_z; + if (dp < mindp) { + mindp = dp; + f_below = faces[j]; + } + } + + // -- fill the val + (*dz)[i] = mesh_->getFaceCentroid(f_above)[2] - mesh_->getFaceCentroid(f_below)[2]; + AMANZI_ASSERT( (*dz)[i] > 0. ); + f_above = f_below; + } +} + +// helper function for collecting dz, depth, and volume for a given column +void EcoSIM::VolDepthDz_(AmanziMesh::Entity_ID column, + Teuchos::Ptr depth, + Teuchos::Ptr dz, + Teuchos::Ptr volume) { + AmanziMesh::Entity_ID f_above = mesh_surf_->getEntityParent(AmanziMesh::Entity_kind::CELL, column); + auto col_iter = mesh_->columns.getCells(column); + ncells_per_col_ = col_iter.size(); + + AmanziGeometry::Point surf_centroid = mesh_->getFaceCentroid(f_above); + AmanziGeometry::Point neg_z(3); + neg_z.set(0.,0.,-1); + + for (std::size_t i=0; i!=col_iter.size(); ++i) { + // depth centroid + (*depth)[i] = surf_centroid[2] - mesh_->getCellCentroid(col_iter[i])[2]; + + // dz + // -- find face_below + //AmanziMesh::Entity_ID_List faces; + //std::vector dirs; + //mesh_->cell_get_faces_and_dirs(col_iter[i], &faces, &dirs); + + const auto& [faces, dirs] = mesh_->getCellFacesAndDirections(col_iter[i]); + + //double vol = mesh_->cell_volume(col_iter[i]); + (*volume)[i] = mesh_->getCellVolume(col_iter[i]); + + // -- mimics implementation of build_columns() in Mesh + double mindp = 999.0; + AmanziMesh::Entity_ID f_below = -1; + for (std::size_t j=0; j!=faces.size(); ++j) { + AmanziGeometry::Point normal = mesh_->getFaceNormal(faces[j]); + if (dirs[j] == -1) normal *= -1; + normal /= AmanziGeometry::norm(normal); + + double dp = -normal * neg_z; + if (dp < mindp) { + mindp = dp; + f_below = faces[j]; + } + } + + // -- fill the val + (*dz)[i] = mesh_->getFaceCentroid(f_above)[2] - mesh_->getFaceCentroid(f_below)[2]; + AMANZI_ASSERT( (*dz)[i] > 0. ); + f_above = f_below; + } +} + +//Copy to EcoSIM +void EcoSIM::CopyToEcoSIM_process(int proc_rank, + BGCProperties& props, + BGCState& state, + BGCAuxiliaryData& aux_data, + const Tag& water_tag) +{ + //This is the copy function for a loop over a single process instead of a single column + //Fill state with ATS variables that are going to be changed by EcoSIM + const Epetra_Vector& porosity = *(*S_->Get(porosity_key_, water_tag).ViewComponent("cell", false))(0); + + //Transport removal + /*const Epetra_MultiVector& mole_fraction= *(S_->GetPtr(mole_fraction_key_, water_tag)->ViewComponent("cell")); + int mole_fraction_num = mole_fraction.NumVectors();*/ + int mole_fraction_num = 1; + + const Epetra_Vector& liquid_saturation = *(*S_->Get(saturation_liquid_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& water_content = *(*S_->Get(water_content_key_, water_tag).ViewComponent("cell", false))(0); + //const Epetra_Vector& relative_permeability = *(*S_->Get(relative_permeability_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& liquid_density = *(*S_->Get(liquid_density_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& rock_density = *(*S_->Get(rock_density_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& cell_volume = *(*S_->Get(cell_volume_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& hydraulic_conductivity = *(*S_->Get(hydraulic_conductivity_key_, water_tag).ViewComponent("cell", false))(0); + //const Epetra_Vector& matric_pressure = *(*S_->Get(matric_pressure_key_, water_tag).ViewComponent("cell", false))(0); + //const Epetra_Vector& bulk_density = *(*S_->Get(bulk_density_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& rooting_depth_fraction = *(*S_->Get(f_root_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& plant_wilting_factor = *(*S_->Get(f_wp_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& temp = *(*S_->Get(T_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& thermal_conductivity = *(*S_->Get(thermal_conductivity_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& capillary_pressure = *(*S_->Get(cap_pres_key_, water_tag).ViewComponent("cell", false))(0); + + //const auto& shortwave_radiation = *S_.Get(sw_key_, water_tag).ViewComponent("cell", false); + const Epetra_Vector& shortwave_radiation = *(*S_->Get(sw_key_, water_tag).ViewComponent("cell", false))(0); + //const Epetra_Vector& longwave_radiation = *(*S_->Get(lw_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& air_temperature = *(*S_->Get(air_temp_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& vapor_pressure_air = *(*S_->Get(vp_air_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& wind_speed= *(*S_->Get(wind_speed_key_, water_tag).ViewComponent("cell", false))(0); + //define these outside of the loop to prevent issues: + const Epetra_Vector* precipitation = nullptr; + const Epetra_Vector* precipitation_snow = nullptr; + + if(p_bool){ + precipitation = &(*(*S_->Get(p_total_key_, water_tag).ViewComponent("cell", false))(0)); + } else { + precipitation = &(*(*S_->Get(p_rain_key_, water_tag).ViewComponent("cell", false))(0)); + precipitation_snow = &(*(*S_->Get(p_snow_key_, water_tag).ViewComponent("cell", false))(0)); + } + + const Epetra_Vector& elevation = *(*S_->Get(elev_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& aspect = *(*S_->Get(aspect_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& slope = *(*S_->Get(slope_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& snow_albedo = *(*S_->Get(snow_albedo_key_, water_tag).ViewComponent("cell", false))(0); + + const Epetra_Vector& LAI = *(*S_->Get(lai_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& SAI = *(*S_->Get(sai_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& vegetation_type = *(*S_->Get(v_type_key_, water_tag).ViewComponent("cell", false))(0); + + const Epetra_Vector& surface_energy_source = *(*S_->Get(surface_energy_source_ecosim_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& subsurface_energy_source = *(*S_->Get(subsurface_energy_source_ecosim_key_, water_tag).ViewComponent("cell", false))(0); + + const Epetra_Vector& surface_water_source = *(*S_->Get(surface_water_source_ecosim_key_, water_tag).ViewComponent("cell", false))(0); + const Epetra_Vector& subsurface_water_source = *(*S_->Get(subsurface_water_source_ecosim_key_, water_tag).ViewComponent("cell", false))(0); + + auto& snow_depth = *S_->GetW(snow_depth_key_,tag_next_,snow_depth_key_).ViewComponent("cell"); + + auto& canopy_longwave_radiation = *S_->GetW(canopy_lw_key_, tag_next_, canopy_lw_key_).ViewComponent("cell"); + auto& canopy_latent_heat = *S_->GetW(canopy_latent_heat_key_, tag_next_, canopy_latent_heat_key_).ViewComponent("cell"); + auto& canopy_sensible_heat = *S_->GetW(canopy_sensible_heat_key_, tag_next_, canopy_sensible_heat_key_).ViewComponent("cell"); + auto& canopy_surface_water = *S_->GetW(canopy_surface_water_key_, tag_next_, canopy_surface_water_key_).ViewComponent("cell"); + auto& transpiration = *S_->GetW(transpiration_key_, tag_next_, transpiration_key_).ViewComponent("cell"); + auto& evaporation_canopy = *S_->GetW(evaporation_canopy_key_, tag_next_, evaporation_canopy_key_).ViewComponent("cell"); + auto& evaporation_ground = *S_->GetW(evaporation_ground_key_, tag_next_, evaporation_ground_key_).ViewComponent("cell"); + auto& evaporation_litter = *S_->GetW(evaporation_litter_key_, tag_next_, evaporation_litter_key_).ViewComponent("cell"); + auto& evaporation_snow = *S_->GetW(evaporation_snow_key_, tag_next_, evaporation_snow_key_).ViewComponent("cell"); + auto& sublimation_snow = *S_->GetW(sublimation_snow_key_, tag_next_, sublimation_snow_key_).ViewComponent("cell"); + + auto col_porosity = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_l_sat = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_l_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_wc = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_relative_permeability = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_mat_p = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_r_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_vol = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_g_sat = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_g_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_i_sat = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_i_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_temp = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_cond = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_h_cond = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_b_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_depth = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_dz = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_wp = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_rf = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_lai = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_sai = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_v_type = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_ss_energy_source = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_ss_water_source = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_depth_c = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_cap_pres = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + + auto col_mole_fraction = Teuchos::rcp(new Epetra_SerialDenseMatrix(mole_fraction_num,ncells_per_col_)); + + //Gather columns on this process: + num_columns_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL,false).NumGlobalElements(); + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + num_columns_global_ptype = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::ALL); + + //Trying to loop over processors now: + int p_rank; + MPI_Comm_rank(MPI_COMM_WORLD, &p_rank); + MPI_Barrier(MPI_COMM_WORLD); + + std::cout << "ATS2EcoSIM rank: " << p_rank <getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + //Now that the arrays are flat we need to be a little more careful about how we load an unload the data + /*for (int column=0; column!=num_columns_local; ++column) { + FieldToColumn_(column, temp, col_temp.ptr()); + + for (int i=0; i < ncells_per_col_; ++i) { + state.temperature.data[column * ncells_per_col_ + i] = (*col_temp)[i]; + state.temperature.data[column * ncells_per_col_ + i] = 222.0; + } + }*/ + + //Loop over columns on this process + for (int column=0; column!=num_columns_local; ++column) { + FieldToColumn_(column,porosity,col_porosity.ptr()); + FieldToColumn_(column,liquid_saturation,col_l_sat.ptr()); + FieldToColumn_(column,water_content,col_wc.ptr()); + //FieldToColumn_(column,relative_permeability,col_relative_permeability.ptr()); + FieldToColumn_(column,liquid_density,col_l_dens.ptr()); + FieldToColumn_(column,rock_density,col_r_dens.ptr()); + FieldToColumn_(column,cell_volume,col_vol.ptr()); + FieldToColumn_(column,hydraulic_conductivity,col_h_cond.ptr()); + //FieldToColumn_(column,bulk_density,col_b_dens.ptr()); + FieldToColumn_(column,plant_wilting_factor,col_wp.ptr()); + FieldToColumn_(column,rooting_depth_fraction,col_rf.ptr()); + FieldToColumn_(column,subsurface_water_source,col_ss_water_source.ptr()); + FieldToColumn_(column,subsurface_energy_source,col_ss_energy_source.ptr()); + //setting matric pressure to capillary pressure here + FieldToColumn_(column,capillary_pressure,col_mat_p.ptr()); + FieldToColumn_(column,temp, col_temp.ptr()); + FieldToColumn_(column,thermal_conductivity,col_cond.ptr()); + //FieldToColumn_(column,capillary_pressure,col_cap_pres.ptr()); + + //MatrixFieldToColumn_(column, mole_fraction, col_mole_fraction.ptr()); + + // This is for computing depth + //ColDepthDz_(column, col_depth.ptr(), col_dz.ptr()); + + //Grabbing the cross sectional area in the z direction + int f = mesh_surf_->getEntityParent(AmanziMesh::Entity_kind::CELL, column); + auto col_iter = mesh_->columns.getCells(column); + std::size_t ncol_cells = col_iter.size(); + + double column_area = mesh_->getFaceArea(f); + //std::cout << "column: " << column << " column_area: " << column_area << std::endl; + props.column_area.data[column] = column_area; + + VolDepthDz_(column, col_depth.ptr(), col_dz.ptr(), col_vol.ptr()); + double sum = 0.0; + for (int i = ncells_per_col_ - 1; i >= 0; --i) { + sum += (*col_dz)[i]; + (*col_depth_c)[i] = sum; + } + + for (int i=0; i < ncells_per_col_; ++i) { + state.liquid_density.data[column * ncells_per_col_ + i] = (*col_l_dens)[i]; + state.rock_density.data[column * ncells_per_col_ + i] = (*col_r_dens)[i]; + state.porosity.data[column * ncells_per_col_ + i] = (*col_porosity)[i]; + state.water_content.data[column * ncells_per_col_ + i] = (*col_wc)[i]; + state.hydraulic_conductivity.data[column * ncells_per_col_ + i] = (*col_h_cond)[i]; + //state.bulk_density.data[column * ncells_per_col_ + i] = (*col_b_dens)[i]; + state.subsurface_water_source.data[column * ncells_per_col_ + i] = (*col_ss_water_source)[i]; + state.subsurface_energy_source.data[column * ncells_per_col_ + i] = (*col_ss_energy_source)[i]; + state.matric_pressure.data[column * ncells_per_col_ + i] = (*col_mat_p)[i]; + state.temperature.data[column * ncells_per_col_ + i] = (*col_temp)[i]; + + props.plant_wilting_factor.data[column * ncells_per_col_ + i] = (*col_wp)[i]; + props.rooting_depth_fraction.data[column * ncells_per_col_ + i] = (*col_rf)[i]; + props.liquid_saturation.data[column * ncells_per_col_ + i] = (*col_l_sat)[i]; + props.relative_permeability.data[column * ncells_per_col_ + i] = (*col_relative_permeability)[i]; + props.volume.data[column * ncells_per_col_ + i] = (*col_vol)[i]; + props.depth.data[column * ncells_per_col_ + i] = (*col_depth)[i]; + props.depth_c.data[column * ncells_per_col_ + i] = (*col_depth_c)[i]; + props.dz.data[column * ncells_per_col_ + i] = (*col_dz)[i]; + + if (has_gas) { + props.gas_saturation.data[column * ncells_per_col_ + i] = (*col_g_sat)[i]; + //state.gas_density.data[column][i] = (*col_g_dens)[i]; + } + + if (has_ice) { + state.ice_density.data[column * ncells_per_col_ + i] = (*col_i_dens)[i]; + props.ice_saturation.data[column * ncells_per_col_ + i] = (*col_i_sat)[i]; + } + + } + //fill surface variables + + state.surface_energy_source.data[column] = surface_energy_source[column]; + state.surface_water_source.data[column] = surface_water_source[column]; + state.snow_depth.data[column] = snow_depth[0][column]; + state.canopy_longwave_radiation.data[column] = canopy_longwave_radiation[0][column]; + state.boundary_latent_heat_flux.data[column] = canopy_latent_heat[0][column]; + state.boundary_sensible_heat_flux.data[column] = canopy_sensible_heat[0][column]; + state.canopy_surface_water.data[column] = canopy_surface_water[0][column]; + state.transpiration.data[column] = transpiration[0][column]; + state.evaporation_canopy.data[column] = evaporation_canopy[0][column]; + state.evaporation_bare_ground.data[column] = evaporation_ground[0][column]; + state.evaporation_litter.data[column] = evaporation_litter[0][column]; + state.evaporation_snow.data[column] = evaporation_snow[0][column]; + state.sublimation_snow.data[column] = sublimation_snow[0][column]; + + props.shortwave_radiation.data[column] = shortwave_radiation[column]; + //props.longwave_radiation.data[column] = longwave_radiation[column]; + props.air_temperature.data[column] = air_temperature[column]; + props.vapor_pressure_air.data[column] = vapor_pressure_air[column]; + props.wind_speed.data[column] = wind_speed[column]; + props.elevation.data[column] = elevation[column]; + props.aspect.data[column] = aspect[column]; + props.slope.data[column] = slope[column]; + props.LAI.data[column] = LAI[column]; + props.SAI.data[column] = SAI[column]; + props.snow_albedo.data[column] = snow_albedo[column]; + props.vegetation_type.data[column] = vegetation_type[column]; + + if(p_bool){ + props.precipitation.data[column] = (*precipitation)[column]; + } else { + props.precipitation.data[column] = (*precipitation)[column]; + props.precipitation_snow.data[column] = (*precipitation_snow)[column]; + } + /*Don't need this loop until transport is implemented + for (int i = 0; i < state.total_component_concentration.columns; i++) { + for (int j = 0; j < state.total_component_concentration.cells; j++) { + for (int k = 0; k < state.total_component_concentration.components; k++) { + } + } + }*/ + } + + //Fill the atmospheric abundances + //NOTE: probably want to add an if statement here to only do this only once + props.atm_n2 = atm_n2_; + props.atm_o2 = atm_o2_; + props.atm_co2 = atm_co2_; + props.atm_ch4 = atm_ch4_; + props.atm_n2o = atm_n2o_; + props.atm_h2 = atm_h2_; + props.atm_nh3 = atm_nh3_; + props.heat_capacity = c_m_; + props.field_capacity = pressure_at_field_capacity; + props.wilting_point = pressure_at_wilting_point; + props.p_bool = p_bool; + props.a_bool = a_bool; + props.pheno_bool = pheno_bool; + + std::cout << "Data from state after setting struct: " << std::endl; + for (int col=0; col!=num_columns_local; ++col) { + if (std::isnan(surface_water_source[col]) || + std::isinf(surface_water_source[col])) { + std::cout << "Process " << p_rank << " found bad value at column " + << col << ": " << surface_water_source[col] << std::endl; + } + } + +} + +void EcoSIM::CopyFromEcoSIM_process(const int column, + const BGCProperties& props, + const BGCState& state, + const BGCAuxiliaryData& aux_data, + const Tag& water_tag) +{ + + //Transport removal + /*Epetra_MultiVector& mole_fraction= *(S_->GetPtrW(mole_fraction_key_, Tags::DEFAULT, "subsurface transport")->ViewComponent("cell",false)); + int mole_fraction_num = mole_fraction.NumVectors();*/ + int mole_fraction_num = 1; + + auto& porosity = *(*S_->GetW(porosity_key_, Tags::DEFAULT, porosity_key_).ViewComponent("cell",false))(0); + auto& liquid_saturation = *(*S_->GetW(saturation_liquid_key_, Tags::DEFAULT, saturation_liquid_key_).ViewComponent("cell",false))(0); + auto& water_content = *(*S_->GetW(water_content_key_, Tags::DEFAULT, water_content_key_).ViewComponent("cell",false))(0); + //auto& suction_head = *(*S_->GetW(suc_key_, Tags::DEFAULT, suc_key_).ViewComponent("cell",false))(0); + //auto& relative_permeability = *(*S_->GetW(relative_permeability_key_, Tags::DEFAULT, relative_permeability_key_).ViewComponent("cell",false))(0); + auto& liquid_density = *(*S_->GetW(liquid_density_key_, Tags::DEFAULT, liquid_density_key_).ViewComponent("cell",false))(0); + auto& rock_density = *(*S_->GetW(rock_density_key_, Tags::DEFAULT, rock_density_key_).ViewComponent("cell",false))(0); + auto& cell_volume = *(*S_->GetW(cell_volume_key_, Tags::DEFAULT, cell_volume_key_).ViewComponent("cell",false))(0); + + auto& surface_energy_source = *(*S_->GetW(surface_energy_source_ecosim_key_, Tags::DEFAULT, name_).ViewComponent("cell", false))(0); + auto& subsurface_energy_source = *(*S_->GetW(subsurface_energy_source_ecosim_key_, Tags::DEFAULT, subsurface_energy_source_ecosim_key_).ViewComponent("cell", false))(0); + + auto& surface_water_source = *(*S_->GetW(surface_water_source_ecosim_key_, Tags::DEFAULT, surface_water_source_ecosim_key_).ViewComponent("cell", false))(0); + auto& subsurface_water_source = *(*S_->GetW(subsurface_water_source_ecosim_key_, Tags::DEFAULT, subsurface_water_source_ecosim_key_).ViewComponent("cell", false))(0); + auto& temp = *(*S_->GetW(T_key_, Tags::DEFAULT, "subsurface energy").ViewComponent("cell",false))(0); + auto& thermal_conductivity = *(*S_->GetW(thermal_conductivity_key_, Tags::DEFAULT, thermal_conductivity_key_).ViewComponent("cell",false))(0); + auto& snow_temperature = *(*S_->GetW(snow_temperature_key_, Tags::DEFAULT, snow_temperature_key_).ViewComponent("cell", false))(0); + + auto& snow_depth = *S_->GetW(snow_depth_key_,tag_next_,snow_depth_key_).ViewComponent("cell"); + auto& canopy_longwave_radiation = *S_->GetW(canopy_lw_key_, tag_next_, canopy_lw_key_).ViewComponent("cell"); + auto& canopy_latent_heat = *S_->GetW(canopy_latent_heat_key_, tag_next_, canopy_latent_heat_key_).ViewComponent("cell"); + auto& canopy_sensible_heat = *S_->GetW(canopy_sensible_heat_key_, tag_next_, canopy_sensible_heat_key_).ViewComponent("cell"); + auto& canopy_surface_water = *S_->GetW(canopy_surface_water_key_, tag_next_, canopy_surface_water_key_).ViewComponent("cell"); + auto& transpiration = *S_->GetW(transpiration_key_, tag_next_, transpiration_key_).ViewComponent("cell"); + auto& evaporation_canopy = *S_->GetW(evaporation_canopy_key_, tag_next_, evaporation_canopy_key_).ViewComponent("cell"); + auto& evaporation_ground = *S_->GetW(evaporation_ground_key_, tag_next_, evaporation_ground_key_).ViewComponent("cell"); + auto& evaporation_litter = *S_->GetW(evaporation_litter_key_, tag_next_, evaporation_litter_key_).ViewComponent("cell"); + auto& evaporation_snow = *S_->GetW(evaporation_snow_key_, tag_next_, evaporation_snow_key_).ViewComponent("cell"); + auto& sublimation_snow = *S_->GetW(sublimation_snow_key_, tag_next_, sublimation_snow_key_).ViewComponent("cell"); + + auto col_porosity = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_l_sat = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_wc = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_suc = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_relative_permeability = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_l_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_r_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_vol = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_g_sat = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_g_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_i_sat = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_i_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_temp = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_cond = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_h_cond = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_b_dens = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_ss_energy_source = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_ss_water_source = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + auto col_snow_temperature = Teuchos::rcp(new Epetra_SerialDenseVector(ncells_per_col_)); + + auto col_mole_fraction = Teuchos::rcp(new Epetra_SerialDenseMatrix(mole_fraction_num,ncells_per_col_)); + + //Gather columns on this process: + num_columns_global = mesh_surf_->getMap(AmanziMesh::Entity_kind::CELL, false).NumGlobalElements(); + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + num_columns_global_ptype = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::ALL); + + //Trying to loop over processors now: + int p_rank; + MPI_Comm_rank(MPI_COMM_WORLD, &p_rank); + MPI_Barrier(MPI_COMM_WORLD); + + std::cout << "Data from struct after pass back: " << std::endl; + for (int col=0; col!=num_columns_local; ++col) { + if (std::isnan(state.surface_water_source.data[col]) || + std::isinf(state.surface_water_source.data[col])) { + std::cout << "Process " << p_rank << " found bad value at column " + << col << ": " << state.surface_water_source.data[col] << std::endl; + } + } + + num_columns_local = mesh_surf_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + double energy_source_tot = state.surface_energy_source.data[0]; + double water_source_tot = state.surface_water_source.data[0]; + double snow_depth_cell = state.snow_depth.data[0]; + + for (int col=0; col!=num_columns_local; ++col) { + for (int i=0; i < ncells_per_col_; ++i) { + (*col_ss_water_source)[i] = state.subsurface_water_source.data[col * ncells_per_col_ + i]; + (*col_ss_energy_source)[i] = state.subsurface_energy_source.data[col * ncells_per_col_ + i]; + (*col_snow_temperature)[i] = state.snow_temperature.data[col * ncells_per_col_ + i]; + } + + ColumnToField_(col, subsurface_water_source, col_ss_water_source.ptr()); + ColumnToField_(col, subsurface_energy_source, col_ss_energy_source.ptr()); + ColumnToField_(col, snow_temperature, col_snow_temperature.ptr()); + + surface_energy_source[col] = state.surface_energy_source.data[col]/(3600.0); + surface_water_source[col] = state.surface_water_source.data[col]/(3600.0); + snow_depth[0][col] = state.snow_depth.data[col]; + canopy_longwave_radiation[0][col] = state.canopy_longwave_radiation.data[col]; + canopy_latent_heat[0][col] = state.boundary_latent_heat_flux.data[col]; + canopy_sensible_heat[0][col] = state.boundary_sensible_heat_flux.data[col]; + canopy_surface_water[0][col] = state.canopy_surface_water.data[col]; + transpiration[0][col] = state.transpiration.data[col]; + evaporation_canopy[0][col] = state.evaporation_canopy.data[col]; + evaporation_ground[0][col] = state.evaporation_bare_ground.data[col]; + evaporation_litter[0][col] = state.evaporation_litter.data[col]; + evaporation_snow[0][col] = state.evaporation_snow.data[col]; + sublimation_snow[0][col] = state.sublimation_snow.data[col]; + } + + std::cout << "(CopyFromEcoSIM) subsurface energy flux: " << std::endl; + + /*for (int col=0; col!=num_columns_local; ++col) { + for (int i=0; i < ncells_per_col_; ++i) { + std::cout << "col: " << col << " cell: " << i << "value: " << subsurface_energy_source[col*ncells_per_col_+i] << std::endl; + } + } + + for (int col=0; col!=num_columns_local; ++col) { + for (int i=0; i < ncells_per_col_; ++i) { + std::cout << "col: " << col << " cell: " << i << "value: " << subsurface_water_source[col*ncells_per_col_+i] << std::endl; + } + }*/ +} + +int EcoSIM::InitializeSingleProcess(int proc) +{ + int num_iterations = 1; + int num_columns = 1; + + num_columns = num_columns_local; + + Teuchos::OSTab tab = vo_->getOSTab(); + + CopyToEcoSIM_process(proc, bgc_props_, bgc_state_, bgc_aux_data_, Tags::DEFAULT); + + bgc_sizes_.num_columns = num_columns; + bgc_sizes_.ncells_per_col_ = ncells_per_col_; + bgc_sizes_.num_components = 1; + + bgc_engine_->Setup(bgc_props_, bgc_state_, bgc_sizes_, num_iterations, num_columns,ncells_per_col_); + CopyFromEcoSIM_process(proc, bgc_props_, bgc_state_, bgc_aux_data_, Tags::DEFAULT); +} + +int EcoSIM::AdvanceSingleProcess(double dt, int proc) +{ + //Function to run EcoSIMs advance function, ecosim is run as an hourly model + // Every hour the data is taken from ATS state copied to ecosim + // The model is run and passed back to ATS. + + int num_iterations = 1; + int num_columns = 1; + + num_columns = num_columns_local; + + // Time tracking variables + current_time_ = S_->get_time(); //Current time + static double last_ecosim_time = 0.0; + int total_days = static_cast(current_time_ / 86400.0); + int current_day = (day0_ + total_days) % 365; + int current_year = year0_ + ((day0_ + total_days)/365); + + bgc_props_.current_day = current_day; + bgc_props_.current_year = current_year; + + Teuchos::OSTab tab = vo_->getOSTab(); + + if (current_time_ - last_ecosim_time >= 3600.0) { + *vo_->os() << "Hour completed at total_time: " << current_time_ + << ", Year: " << current_year << ", Day: " << current_day << std::endl; + *vo_->os() << "Running EcoSIM Advance: " << std::endl; + + CopyToEcoSIM_process(proc, bgc_props_, bgc_state_, bgc_aux_data_, Tags::DEFAULT); + + bgc_engine_->Advance(dt, bgc_props_, bgc_state_, bgc_sizes_, num_iterations, num_columns); + + CopyFromEcoSIM_process(proc, bgc_props_, bgc_state_, bgc_aux_data_, Tags::DEFAULT); + + last_ecosim_time = current_time_; + } + + return num_iterations; +} + +double** ConvertTo2DArray(BGCMatrixDouble* matrix) { + double** data_2d = new double*[matrix->cells]; + for (int i = 0; i < matrix->cells; ++i) { + data_2d[i] = &(matrix->data[i * matrix->capacity_columns]); + } + return data_2d; +} + + +} // namespace EcoSIM +} // namespace Amanzi diff --git a/src/pks/ecosim/EcoSIM_ATS_interface.hh b/src/pks/ecosim/EcoSIM_ATS_interface.hh new file mode 100644 index 000000000..b7c49bb02 --- /dev/null +++ b/src/pks/ecosim/EcoSIM_ATS_interface.hh @@ -0,0 +1,363 @@ +/*-------------------------------------------------------------------------- + ATS + + License: see $ATS_DIR/COPYRIGHT + Author: Andrew Graus (agraus@lbl.gov) + + --------------------------------------------------------------------------*/ +/*! + +This PK couples ATS to the BGC code EcoSIM. This PK essentially takes over the +surface balance aspects of ATS and replaces them with those from EcoSIM. + +In addition this code takes data from ATS state and loads it into a struct that +is fortran readable (BGCContainer). The data structures and methods were adapted +from those used in the Alquimia code, additionally the general code structure +Engine code are adapted from Alquimia as well. + +Structures for looping over cells of columns were adapted from ATS's simpleBGC code + +`"PK type`" = `"EcoSIM for ATS`" + +.. _pk-ecosim-spec: +.. admonition:: pk-ecosim-spec + + * `"engine`" ``[string]`` **EcoSIM** Engine name - inspired by Alquimias options if this + code is adapted to used to drive other BGC codes. + + * `"heat capacity [MJ mol^-1 K^-1]`" ``[double]]`` **0.02** heat capacity of the soil layers + + * `"Field Capacity [MPa]`" ``[double]]`` **-0.033** pressure at field capacity + + * `"Wilting Point [MPa]`" ``[double]]`` **-1.5** pressure at wilting point + + * `"initial time step [s]`" ``[double]]`` **3600.0** EcoSIM is an hourly model so the + standard is to run it after the end of hour 1 + + * `"EcoSIM Precipitation`" ``[bool]`` This allows EcoSIM to partition the precipitation + itself. If false it will expect the precipitation forcing to be already divided into + rain and snow as in ATS. + + * `"Prescribe Albedo`" ``[bool]`` determines if the code will use EcoSIM's internal + albedo calculation, or if it will be prescirbed from data. + + * `"Prescribe Phenology`" ``[bool]`` Determines if the coupling will use the prescribed + phenology methodology where LAI and PFT are input. This will eventually be complemented + with a full phenology method, where plant parameters are directly, but will remain in + the code as an option. + + * `"starting day of year [0-364]`" ``[int]`` day of the year, needed for EcoSIMs internal + radiation and biogeochemical processes. + + * `"Starting year`" ``[int]`` Year also needed for internal EcoSIM computations + + * `"Number of PFTs [1-5] "`" ``[int]`` Number of PFTs allowed in every column. 5 + is the maximum number allowed + + * `"domain name`" ``[string]`` **domain** + + * `"surface domain name`" ``[string]`` **surface** + + EVALUATORS: + + - `"Bulk Density`" `[ ]` + - `"Hydraulic Conductivity `[Pa]` + - `"Matric Pressure`" `[Pa]` + + DEPENDENCIES + //Sources + `"surface water source ecosim`" **surface-ecosim_water_source** + `"surface energy source ecosim`" **surface-ecosim_source** + `"subsurface water source ecosim`" **ecosim_water_source** + `"surface water source ecosim`" **surface-ecosim_water_source** + + //surface balance variables + `"incoming shortwave radiation`" **surface-incoming_shortwave_radiation** + `"incoming longwave radiation`" **surface-incoming_longwave_radiation** + `"air temperature`" **surface-air_temperature** + `"vapor pressure air`" **surface-vapor_pressure_air** + `"wind speed`" **surface-wind_speed** + `"precipitation rain`" **surface-precipitation_rain** + `"precipitation snow`" **surface-precipitation_snow** + `"precipitation total`" **surface-precipitation_total** + `"snow depth`" **surface-snow_depth** + `"snow_albedo`" **surface-snow_albedo** + `"snow temperature`" **surface-snow_temperature** + `"canopy longwave radiation`" **surface-canopy_longwave_radiation** + `"canopy latent heat`" **surface-canopy_latent_heat** + `"canopy sensible heat`" **surface-canopy_sensible_heat** + `"canopy surface water`" **surface-canopy_surface_water** + `"evapotranspiration`" **surface-evapotranspiration** + `"evaporation ground`" **surface-evaporation_ground** + `"evaporation litter`" **surface-evaporation_litter** + `"evaporation snow`" **surface-evaporation_snow** + `"sublimation snow`" **surface-sublimation_snow** + `"LAI`" **surface-LAI** + `"SAI`" **surface-SAI** + `"vegetation type`" **surface-vegetation_type** + + //General Flow Transport Energy + `"mole fraction`" **mole_fraction** + `"porosity`" **porosity** + `"saturation liquid`" **saturation_liquid** + `"saturation gas`" **saturation_gas** + `"saturation ice`" **saturation_ice** + `"water content`" **water_content** + `"mass density liquid`" **mass_density_liquid** + `"mass density ice`" **mass_density_ice** + `"mass density gas`" **mass_density_gas** + `"density rock`" **density_rock** + `"temperature`" **temperature** + `"thermal conductivity`" **thermal_conductivity** + `"cell volume`" **cell_volume** + + */ + + +#ifndef PKS_ECOSIM_HH_ +#define PKS_ECOSIM_HH_ + +#include +#include +#include + +#include "Epetra_MultiVector.h" +#include "Teuchos_ParameterList.hpp" +#include "Teuchos_RCP.hpp" +#include "Epetra_SerialDenseVector.h" +#include "Epetra_SerialDenseMatrix.h" + +#include "VerboseObject.hh" +#include "TreeVector.hh" + +#include "Key.hh" +#include "Mesh.hh" +#include "State.hh" +#include "BGCEngine.hh" +#include "PK_Factory.hh" +#include "PK_Physical_Default.hh" +#include "PK_Physical.hh" +#include "MeshPartition.hh" + +namespace Amanzi { +namespace EcoSIM { + +//using namespace Amanzi::Flow; + +class EcoSIM : public PK_Physical_Default { + + public: + + //Unclear if the constructor is neccessary + EcoSIM(Teuchos::ParameterList& pk_tree, + const Teuchos::RCP& plist, + const Teuchos::RCP& S, + const Teuchos::RCP& solution); + // Virtual destructor + ~EcoSIM(); + + // is a PK + // -- Setup data + //virtual void Setup(const Teuchos::Ptr&S); + virtual void Setup() final; + + // -- initalize owned (dependent) variables + //virtual void Initialize(const Teuchos::Ptr& S); + virtual void Initialize() final; + + // --provide timestep size + virtual double get_dt() final { + return dt_; + } + + virtual void set_dt(double dt) final { + dt_ = dt; + } + + // -- commit the model + //virtual void CommitStep(double t_old, double t_new, const Teuchos::RCP& S); + virtual void CommitStep(double t_old, double t_new, const Tag& tag) final; + + // -- Update diagnostics for vis. + //virtual void CalculateDiagnostics(const Teuchos::RCP& S) {} + //virtual void CalculateDiagnostics(const Tag& tag) override; + + // -- advance the model + virtual bool AdvanceStep(double t_old, double t_new, bool reinit) final; + + virtual std::string name(){return "EcoSIM for ATS";}; + + //This is not in the Alquimia_PK, for whatever reason it is defined in + //The Chemistry_PK even though it isn't used there, and then included + Teuchos::RCP bgc_engine() { return bgc_engine_; } + + private: + + //Helper functions from Alquimia + void CopyToEcoSIM(int column, + BGCProperties& props, + BGCState& state, + BGCAuxiliaryData& aux_data, + const Tag& water_tag = Tags::DEFAULT); + + void CopyFromEcoSIM(const int cell, + const BGCProperties& props, + const BGCState& state, + const BGCAuxiliaryData& aux_data, + const Tag& water_tag = Tags::DEFAULT); + + //Helper functions from Alquimia + void CopyToEcoSIM_process(int proc, + BGCProperties& props, + BGCState& state, + BGCAuxiliaryData& aux_data, + const Tag& water_tag = Tags::DEFAULT); + + void CopyFromEcoSIM_process(const int proc, + const BGCProperties& props, + const BGCState& state, + const BGCAuxiliaryData& aux_data, + const Tag& water_tag = Tags::DEFAULT); + + int InitializeSingleProcess(int proc); + + int AdvanceSingleProcess(double dt, int proc); + + void ComputeNextTimeStep(); + + protected: + double dt_; + double c_m_; + Teuchos::RCP mesh_surf_; //might need this? + Key domain_surface_; + std::string passwd_ = "state"; + + //The helper functions from BGC are protected not private (unclear why) + //I don't think I need this here, probably in the engine + void FieldToColumn_(AmanziMesh::Entity_ID column, const Epetra_Vector& vec, + Teuchos::Ptr col_vec); + + void FieldToColumn_(AmanziMesh::Entity_ID column, Teuchos::Ptr vec, + Teuchos::Ptr col_vec); + + void MatrixFieldToColumn_(AmanziMesh::Entity_ID column, const Epetra_MultiVector& m_arr, + Teuchos::Ptr col_arr); + + //void FieldToColumn_(AmanziMesh::Entity_ID column, const Epetra_Vector& vec, + // double* col_vec); + void ColDepthDz_(AmanziMesh::Entity_ID column, + Teuchos::Ptr depth, + Teuchos::Ptr dz); + + void VolDepthDz_(AmanziMesh::Entity_ID column, + Teuchos::Ptr depth, + Teuchos::Ptr dz, + Teuchos::Ptr volume); + + void ColumnToField_(AmanziMesh::Entity_ID column, Epetra_Vector& vec, + Teuchos::Ptr col_vec); + + void ColumnToField_(AmanziMesh::Entity_ID column, Teuchos::Ptr vec, + Teuchos::Ptr col_vec); + + void MatrixColumnToField_(AmanziMesh::Entity_ID column, Epetra_MultiVector& m_arr, + Teuchos::Ptr col_arr); + + int number_aqueous_components_; + int ncells_per_col_; + int num_columns_; + int num_columns_local; + int num_columns_global; + int num_columns_global_ptype; + int day0_, year0_, curr_day_, curr_year_; + double saved_time_; + double current_time_; + double t_ecosim = 0.0; + + // keys + Key mole_fraction_key_; + Key porosity_key_; + Key saturation_liquid_key_; + Key saturation_gas_key_; + Key saturation_ice_key_; + Key elev_key_; + Key water_content_key_; + Key relative_permeability_key_; + Key liquid_density_key_; + Key ice_density_key_; + Key gas_density_key_; + Key gas_density_key_test_; + Key rock_density_key_; + Key T_key_; + Key thermal_conductivity_key_; + Key cell_volume_key_; + Key ecosim_aux_data_key_; + Key bulk_density_key_; + Key hydraulic_conductivity_key_; + Key sw_key_; + Key lw_key_; + Key air_temp_key_; + Key vp_air_key_; + Key wind_speed_key_; + Key p_rain_key_; + Key p_snow_key_; + Key p_total_key_; + Key f_wp_key_; + Key f_root_key_; + Key matric_pressure_key_; + Key aspect_key_; + Key slope_key_; + Key lai_key_; + Key sai_key_; + Key v_type_key_; + Key surface_energy_source_key_; + Key subsurface_energy_source_key_; + Key surface_water_source_key_; + Key subsurface_water_source_key_; + Key surface_energy_source_ecosim_key_; + Key surface_water_source_ecosim_key_; + Key subsurface_energy_source_ecosim_key_; + Key subsurface_water_source_ecosim_key_; + Key snow_depth_key_; + Key snow_albedo_key_; + Key canopy_lw_key_; + Key canopy_latent_heat_key_; + Key canopy_sensible_heat_key_; + Key canopy_surface_water_key_; + Key transpiration_key_; + Key evaporation_canopy_key_; + Key evaporation_ground_key_; + Key evaporation_litter_key_; + Key evaporation_snow_key_; + Key sublimation_snow_key_; + Key snow_temperature_key_; + Key cap_pres_key_; + + Teuchos::RCP bgc_engine_; + + double atm_n2_, atm_o2_, atm_co2_, atm_ch4_, atm_n2o_, atm_h2_, atm_nh3_; + double pressure_at_field_capacity, pressure_at_wilting_point; + + private: + BGCState bgc_state_; + BGCProperties bgc_props_; + BGCAuxiliaryData bgc_aux_data_; + BGCSizes bgc_sizes_; + + Teuchos::RCP column_vol_save; + Teuchos::RCP column_wc_save; + + bool bgc_initialized_; + bool has_energy, has_gas, has_ice, p_bool, a_bool, pheno_bool; + std::vector component_names_; + int num_components; + + private: + //factory registration + static RegisteredPKFactory reg_; +}; + +} // namespace EcoSIM +} // namespace Amanzi + +#endif diff --git a/src/pks/ecosim/EcoSIM_ATS_interface_reg.hh b/src/pks/ecosim/EcoSIM_ATS_interface_reg.hh new file mode 100644 index 000000000..8179c2c00 --- /dev/null +++ b/src/pks/ecosim/EcoSIM_ATS_interface_reg.hh @@ -0,0 +1,18 @@ +/* ------------------------------------------------------------------------- + * ATS + * + * License: see $ATS_DIR/COPYRIGHT + * Author: Ethan Coon + * + * ------------------------------------------------------------------------- */ + +#include "EcoSIM_ATS_interface.hh" + +namespace Amanzi { +namespace EcoSIM { + +RegisteredPKFactory EcoSIM::reg_("EcoSIM for ATS"); + + +} // namespace +} // namespace diff --git a/src/pks/ecosim/constitutive_relations/CMakeLists.txt b/src/pks/ecosim/constitutive_relations/CMakeLists.txt new file mode 100644 index 000000000..62a7ba000 --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/CMakeLists.txt @@ -0,0 +1,73 @@ +# -*- mode: cmake -*- +# +# ATS +# Constitutive relations for EcoSIM +# + +#add_subdirectory(bulk_density) +#add_subdirectory(hydraulic_conductivity) + +# collect all sources +set(ats_ecosim_relations_src_files +# bulk_density/bulk_density_evaluator.cc +# bulk_density/bulk_density_model.cc + hydraulic_conductivity/hydraulic_conductivity_evaluator.cc + hydraulic_conductivity/hydraulic_conductivity_model.cc +# matric_pressure/matric_pressure_evaluator.cc +# matric_pressure/matric_pressure_model.cc +) + +set(ats_ecosim_relations_inc_files +# bulk_density/bulk_density_evaluator.hh +# bulk_density/bulk_density_model.hh + hydraulic_conductivity/hydraulic_conductivity_evaluator.hh + hydraulic_conductivity/hydraulic_conductivity_model.hh +# matric_pressure/matric_pressure_evaluator.hh +# matric_pressure/matric_pressure_model.hh +) + +set(ats_ecosim_relations_link_libs + ${Teuchos_LIBRARIES} + ${Epetra_LIBRARIES} + error_handling + atk + mesh + data_structures + whetstone + operators + solvers + time_integration + state + pks + chemistry_pk + ats_pks + ats_eos + ats_operators +) + +# make the library +add_amanzi_library(ats_ecosim_relations + SOURCE ${ats_ecosim_relations_src_files} + HEADERS ${ats_ecosim_relations_inc_files} + LINK_LIBS ${ats_ecosim_relations_link_libs}) + +#register_evaluator_with_factory( +# HEADERFILE bulk_density/bulk_density_evaluator_reg.hh +# LISTNAME ATS_ECOSIM_RELATIONS_REG +#) + +register_evaluator_with_factory( + HEADERFILE hydraulic_conductivity/hydraulic_conductivity_evaluator_reg.hh + LISTNAME ATS_ECOSIM_RELATIONS_REG +) + +#register_evaluator_with_factory( +# HEADERFILE matric_pressure/matric_pressure_evaluator_reg.hh +# LISTNAME ATS_ECOSIM_RELATIONS_REG +#) + +generate_evaluators_registration_header( + HEADERFILE ats_ecosim_relations_registration.hh + LISTNAME ATS_ECOSIM_RELATIONS_REG + INSTALL True + ) diff --git a/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity.py b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity.py new file mode 100644 index 000000000..64a913725 --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity.py @@ -0,0 +1,20 @@ +"""Richards water content evaluator: the standard form as a function of liquid saturation.""" + +import sys, os +sys.path.append(os.path.join(os.environ['ATS_SRC_DIR'], "tools", "evaluator_generator")) +from evaluator_generator import generate_evaluator + +deps = [("permeability","k"), + ("mass_density_liquid", "rho"), + ("viscosity_liquid", "mu"), + ] +params = [("g", "double", "gravitational constant")] + +import sympy +k, rho, mu = sympy.var("k,rho,mu") +g_ = sympy.var("g_") +expression = (k * rho * mu) / mu; + +generate_evaluator("hydraulic_conductivity", "ecosim", + "hydraulic conductivity", "hydraulic_conductivity", + deps, params, expression=expression, doc=__doc__) diff --git a/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator.cc b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator.cc new file mode 100644 index 000000000..986e0fd39 --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator.cc @@ -0,0 +1,156 @@ +/* + The hydraulic conductivity evaluator is an algebraic evaluator of a given model. +Richards water content evaluator: the standard form as a function of liquid saturation. + Generated via evaluator_generator. +*/ + +#include "hydraulic_conductivity_evaluator.hh" +#include "hydraulic_conductivity_model.hh" + +namespace Amanzi { +namespace Ecosim { +namespace Relations { + +// Constructor from ParameterList +HydraulicConductivityEvaluator::HydraulicConductivityEvaluator(Teuchos::ParameterList& plist) : + EvaluatorSecondaryMonotypeCV(plist) +{ + Teuchos::ParameterList& sublist = plist_.sublist("hydraulic_conductivity parameters"); + model_ = Teuchos::rcp(new HydraulicConductivityModel(sublist)); + InitializeFromPlist_(); +} + + +// Copy constructor +//Don't seem to need this +/*HydraulicConductivityEvaluator::HydraulicConductivityEvaluator(const HydraulicConductivityEvaluator& other) : + EvaluatorSecondaryMonotypeCV(other), + k_key_(other.k_key_), + rho_key_(other.rho_key_), + mu_key_(other.mu_key_), + model_(other.model_) {}*/ + + +// Virtual copy constructor +Teuchos::RCP +HydraulicConductivityEvaluator::Clone() const +{ + return Teuchos::rcp(new HydraulicConductivityEvaluator(*this)); +} + + +// Initialize by setting up dependencies +void +HydraulicConductivityEvaluator::InitializeFromPlist_() +{ + // Set up my dependencies + // - defaults to prefixed via domain + //Key domain_name = Keys::getDomain(my_key_); + Key domain_name = Keys::getDomain(my_keys_.front().first); + Tag tag = my_keys_.front().second; + + // - pull Keys from plist + // dependency: permeability + k_key_ = Keys::readKey(plist_, domain_name, "permeability", "permeability"); + dependencies_.insert(KeyTag{ k_key_, tag }); + + // dependency: mass_density_liquid + rho_key_ = Keys::readKey(plist_, domain_name, "mass density liquid", "mass_density_liquid"); + dependencies_.insert(KeyTag{ rho_key_, tag}); + + // dependency: viscosity_liquid + mu_key_ = Keys::readKey(plist_, domain_name, "viscosity liquid", "viscosity_liquid"); + dependencies_.insert(KeyTag{ mu_key_, tag}); +} + + +void +HydraulicConductivityEvaluator::Evaluate_(const State& S, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + Teuchos::RCP k = S.GetPtr(k_key_, tag); + Teuchos::RCP rho = S.GetPtr(rho_key_, tag); + Teuchos::RCP mu = S.GetPtr(mu_key_, tag); + + const AmanziGeometry::Point& gravity = S.Get("gravity", Tags::DEFAULT); + double gz = -gravity[2]; + + for (CompositeVector::name_iterator comp=result[0]->begin(); + comp!=result[0]->end(); ++comp) { + const Epetra_MultiVector& k_v = *k->ViewComponent(*comp, false); + const Epetra_MultiVector& rho_v = *rho->ViewComponent(*comp, false); + const Epetra_MultiVector& mu_v = *mu->ViewComponent(*comp, false); + Epetra_MultiVector& result_v = *result[0]->ViewComponent(*comp,false); + + int ncomp = result[0]->size(*comp, false); + for (int i=0; i!=ncomp; ++i) { + result_v[0][i] = model_->HydraulicConductivity(k_v[0][i], rho_v[0][i], mu_v[0][i],gz); + } + } +} + + +void +HydraulicConductivityEvaluator::EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) +{ + Tag tag = my_keys_.front().second; + Teuchos::RCP k = S.GetPtr(k_key_, tag); + Teuchos::RCP rho = S.GetPtr(rho_key_, tag); + Teuchos::RCP mu = S.GetPtr(mu_key_, tag); + + const AmanziGeometry::Point& gravity = S.Get("gravity", Tags::DEFAULT); + double gz = -gravity[2]; + + if (wrt_key == k_key_) { + for (CompositeVector::name_iterator comp=result[0]->begin(); + comp!=result[0]->end(); ++comp) { + const Epetra_MultiVector& k_v = *k->ViewComponent(*comp, false); + const Epetra_MultiVector& rho_v = *rho->ViewComponent(*comp, false); + const Epetra_MultiVector& mu_v = *mu->ViewComponent(*comp, false); + Epetra_MultiVector& result_v = *result[0]->ViewComponent(*comp,false); + + int ncomp = result[0]->size(*comp, false); + for (int i=0; i!=ncomp; ++i) { + result_v[0][i] = model_->DHydraulicConductivityDPermeability(k_v[0][i], rho_v[0][i], mu_v[0][i],gz); + } + } + + } else if (wrt_key == rho_key_) { + for (CompositeVector::name_iterator comp=result[0]->begin(); + comp!=result[0]->end(); ++comp) { + const Epetra_MultiVector& k_v = *k->ViewComponent(*comp, false); + const Epetra_MultiVector& rho_v = *rho->ViewComponent(*comp, false); + const Epetra_MultiVector& mu_v = *mu->ViewComponent(*comp, false); + Epetra_MultiVector& result_v = *result[0]->ViewComponent(*comp,false); + + int ncomp = result[0]->size(*comp, false); + for (int i=0; i!=ncomp; ++i) { + result_v[0][i] = model_->DHydraulicConductivityDMassDensityLiquid(k_v[0][i], rho_v[0][i], mu_v[0][i],gz); + } + } + + } else if (wrt_key == mu_key_) { + for (CompositeVector::name_iterator comp=result[0]->begin(); + comp!=result[0]->end(); ++comp) { + const Epetra_MultiVector& k_v = *k->ViewComponent(*comp, false); + const Epetra_MultiVector& rho_v = *rho->ViewComponent(*comp, false); + const Epetra_MultiVector& mu_v = *mu->ViewComponent(*comp, false); + Epetra_MultiVector& result_v = *result[0]->ViewComponent(*comp,false); + + int ncomp = result[0]->size(*comp, false); + for (int i=0; i!=ncomp; ++i) { + result_v[0][i] = model_->DHydraulicConductivityDViscosityLiquid(k_v[0][i], rho_v[0][i], mu_v[0][i],gz); + } + } + + } else { + AMANZI_ASSERT(0); + } +} + + +} //namespace +} //namespace +} //namespace diff --git a/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator.hh b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator.hh new file mode 100644 index 000000000..c9255d435 --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator.hh @@ -0,0 +1,54 @@ +/* + The hydraulic conductivity evaluator is an algebraic evaluator of a given model. + + Generated via evaluator_generator with: +Richards water content evaluator: the standard form as a function of liquid saturation. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +#ifndef AMANZI_ECOSIM_HYDRAULIC_CONDUCTIVITY_EVALUATOR_HH_ +#define AMANZI_ECOSIM_HYDRAULIC_CONDUCTIVITY_EVALUATOR_HH_ + +#include "Factory.hh" +#include "EvaluatorSecondaryMonotype.hh" + +namespace Amanzi { +namespace Ecosim { +namespace Relations { + +class HydraulicConductivityModel; + +class HydraulicConductivityEvaluator : public EvaluatorSecondaryMonotypeCV { + public: + explicit HydraulicConductivityEvaluator(Teuchos::ParameterList& plist); + HydraulicConductivityEvaluator(const HydraulicConductivityEvaluator& other) = default; + virtual Teuchos::RCP Clone() const override; + + Teuchos::RCP get_model() { return model_; } + + protected: + // Required methods from EvaluatorSecondaryMonotypeCV + virtual void Evaluate_(const State& S, + const std::vector& result) override; + virtual void EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) override; + + void InitializeFromPlist_(); + + Key k_key_; + Key rho_key_; + Key mu_key_; + + Teuchos::RCP model_; + + private: + static Utils::RegisteredFactory reg_; + +}; + +} //namespace +} //namespace +} //namespace + +#endif diff --git a/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator_reg.hh b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator_reg.hh new file mode 100644 index 000000000..63570e4a7 --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_evaluator_reg.hh @@ -0,0 +1,11 @@ +#include "hydraulic_conductivity_evaluator.hh" + +namespace Amanzi { +namespace Ecosim { +namespace Relations { + +Utils::RegisteredFactory HydraulicConductivityEvaluator::reg_("hydraulic conductivity"); + +} //namespace +} //namespace +} //namespace diff --git a/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_model.cc b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_model.cc new file mode 100644 index 000000000..9696f36aa --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_model.cc @@ -0,0 +1,60 @@ +/* + The hydraulic conductivity model is an algebraic model with dependencies. + + Generated via evaluator_generator with: +Richards water content evaluator: the standard form as a function of liquid saturation. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +#include "Teuchos_ParameterList.hpp" +#include "dbc.hh" +#include "hydraulic_conductivity_model.hh" + +namespace Amanzi { +namespace Ecosim { +namespace Relations { + +// Constructor from ParameterList +HydraulicConductivityModel::HydraulicConductivityModel(Teuchos::ParameterList& plist) +{ + InitializeFromPlist_(plist); +} + + +// Initialize parameters +void +HydraulicConductivityModel::InitializeFromPlist_(Teuchos::ParameterList& plist) +{ + //g_ = plist.get("gravitational constant"); +} + + +// main method +double +HydraulicConductivityModel::HydraulicConductivity(double k, double rho, double mu, double gz) const +{ + return k*rho*gz/mu; +} + +double +HydraulicConductivityModel::DHydraulicConductivityDPermeability(double k, double rho, double mu, double gz) const +{ + return rho*gz/mu; +} + +double +HydraulicConductivityModel::DHydraulicConductivityDMassDensityLiquid(double k, double rho, double mu, double gz) const +{ + return k*gz/mu; +} + +double +HydraulicConductivityModel::DHydraulicConductivityDViscosityLiquid(double k, double rho, double mu, double gz) const +{ + return -1.0*k*rho*gz/pow(mu,2); +} + +} //namespace +} //namespace +} //namespace diff --git a/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_model.hh b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_model.hh new file mode 100644 index 000000000..5209ef609 --- /dev/null +++ b/src/pks/ecosim/constitutive_relations/hydraulic_conductivity/hydraulic_conductivity_model.hh @@ -0,0 +1,42 @@ +/* + The hydraulic conductivity model is an algebraic model with dependencies. + + Generated via evaluator_generator with: +Richards water content evaluator: the standard form as a function of liquid saturation. + + Authors: Ethan Coon (ecoon@lanl.gov) +*/ + +#ifndef AMANZI_ECOSIM_HYDRAULIC_CONDUCTIVITY_MODEL_HH_ +#define AMANZI_ECOSIM_HYDRAULIC_CONDUCTIVITY_MODEL_HH_ + +namespace Amanzi { +namespace Ecosim { +namespace Relations { + +class HydraulicConductivityModel { + + public: + explicit + HydraulicConductivityModel(Teuchos::ParameterList& plist); + + double HydraulicConductivity(double k, double rho, double mu, double gz) const; + + double DHydraulicConductivityDPermeability(double k, double rho, double mu, double gz) const; + double DHydraulicConductivityDMassDensityLiquid(double k, double rho, double mu, double gz) const; + double DHydraulicConductivityDViscosityLiquid(double k, double rho, double mu, double gz) const; + + protected: + void InitializeFromPlist_(Teuchos::ParameterList& plist); + + protected: + + double g_; + +}; + +} //namespace +} //namespace +} //namespace + +#endif diff --git a/src/pks/ecosim/data/BGC_containers.cc b/src/pks/ecosim/data/BGC_containers.cc new file mode 100644 index 000000000..ae5c04c3f --- /dev/null +++ b/src/pks/ecosim/data/BGC_containers.cc @@ -0,0 +1,50 @@ +/* -*- mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */ + +/* +** Alquimia Copyright (c) 2013-2016, The Regents of the University of California, +** through Lawrence Berkeley National Laboratory (subject to receipt of any +** required approvals from the U.S. Dept. of Energy). All rights reserved. +** +** Alquimia is available under a BSD license. See LICENSE.txt for more +** information. +** +** If you have questions about your rights to use or distribute this software, +** please contact Berkeley Lab's Technology Transfer and Intellectual Property +** Management at TTD@lbl.gov referring to Alquimia (LBNL Ref. 2013-119). +** +** NOTICE. This software was developed under funding from the U.S. Department +** of Energy. As such, the U.S. Government has been granted for itself and +** others acting on its behalf a paid-up, nonexclusive, irrevocable, worldwide +** license in the Software to reproduce, prepare derivative works, and perform +** publicly and display publicly. Beginning five (5) years after the date +** permission to assert copyright is obtained from the U.S. Department of Energy, +** and subject to any subsequent five (5) year renewals, the U.S. Government is +** granted for itself and others acting on its behalf a paid-up, nonexclusive, +** irrevocable, worldwide license in the Software to reproduce, prepare derivative +** works, distribute copies to the public, perform publicly and display publicly, +** and to permit others to do so. +** +** Authors: Benjamin Andre +*/ + +#include "BGC_containers.hh" + +#include "../ecosim_interface.h" + +//String lengths +const int kBGCMaxStringLength = 512; +const int kBGCMaxWordLength = 32; + +/* Its kind of silly to have this as a separate code at this point, but this +will in theory become the switching code if we add other bgc +codes, in alquimia this switches between CrunchFlow and PFloTran*/ + +void CreateBGCInterface(const char* const engine_name, BGCInterface* interface) + { + + interface->DataTest = &ecosim_datatest; + interface->Setup = &ecosim_setup; + interface->Shutdown = &ecosim_shutdown; + interface->Advance = &ecosim_advance; + + } /* end CreateBGCInterface() */ diff --git a/src/pks/ecosim/data/BGC_containers.hh b/src/pks/ecosim/data/BGC_containers.hh new file mode 100644 index 000000000..a2dc1b4e7 --- /dev/null +++ b/src/pks/ecosim/data/BGC_containers.hh @@ -0,0 +1,217 @@ +/* -*- mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */ + +/* +** Alquimia Copyright (c) 2013-2016, The Regents of the University of California, +** through Lawrence Berkeley National Laboratory (subject to receipt of any +** required approvals from the U.S. Dept. of Energy). All rights reserved. +** +** Alquimia is available under a BSD license. See LICENSE.txt for more +** information. +** +** If you have questions about your rights to use or distribute this software, +** please contact Berkeley Lab's Technology Transfer and Intellectual Property +** Management at TTD@lbl.gov referring to Alquimia (LBNL Ref. 2013-119). +** +** NOTICE. This software was developed under funding from the U.S. Department +** of Energy. As such, the U.S. Government has been granted for itself and +** others acting on its behalf a paid-up, nonexclusive, irrevocable, worldwide +** license in the Software to reproduce, prepare derivative works, and perform +** publicly and display publicly. Beginning five (5) years after the date +** permission to assert copyright is obtained from the U.S. Department of Energy, +** and subject to any subsequent five (5) year renewals, the U.S. Government is +** granted for itself and others acting on its behalf a paid-up, nonexclusive, +** irrevocable, worldwide license in the Software to reproduce, prepare derivative +** works, distribute copies to the public, perform publicly and display publicly, +** and to permit others to do so. +** +** Authors: Benjamin Andre +*/ + +#ifndef BGC_CONTAINERS_H_ +#define BGC_CONTAINERS_H_ + +/******************************************************************************* + ** + ** C implementation of the alquimia containers. + ** + ** These are passed directly into the fortran routines. The + ** signatures must match exactly with the fortran side of things. + ** + ******************************************************************************/ + +#include +#include + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +//Defining String Lengths +extern const int kBGCMaxStringLength; +extern const int kBGCMaxWordLength; + + typedef struct { + int size, capacity; + double* data; + } BGCVectorDouble; + + typedef struct { + int size, capacity; + int* data; + } BGCVectorInt; + + typedef struct { + /* NOTE: this is a vector of strings */ + int size, capacity; + char** data; + } BGCVectorString; + + typedef struct { + int cells, columns, capacity_cells, capacity_columns; + double* data; + } BGCMatrixDouble; + + typedef struct { + int cells, columns, capacity_cells, capacity_columns; + int* data; + } BGCMatrixInt; + + typedef struct { + int cells, columns, capacity; + char* data; + } BGCMatrixString; + + typedef struct { + int cells, columns, components, capacity_cells, capacity_columns, capacity_components; + double*** data; + } BGCTensorDouble; + + typedef struct { + int cells, columns, components, capacity_cells, capacity_columns, capacity_components; + int*** data; + } BGCTensorInt; + + typedef struct { + int cells, columns, components, capacity; + char*** data; + } BGCTensorString; + + typedef struct { + int ncells_per_col_; + int num_components; + int num_columns; + } BGCSizes; + + typedef struct { + BGCMatrixDouble liquid_density; + BGCMatrixDouble gas_density; + BGCMatrixDouble ice_density; + BGCMatrixDouble rock_density; + BGCMatrixDouble porosity; + BGCMatrixDouble water_content; + BGCMatrixDouble matric_pressure; + BGCMatrixDouble temperature; + BGCMatrixDouble hydraulic_conductivity; + BGCMatrixDouble bulk_density; + BGCMatrixDouble subsurface_water_source; + BGCMatrixDouble subsurface_energy_source; + BGCVectorDouble surface_energy_source; + BGCVectorDouble surface_water_source; + BGCVectorDouble snow_depth; + BGCVectorDouble canopy_longwave_radiation; + BGCVectorDouble boundary_latent_heat_flux; + BGCVectorDouble boundary_sensible_heat_flux; + BGCVectorDouble canopy_surface_water; + BGCVectorDouble transpiration; + BGCVectorDouble evaporation_canopy; + BGCVectorDouble evaporation_bare_ground; + BGCVectorDouble evaporation_litter; + BGCVectorDouble evaporation_snow; + BGCVectorDouble sublimation_snow; + BGCMatrixDouble snow_temperature; + BGCTensorDouble total_component_concentration; + } BGCState; + + typedef struct { + BGCMatrixDouble liquid_saturation; + BGCMatrixDouble gas_saturation; + BGCMatrixDouble ice_saturation; + BGCMatrixDouble relative_permeability; + BGCMatrixDouble thermal_conductivity; + BGCMatrixDouble volume; + BGCMatrixDouble depth; + BGCMatrixDouble depth_c; + BGCMatrixDouble dz; + BGCMatrixDouble plant_wilting_factor; + BGCMatrixDouble rooting_depth_fraction; + BGCVectorDouble column_area; + BGCVectorDouble shortwave_radiation; + BGCVectorDouble longwave_radiation; + BGCVectorDouble air_temperature; + BGCVectorDouble vapor_pressure_air; + BGCVectorDouble wind_speed; + BGCVectorDouble precipitation; + BGCVectorDouble precipitation_snow; + BGCVectorDouble elevation; + BGCVectorDouble aspect; + BGCVectorDouble slope; + BGCVectorDouble LAI; + BGCVectorDouble SAI; + BGCVectorDouble vegetation_type; + BGCVectorDouble snow_albedo; + double atm_n2; + double atm_o2; + double atm_co2; + double atm_ch4; + double atm_n2o; + double atm_h2; + double atm_nh3; + double heat_capacity; + double field_capacity; + double wilting_point; + int current_day; + int current_year; + bool p_bool; + bool a_bool; + bool pheno_bool; + } BGCProperties; + + typedef struct { + BGCVectorInt aux_ints; /* [-] */ + BGCVectorDouble aux_doubles; /* [-] */ + } BGCAuxiliaryData; + + typedef struct { + /* read data files/structures, initialize memory, basis management + (includes reading database, swapping basis, etc.) */ + void (*DataTest)(); + + void (*Setup)( + BGCProperties* properties, + BGCState* state, + BGCSizes* sizes, + int num_iterations, + int num_columns, + int ncells_per_col_); + + /* gracefully shutdown the engine, cleanup memory */ + void (*Shutdown)(); + + /* take one (or more?) reaction steps in operator split mode */ + void (*Advance)( + double delta_t, + BGCProperties* properties, + BGCState* state, + BGCSizes* sizes, + int num_iterations, + int num_columns); + + } BGCInterface; + + void CreateBGCInterface(const char* const engine_name, BGCInterface* interface); + +#ifdef __cplusplus +} +#endif /* __cplusplus */ + +#endif /* ALQUIMIA_CONTAINERS_H_ */ diff --git a/src/pks/ecosim/data/BGC_memory.cc b/src/pks/ecosim/data/BGC_memory.cc new file mode 100644 index 000000000..ca785cbd0 --- /dev/null +++ b/src/pks/ecosim/data/BGC_memory.cc @@ -0,0 +1,834 @@ +/* -*- mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */ + +/* +** Alquimia Copyright (c) 2013-2016, The Regents of the University of California, +** through Lawrence Berkeley National Laboratory (subject to receipt of any +** required approvals from the U.S. Dept. of Energy). All rights reserved. +** +** Alquimia is available under a BSD license. See LICENSE.txt for more +** information. +** +** If you have questions about your rights to use or distribute this software, +** please contact Berkeley Lab's Technology Transfer and Intellectual Property +** Management at TTD@lbl.gov referring to Alquimia (LBNL Ref. 2013-119). +** +** NOTICE. This software was developed under funding from the U.S. Department +** of Energy. As such, the U.S. Government has been granted for itself and +** others acting on its behalf a paid-up, nonexclusive, irrevocable, worldwide +** license in the Software to reproduce, prepare derivative works, and perform +** publicly and display publicly. Beginning five (5) years after the date +** permission to assert copyright is obtained from the U.S. Department of Energy, +** and subject to any subsequent five (5) year renewals, the U.S. Government is +** granted for itself and others acting on its behalf a paid-up, nonexclusive, +** irrevocable, worldwide license in the Software to reproduce, prepare derivative +** works, distribute copies to the public, perform publicly and display publicly, +** and to permit others to do so. +** +** Authors: Benjamin Andre +*/ + + +/******************************************************************************* + ** + ** Alquimia C memory utilities to handle memory management + ** + ** Notes: + ** + ** - calloc/malloc always return NULL pointers if they fail, so + ** there is no need to pre-assign NULL for the pointers we are + ** allocating here. For zero size or zero members, the returned + ** pointer should be NULL or something that can be freed.... + ** + ** - free just releases the memory, it does not change the value + ** of the pointer. After free, the pointer is no longer valid, so + ** we set it to NULL. + ** + *******************************************************************************/ + +#include +#include "BGC_memory.hh" +#include "BGC_containers.hh" + +// Returns the nearest power of 2 greater than or equal to n, or 0 if n == 0. +static inline int nearest_power_of_2(int n) +{ + if (n == 0) return 0; + int twop = 1; + while (twop < n) + twop *= 2; + return twop; +} + +/******************************************************************************* + ** + ** BGC Vectors + ** + *******************************************************************************/ +void AllocateBGCVectorDouble(const int size, BGCVectorDouble* vector) { + if (size > 0) { + vector->size = size; + vector->capacity = nearest_power_of_2(size); + vector->data = (double*) calloc((size_t)vector->capacity, sizeof(double)); + //ALQUIMIA_ASSERT(NULL != vector->data); + } else { + vector->size = 0; + vector->capacity = 0; + vector->data = NULL; + } +} /* end AllocateBGCVectorDouble() */ + +void FreeBGCVectorDouble(BGCVectorDouble* vector) { + if (vector != NULL) { + free(vector->data); + vector->data = NULL; + vector->size = 0; + vector->capacity = 0; + } +} /* end FreeBGCVectorDouble() */ + +void AllocateBGCVectorInt(const int size, BGCVectorInt* vector) { + if (size > 0) { + vector->size = size; + vector->capacity = nearest_power_of_2(size); + vector->data = (int*) calloc((size_t)vector->capacity, sizeof(int)); + //ALQUIMIA_ASSERT(NULL != vector->data); + } else { + vector->size = 0; + vector->capacity = 0; + vector->data = NULL; + } +} /* end AllocateBGCVectorInt() */ + +void FreeBGCVectorInt(BGCVectorInt* vector) { + if (vector != NULL) { + free(vector->data); + vector->data = NULL; + vector->size = 0; + vector->capacity = 0; + } +} /* end FreeBGCVectorInt() */ + +void AllocateBGCVectorString(const int size, BGCVectorString* vector) { + int i; + if (size > 0) { + vector->size = size; + vector->capacity = nearest_power_of_2(size); + vector->data = (char**) calloc((size_t)vector->capacity, sizeof(char*)); + //ALQUIMIA_ASSERT(NULL != vector->data); + for (i = 0; i < vector->size; ++i) { + vector->data[i] = (char*) calloc((size_t)kBGCMaxStringLength, sizeof(char)); + //ALQUIMIA_ASSERT(NULL != vector->data[i]); + } + } else { + vector->size = 0; + vector->capacity = 0; + vector->data = NULL; + } +} /* end AllocateBGCVectorString() */ + +void FreeBGCVectorString(BGCVectorString* vector) { + int i; + if (vector != NULL) { + for (i = 0; i < vector->size; ++i) { + free(vector->data[i]); + } + free(vector->data); + vector->data = NULL; + vector->size = 0; + vector->capacity = 0; + } +} /* end FreeBGCVectorString() */ + +/******************************************************************************* + ** + ** BGC Matrix + ** + *******************************************************************************/ +void AllocateBGCMatrixDouble(const int cells, const int columns, BGCMatrixDouble* matrix) { + if ((cells > 0 ) || (columns > 0)){ + matrix->cells = cells; + matrix->columns = columns; + matrix->capacity_cells = nearest_power_of_2(cells); + matrix->capacity_columns = nearest_power_of_2(columns); + matrix->data = (double*) calloc((size_t)matrix->capacity_cells * matrix->capacity_columns, sizeof(double)); + + //for (int i = 0; i < matrix->columns; ++i) { + // matrix->data[i] = (double*) calloc((size_t)matrix->capacity_cells, sizeof(double)); + //} + //ALQUIMIA_ASSERT(NULL != matrix->data); + } else { + matrix->cells= 0; + matrix->columns= 0; + matrix->capacity_cells= 0; + matrix->capacity_columns= 0; + matrix->data = NULL; + } +} /* end AllocateBGCmatrixDouble() */ + +void FreeBGCMatrixDouble(BGCMatrixDouble* matrix) { + if (matrix != NULL) { + free(matrix->data); + matrix->data = NULL; + matrix->cells= 0; + matrix->columns= 0; + matrix->capacity_cells= 0; + matrix->capacity_columns= 0; + } +} /* end FreeBGCmatrixDouble() */ + +void AllocateBGCMatrixInt(const int cells, const int columns, BGCMatrixInt* matrix) { + if ((cells> 0) || (columns> 0)) { + matrix->cells= cells; + matrix->columns= columns; + matrix->capacity_cells= nearest_power_of_2(cells); + matrix->capacity_columns= nearest_power_of_2(columns); + matrix->data = (int*) calloc((size_t)matrix->capacity_cells * matrix->capacity_columns, sizeof(int)); + //for (int i = 0; i < matrix->columns; ++i) { + // matrix->data[i] = (int*) calloc((size_t)matrix->capacity_cells, sizeof(int)); + //} + //ALQUIMIA_ASSERT(NULL != matrix->data); + } else { + matrix->cells= 0; + matrix->columns= 0; + matrix->capacity_cells= 0; + matrix->capacity_columns= 0; + matrix->data = NULL; + } +} /* end AllocateBGCMatrixInt() */ + +void FreeBGCMatrixInt(BGCMatrixInt* matrix) { + if (matrix != NULL) { + free(matrix->data); + matrix->data = NULL; + matrix->cells= 0; + matrix->columns= 0; + matrix->capacity_columns= 0; + matrix->capacity_cells= 0; + } +} /* end FreeBGCMatrixInt() */ + +/******************************************************************************* + ** + ** BGC Tensor + ** + *******************************************************************************/ + +void AllocateBGCTensorDouble(const int cells, const int columns, const int components, BGCTensorDouble* tensor) { + if ((cells> 0 ) || (columns> 0) || (components > 0)){ + tensor->cells = cells; + tensor->columns = columns; + tensor->components = components; + + tensor->capacity_cells= nearest_power_of_2(cells); + tensor->capacity_columns= nearest_power_of_2(columns); + tensor->capacity_components = nearest_power_of_2(components); + + tensor->data = (double***) calloc((size_t)tensor->capacity_columns, sizeof(double**)); + for (int i = 0; i < tensor->columns; ++i) { + tensor->data[i] = (double**) calloc((size_t)tensor->capacity_cells, sizeof(double*)); + for (int j = 0; j < tensor->cells; ++j) { + tensor->data[i][j] = (double*) calloc((size_t)tensor->capacity_components, sizeof(double)); + } + } + //ALQUIMIA_ASSERT(NULL != matrix->data); + } else { + tensor->cells= 0; + tensor->columns= 0; + tensor->components = 0; + tensor->capacity_cells = 0; + tensor->capacity_columns = 0; + tensor->capacity_components = 0; + tensor->data = NULL; + } +} /* end AllocateBGCmatrixDouble() */ + +void FreeBGCTensorDouble(BGCTensorDouble* tensor) { + if (tensor != NULL) { + free(tensor->data); + tensor->data = NULL; + tensor->cells= 0; + tensor->columns = 0; + tensor->capacity_cells = 0; + tensor->capacity_columns = 0; + tensor->capacity_components = 0; + } +} /* end FreeBGCmatrixDouble() */ + +void AllocateBGCTensorInt(const int cells, const int columns, const int components, BGCTensorInt* tensor) { + if ((cells> 0 ) || (columns> 0) || (components > 0)){ + tensor->cells= cells; + tensor->columns= columns; + tensor->components = components; + + tensor->capacity_cells= nearest_power_of_2(cells); + tensor->capacity_columns= nearest_power_of_2(columns); + tensor->capacity_components = nearest_power_of_2(components); + + tensor->data = (int***) calloc((size_t)tensor->capacity_columns, sizeof(int**)); + for (int i = 0; i < tensor->columns; ++i) { + tensor->data[i] = (int**) calloc((size_t)tensor->capacity_cells, sizeof(int*)); + for (int j = 0; j < tensor->cells; ++j) { + tensor->data[i][j] = (int*) calloc((size_t)tensor->capacity_components, sizeof(int)); + } + } + //ALQUIMIA_ASSERT(NULL != matrix->data); + } else { + tensor->cells= 0; + tensor->columns= 0; + tensor->components = 0; + tensor->capacity_cells= 0; + tensor->capacity_columns= 0; + tensor->capacity_components = 0; + tensor->data = NULL; + } +} /* end AllocateBGCmatrixint() */ + +void FreeBGCTensorInt(BGCTensorInt* tensor) { + if (tensor != NULL) { + free(tensor->data); + tensor->data = NULL; + tensor->cells= 0; + tensor->columns= 0; + tensor->capacity_cells= 0; + tensor->capacity_columns= 0; + tensor->capacity_components = 0; + } +} /* end FreeBGCmatrixint() */ + +/******************************************************************************* + ** + ** State + ** + *******************************************************************************/ +/*Note that sizes for all the datasets that are single vectors should just be +the size of the column, need to test +For reference the old function call was: +void AllocateBGCState(const BGCSizes* const sizes, + BGCState* state)*/ + + void AllocateBGCState(BGCSizes* sizes, BGCState* state, + int ncells_per_col_, int num_components, int num_columns) { + sizes->ncells_per_col_ = ncells_per_col_; + sizes->num_components = num_components; + sizes->num_columns = num_columns; + + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->liquid_density)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->gas_density)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->ice_density)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->rock_density)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->porosity)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->water_content)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->matric_pressure)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->temperature)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->hydraulic_conductivity)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->bulk_density)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->subsurface_energy_source)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(state->subsurface_water_source)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->surface_water_source)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->surface_energy_source)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->snow_depth)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->canopy_longwave_radiation)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->boundary_latent_heat_flux)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->boundary_sensible_heat_flux)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->canopy_surface_water)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->transpiration)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->evaporation_canopy)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->evaporation_bare_ground)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->evaporation_litter)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->evaporation_snow)); + AllocateBGCVectorDouble(sizes->num_columns, &(state->sublimation_snow)); + AllocateBGCMatrixDouble(sizes->num_columns, sizes->num_columns, &(state->snow_temperature)); + AllocateBGCTensorDouble(sizes->ncells_per_col_, sizes->num_columns, sizes->num_components, &(state->total_component_concentration)); + //ALQUIMIA_ASSERT(state->total_mobile.data != NULL); + } /* end AllocateBGCState() */ + + void FreeBGCState(BGCState* state) { + if (state != NULL) { + FreeBGCMatrixDouble(&(state->liquid_density)); + FreeBGCMatrixDouble(&(state->gas_density)); + FreeBGCMatrixDouble(&(state->ice_density)); + FreeBGCMatrixDouble(&(state->rock_density)); + FreeBGCMatrixDouble(&(state->porosity)); + FreeBGCMatrixDouble(&(state->water_content)); + FreeBGCMatrixDouble(&(state->matric_pressure)); + FreeBGCMatrixDouble(&(state->temperature)); + FreeBGCMatrixDouble(&(state->hydraulic_conductivity)); + FreeBGCMatrixDouble(&(state->bulk_density)); + FreeBGCMatrixDouble(&(state->subsurface_energy_source)); + FreeBGCMatrixDouble(&(state->subsurface_water_source)); + FreeBGCVectorDouble(&(state->surface_energy_source)); + FreeBGCVectorDouble(&(state->surface_water_source)); + FreeBGCVectorDouble(&(state->snow_depth)); + FreeBGCVectorDouble(&(state->canopy_longwave_radiation)); + FreeBGCVectorDouble(&(state->boundary_latent_heat_flux)); + FreeBGCVectorDouble(&(state->boundary_sensible_heat_flux)); + FreeBGCVectorDouble(&(state->canopy_surface_water)); + FreeBGCVectorDouble(&(state->transpiration)); + FreeBGCVectorDouble(&(state->evaporation_canopy)); + FreeBGCVectorDouble(&(state->evaporation_bare_ground)); + FreeBGCVectorDouble(&(state->evaporation_litter)); + FreeBGCVectorDouble(&(state->evaporation_snow)); + FreeBGCVectorDouble(&(state->sublimation_snow)); + FreeBGCMatrixDouble(&(state->snow_temperature)); + FreeBGCTensorDouble(&(state->total_component_concentration)); + } + } /* end FreeAlquimiaState() */ + + /******************************************************************************* + ** + ** Auxiliary Data + ** + *******************************************************************************/ + /* + void AllocateBGCAuxiliaryData(const BGCSizes* const sizes, BGCAuxiliaryData* aux_data, + int ncells_per_col_) { + AllocateBGCMatrixInt(sizes->ncells_per_col_, + &(aux_data->aux_ints)); + + AllocateBGCMatrixDouble(sizes->ncells_per_col_, + &(aux_data->aux_doubles)); + + } // end AllocateAlquimiaAuxiliaryData() + + void FreeBGCAuxiliaryData(BGCAuxiliaryData* aux_data) { + if (aux_data != NULL) { + FreeBGCMatrixInt(&(aux_data->aux_ints)); + FreeBGCMatrixDouble(&(aux_data->aux_doubles)); + } + } // end FreeAlquimiaAuxiliaryData() + */ + + /******************************************************************************* + ** + ** Properties + ** + *******************************************************************************/ + + void AllocateBGCProperties(BGCSizes* sizes, BGCProperties* properties, + int ncells_per_col_, int num_columns) { + + sizes->ncells_per_col_ = ncells_per_col_; + sizes->num_columns = num_columns; + //sizes->num_components = num_components; + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->liquid_saturation)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->gas_saturation)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->ice_saturation)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->relative_permeability)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->thermal_conductivity)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->volume)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->depth)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->depth_c)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->dz)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->plant_wilting_factor)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_, sizes->num_columns, &(properties->rooting_depth_fraction)); + + AllocateBGCVectorDouble(sizes->num_columns, &(properties->column_area)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->shortwave_radiation)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->longwave_radiation)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->air_temperature)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->vapor_pressure_air)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->wind_speed)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->precipitation)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->precipitation_snow)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->elevation)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->aspect)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->slope)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->LAI)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->SAI)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->vegetation_type)); + AllocateBGCVectorDouble(sizes->num_columns, &(properties->snow_albedo)); + } /* end AllocateAlquimiaProperties() */ + + void FreeBGCProperties(BGCProperties* properties) { + if (properties != NULL) { + FreeBGCMatrixDouble(&(properties->liquid_saturation)); + FreeBGCMatrixDouble(&(properties->gas_saturation)); + FreeBGCMatrixDouble(&(properties->ice_saturation)); + FreeBGCMatrixDouble(&(properties->relative_permeability)); + FreeBGCMatrixDouble(&(properties->thermal_conductivity)); + FreeBGCMatrixDouble(&(properties->volume)); + FreeBGCMatrixDouble(&(properties->depth)); + FreeBGCMatrixDouble(&(properties->dz)); + FreeBGCMatrixDouble(&(properties->plant_wilting_factor)); + FreeBGCMatrixDouble(&(properties->rooting_depth_fraction)); + + FreeBGCVectorDouble(&(properties->column_area)); + FreeBGCVectorDouble(&(properties->shortwave_radiation)); + FreeBGCVectorDouble(&(properties->longwave_radiation)); + FreeBGCVectorDouble(&(properties->air_temperature)); + FreeBGCVectorDouble(&(properties->vapor_pressure_air)); + FreeBGCVectorDouble(&(properties->wind_speed)); + FreeBGCVectorDouble(&(properties->precipitation)); + FreeBGCVectorDouble(&(properties->precipitation_snow)); + FreeBGCVectorDouble(&(properties->elevation)); + FreeBGCVectorDouble(&(properties->aspect)); + FreeBGCVectorDouble(&(properties->slope)); + FreeBGCVectorDouble(&(properties->LAI)); + FreeBGCVectorDouble(&(properties->SAI)); + FreeBGCVectorDouble(&(properties->vegetation_type)); + FreeBGCVectorDouble(&(properties->snow_albedo)); + } + } + +/* OLD VERSION OF THE DATA MODULES +void AllocateBGCState(BGCSizes* sizes, BGCState* state, + int ncells_per_col_, int num_components) { + sizes->ncells_per_col_ = ncells_per_col_; + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->liquid_density)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->gas_density)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->ice_density)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->porosity)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->water_content)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->suction_head)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->temperature)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->hydraulic_conductivity)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(state->bulk_density)); + AllocateBGCMatrixDouble(sizes->ncells_per_col_,sizes->num_components, &(state->total_component_concentration)); + //ALQUIMIA_ASSERT(state->total_mobile.data != NULL); + +} + +void FreeBGCState(BGCState* state) { + if (state != NULL) { + FreeBGCVectorDouble(&(state->liquid_density)); + FreeBGCVectorDouble(&(state->gas_density)); + FreeBGCVectorDouble(&(state->ice_density)); + FreeBGCVectorDouble(&(state->porosity)); + FreeBGCVectorDouble(&(state->water_content)); + FreeBGCVectorDouble(&(state->suction_head)); + FreeBGCVectorDouble(&(state->temperature)); + FreeBGCVectorDouble(&(state->hydraulic_conductivity)); + FreeBGCVectorDouble(&(state->bulk_density)); + FreeBGCMatrixDouble(&(state->total_component_concentration)); + } +} + +void AllocateBGCAuxiliaryData(const BGCSizes* const sizes, BGCAuxiliaryData* aux_data, + int ncells_per_col_) { + AllocateBGCVectorInt(sizes->ncells_per_col_, + &(aux_data->aux_ints)); + + AllocateBGCVectorDouble(sizes->ncells_per_col_, + &(aux_data->aux_doubles)); + +} + +void FreeBGCAuxiliaryData(BGCAuxiliaryData* aux_data) { + if (aux_data != NULL) { + FreeBGCVectorInt(&(aux_data->aux_ints)); + FreeBGCVectorDouble(&(aux_data->aux_doubles)); + } +} + +void AllocateBGCProperties(BGCSizes* sizes, BGCProperties* properties, + int ncells_per_col_) { + sizes->ncells_per_col_ = ncells_per_col_; + + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->liquid_saturation)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->gas_saturation)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->ice_saturation)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->relative_permeability)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->thermal_conductivity)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->volume)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->depth)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->dz)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->plant_wilting_factor)); + AllocateBGCVectorDouble(sizes->ncells_per_col_, &(properties->rooting_depth_fraction)); +} + +void FreeBGCProperties(BGCProperties* properties) { + if (properties != NULL) { + FreeBGCVectorDouble(&(properties->liquid_saturation)); + FreeBGCVectorDouble(&(properties->gas_saturation)); + FreeBGCVectorDouble(&(properties->ice_saturation)); + FreeBGCVectorDouble(&(properties->relative_permeability)); + FreeBGCVectorDouble(&(properties->thermal_conductivity)); + FreeBGCVectorDouble(&(properties->volume)); + FreeBGCVectorDouble(&(properties->depth)); + FreeBGCVectorDouble(&(properties->dz)); + FreeBGCVectorDouble(&(properties->plant_wilting_factor)); + FreeBGCVectorDouble(&(properties->rooting_depth_fraction)); + } +} + +/******************************************************************************* + ** + ** Problem Meta Data + ** + *******************************************************************************/ +/* +void AllocateAlquimiaProblemMetaData(const AlquimiaSizes* const sizes, + AlquimiaProblemMetaData* meta_data) { + + AllocateAlquimiaVectorString(sizes->num_primary, &(meta_data->primary_names)); + ALQUIMIA_ASSERT(meta_data->primary_names.data != NULL); + + AllocateAlquimiaVectorInt(sizes->num_primary, &(meta_data->positivity)); + memset(meta_data->positivity.data, 0, sizeof(int) * sizes->num_primary); + + AllocateAlquimiaVectorString(sizes->num_minerals, + &(meta_data->mineral_names)); + + AllocateAlquimiaVectorString(sizes->num_surface_sites, + &(meta_data->surface_site_names)); + + AllocateAlquimiaVectorString(sizes->num_ion_exchange_sites, + &(meta_data->ion_exchange_names)); + + AllocateAlquimiaVectorString(sizes->num_isotherm_species, + &(meta_data->isotherm_species_names)); + + AllocateAlquimiaVectorString(sizes->num_aqueous_kinetics, + &(meta_data->aqueous_kinetic_names)); + +} //end AllocateAlquimiaProblemMetaData() + +void FreeAlquimiaProblemMetaData(AlquimiaProblemMetaData* meta_data) { + + if (meta_data != NULL) { + FreeAlquimiaVectorString(&(meta_data->primary_names)); + FreeAlquimiaVectorInt(&(meta_data->positivity)); + FreeAlquimiaVectorString(&(meta_data->mineral_names)); + FreeAlquimiaVectorString(&(meta_data->surface_site_names)); + FreeAlquimiaVectorString(&(meta_data->ion_exchange_names)); + FreeAlquimiaVectorString(&(meta_data->isotherm_species_names)); + FreeAlquimiaVectorString(&(meta_data->aqueous_kinetic_names)); + } +} end FreeAlquimiaProblemMetaData() */ + +/******************************************************************************* + ** + ** Auxiliary Output Data + ** + *******************************************************************************/ +/* +void AllocateAlquimiaAuxiliaryOutputData(const AlquimiaSizes* const sizes, + AlquimiaAuxiliaryOutputData* aux_output) { + aux_output->pH = -999.9; + AllocateAlquimiaVectorDouble(sizes->num_minerals, + &(aux_output->mineral_saturation_index)); + + AllocateAlquimiaVectorDouble(sizes->num_aqueous_kinetics, + &(aux_output->aqueous_kinetic_rate)); + + AllocateAlquimiaVectorDouble(sizes->num_minerals, + &(aux_output->mineral_reaction_rate)); + + AllocateAlquimiaVectorDouble(sizes->num_primary, + &(aux_output->primary_free_ion_concentration)); + AllocateAlquimiaVectorDouble(sizes->num_primary, + &(aux_output->primary_activity_coeff)); + + AllocateAlquimiaVectorDouble(sizes->num_aqueous_complexes, + &(aux_output->secondary_free_ion_concentration)); + AllocateAlquimiaVectorDouble(sizes->num_aqueous_complexes, + &(aux_output->secondary_activity_coeff)); + +} // end AllocateAlquimiaAuxiliaryOutputData() + +void FreeAlquimiaAuxiliaryOutputData(AlquimiaAuxiliaryOutputData* aux_output) { + if (aux_output != NULL) { + FreeAlquimiaVectorDouble(&(aux_output->aqueous_kinetic_rate)); + FreeAlquimiaVectorDouble(&(aux_output->mineral_saturation_index)); + FreeAlquimiaVectorDouble(&(aux_output->mineral_reaction_rate)); + FreeAlquimiaVectorDouble(&(aux_output->primary_free_ion_concentration)); + FreeAlquimiaVectorDouble(&(aux_output->primary_activity_coeff)); + FreeAlquimiaVectorDouble(&(aux_output->secondary_free_ion_concentration)); + FreeAlquimiaVectorDouble(&(aux_output->secondary_activity_coeff)); + } +} end FreeAlquimiaAuxiliaryOutputData() */ + +/******************************************************************************* + ** + ** Engine Status + ** + *******************************************************************************/ +/* +void AllocateAlquimiaEngineStatus(AlquimiaEngineStatus* status) { + + status->message = (char*) calloc((size_t)kAlquimiaMaxStringLength, sizeof(char)); + if (NULL == status->message) { + // TODO(bja): error handling + } +} // end AllocateAlquimiaEngineStatus() + +void FreeAlquimiaEngineStatus(AlquimiaEngineStatus* status) { + if (status != NULL) { + free(status->message); + } + status->message = NULL; + +} end FreeAlquimiaEngineStatus() */ + + + +/******************************************************************************* + ** + ** Geochemical conditions/constraints + ** + *******************************************************************************/ +/* +void AllocateAlquimiaGeochemicalConditionVector(const int num_conditions, + AlquimiaGeochemicalConditionVector* condition_list) { + // NOTE: we are only allocating pointers to N conditions here, not + the actual conditions themselves. + fprintf(stdout, " AllocateAlquimiaGeochemicalConditionList() : %d\n", + num_conditions); + condition_list->size = num_conditions; + condition_list->capacity = nearest_power_of_2(num_conditions); + + if (condition_list->size > 0) { + condition_list->data = (AlquimiaGeochemicalCondition*) + calloc((size_t)condition_list->capacity, + sizeof(AlquimiaGeochemicalCondition)); + } +} // end AllocateAlquimiaGeochemicalConditionVector() + +void AllocateAlquimiaGeochemicalCondition(const int size_name, + const int num_aqueous_constraints, + const int num_mineral_constraints, + AlquimiaGeochemicalCondition* condition) { + // NOTE: we are only allocating pointers to N constraints here, not + the actual condstraints themselves. + if (condition != NULL) { + // size_name + 1 to include the null character! + condition->name = (char*) calloc((size_t)size_name+1, sizeof(char)); + AllocateAlquimiaAqueousConstraintVector(num_aqueous_constraints, &condition->aqueous_constraints); + AllocateAlquimiaMineralConstraintVector(num_mineral_constraints, &condition->mineral_constraints); + } +} // end AllocateAlquimiaGeochemicalCondition() + +void AllocateAlquimiaAqueousConstraint(AlquimiaAqueousConstraint* constraint) { + constraint->primary_species_name = + (char*) calloc((size_t)kAlquimiaMaxStringLength, sizeof(char)); + constraint->constraint_type = + (char*) calloc((size_t)kAlquimiaMaxStringLength, sizeof(char)); + constraint->associated_species = + (char*) calloc((size_t)kAlquimiaMaxStringLength, sizeof(char)); + constraint->value = 0.0; +} // end AllocateAlquimiaAqueousConstraint() + +void AllocateAlquimiaAqueousConstraintVector(int num_constraints, + AlquimiaAqueousConstraintVector* constraint_list) { + constraint_list->size = num_constraints; + constraint_list->capacity = nearest_power_of_2(num_constraints); + if (constraint_list->size > 0) { + constraint_list->data = (AlquimiaAqueousConstraint*) + calloc((size_t)constraint_list->capacity, + sizeof(AlquimiaAqueousConstraint)); + } + else + constraint_list->data = NULL; +} + +void AllocateAlquimiaMineralConstraint(AlquimiaMineralConstraint* constraint) { + constraint->mineral_name = + (char*) calloc((size_t)kAlquimiaMaxStringLength, sizeof(char)); + constraint->volume_fraction = -1.0; + constraint->specific_surface_area = -1.0; +} // end AllocateAlquimiaMineralConstraint() + +void AllocateAlquimiaMineralConstraintVector(int num_constraints, + AlquimiaMineralConstraintVector* constraint_list) { + constraint_list->size = num_constraints; + constraint_list->capacity = nearest_power_of_2(num_constraints); + if (constraint_list->size > 0) { + constraint_list->data = (AlquimiaMineralConstraint*) + calloc((size_t)constraint_list->capacity, + sizeof(AlquimiaMineralConstraint)); + } + else + constraint_list->data = NULL; +} + +void FreeAlquimiaGeochemicalConditionVector(AlquimiaGeochemicalConditionVector* condition_list) { + int i; + if (condition_list != NULL) { + for (i = 0; i < condition_list->size; ++i) { + FreeAlquimiaGeochemicalCondition(&(condition_list->data[i])); + } + if (condition_list->data != NULL) { + free(condition_list->data); + condition_list->data = NULL; + } + condition_list->size = 0; + condition_list->capacity = 0; + } +} // end FreeAlquimiaGeochemicalConditionList() + +void FreeAlquimiaGeochemicalCondition(AlquimiaGeochemicalCondition* condition) { + if (condition != NULL) { + if (condition->name != NULL) { + free(condition->name); + condition->name = NULL; + } + FreeAlquimiaAqueousConstraintVector(&(condition->aqueous_constraints)); + FreeAlquimiaMineralConstraintVector(&(condition->mineral_constraints)); + } +} // end FreeAlquimiaGeochemicalCondition() + +void FreeAlquimiaAqueousConstraintVector(AlquimiaAqueousConstraintVector* vector) { + int i; + if (vector != NULL) { + for (i = 0; i < vector->size; ++i) { + FreeAlquimiaAqueousConstraint(&vector->data[i]); + } + if (vector->data != NULL) { + free(vector->data); + vector->data = NULL; + } + vector->size = 0; + vector->capacity = 0; + } +} // end FreeAlquimiaAqueousConstraintVector() + +void FreeAlquimiaAqueousConstraint(AlquimiaAqueousConstraint* constraint) { + free(constraint->primary_species_name); + constraint->primary_species_name = NULL; + free(constraint->constraint_type); + constraint->constraint_type = NULL; + free(constraint->associated_species); + constraint->associated_species = NULL; +} // end FreeAlquimiaAqueousConstraint() + +void FreeAlquimiaMineralConstraintVector(AlquimiaMineralConstraintVector* vector) { + int i; + if (vector != NULL) { + for (i = 0; i < vector->size; ++i) { + FreeAlquimiaMineralConstraint(&vector->data[i]); + } + free(vector->data); + vector->data = NULL; + vector->size = 0; + vector->capacity = 0; + } +} // end FreeAlquimiaMineralConstraintVector() + +void FreeAlquimiaMineralConstraint(AlquimiaMineralConstraint* constraint) { + free(constraint->mineral_name); + constraint->mineral_name = NULL; +} // end FreeAlquimiaMineralConstraint() */ + + +/******************************************************************************* + ** + ** Data convenience struct + ** + *******************************************************************************/ +/* +void AllocateAlquimiaData(AlquimiaData* data) { + AllocateAlquimiaState(&data->sizes, &data->state); + AllocateAlquimiaProperties(&data->sizes, &data->properties); + AllocateAlquimiaAuxiliaryData(&data->sizes, &data->aux_data); + AllocateAlquimiaProblemMetaData(&data->sizes, &data->meta_data); + AllocateAlquimiaAuxiliaryOutputData(&data->sizes, &data->aux_output); +} // end AllocateAlquimiaData() + + +void FreeAlquimiaData(AlquimiaData* data) { + FreeAlquimiaState(&data->state); + FreeAlquimiaProperties(&data->properties); + FreeAlquimiaAuxiliaryData(&data->aux_data); + FreeAlquimiaProblemMetaData(&data->meta_data); + FreeAlquimiaAuxiliaryOutputData(&data->aux_output); +} // end FreeAlquimiaData() */ diff --git a/src/pks/ecosim/data/BGC_memory.hh b/src/pks/ecosim/data/BGC_memory.hh new file mode 100644 index 000000000..b620fdc22 --- /dev/null +++ b/src/pks/ecosim/data/BGC_memory.hh @@ -0,0 +1,134 @@ +/* -*- mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */ + +/* +** Alquimia Copyright (c) 2013-2016, The Regents of the University of California, +** through Lawrence Berkeley National Laboratory (subject to receipt of any +** required approvals from the U.S. Dept. of Energy). All rights reserved. +** +** Alquimia is available under a BSD license. See LICENSE.txt for more +** information. +** +** If you have questions about your rights to use or distribute this software, +** please contact Berkeley Lab's Technology Transfer and Intellectual Property +** Management at TTD@lbl.gov referring to Alquimia (LBNL Ref. 2013-119). +** +** NOTICE. This software was developed under funding from the U.S. Department +** of Energy. As such, the U.S. Government has been granted for itself and +** others acting on its behalf a paid-up, nonexclusive, irrevocable, worldwide +** license in the Software to reproduce, prepare derivative works, and perform +** publicly and display publicly. Beginning five (5) years after the date +** permission to assert copyright is obtained from the U.S. Department of Energy, +** and subject to any subsequent five (5) year renewals, the U.S. Government is +** granted for itself and others acting on its behalf a paid-up, nonexclusive, +** irrevocable, worldwide license in the Software to reproduce, prepare derivative +** works, distribute copies to the public, perform publicly and display publicly, +** and to permit others to do so. +** +** Authors: Benjamin Andre +*/ + +#ifndef BGC_C_MEMORY_H_ +#define BGC_C_MEMORY_H_ + +//#include "alquimia/alquimia_interface.h" +//#include "alquimia/alquimia_containers.h" +#include "BGC_containers.hh" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + /* Alquimia Vectors */ + void AllocateBGCVectorDouble(const int size, BGCVectorDouble* vector); + void FreeBGCVectorDouble(BGCVectorDouble* vector); + + void AllocateBGCVectorInt(const int size, BGCVectorInt* vector); + void FreeBGCVectorInt(BGCVectorInt* vector); + + void AllocateBGCVectorString(const int size, BGCVectorString* vector); + void FreeBGCVectorString(BGCVectorString* vector); + + /* Matrix */ + void AllocateBGCMatrixDouble(const int cells, const int columns, BGCMatrixDouble* matrix); + void FreeBGCMatrixDouble(BGCMatrixDouble* matrix); + + void AllocateBGCMatrixInt(const int cells, const int columns, BGCMatrixInt* matrix); + void FreeBGCMatrixInt(BGCMatrixInt* matrix); + + + void AllocateBGCMatrixString(const int cells, const int columns, BGCMatrixString* matrix); + void FreeBGCMatrixString(BGCMatrixString* matrix); + + void AllocateBGCTensorDouble(const int cells, const int columns, BGCTensorDouble* tensor); + void FreeBGCMatrixDouble(BGCMatrixDouble* tensor); + + void AllocateBGCTensorInt(const int cells, const int columns, BGCTensorInt* tensor); + void FreeBGCTensorInt(BGCTensorInt* tensor); + + /* State */ + void AllocateBGCState(BGCSizes* sizes, + BGCState* state, + int ncells_per_col_, + int num_components, + int num_columns); + void FreeBGCState(BGCState* state); + + /* Auxiliary Data + void AllocateBGCAuxiliaryData(const BGCSizes* const sizes, + BGCAuxiliaryData* aux_data, + int ncells_per_col_); + void FreeBGCAuxiliaryData(BGCAuxiliaryData* aux_data); + */ + /* Properties */ + void AllocateBGCProperties(BGCSizes* sizes, + BGCProperties* properties, + int ncells_per_col_, + int num_columns); + void FreeBGCProperties(BGCProperties* properties); + + // Problem Meta Data + /*void AllocateAlquimiaProblemMetaData(const AlquimiaSizes* const sizes, + AlquimiaProblemMetaData* meta_data); + + void FreeAlquimiaProblemMetaData(AlquimiaProblemMetaData* metda_data); + + // Status + void AllocateAlquimiaEngineStatus(AlquimiaEngineStatus* status); + + void FreeAlquimiaEngineStatus(AlquimiaEngineStatus* status); + + // Auxiliary Output Data + void AllocateAlquimiaAuxiliaryOutputData(const AlquimiaSizes* const sizes, + AlquimiaAuxiliaryOutputData* aux_output); + void FreeAlquimiaAuxiliaryOutputData(AlquimiaAuxiliaryOutputData* aux_output); + + // Geochemical conditions/constraints + void AllocateAlquimiaGeochemicalConditionVector(const int num_conditions, + AlquimiaGeochemicalConditionVector* condition_list); + void AllocateAlquimiaGeochemicalCondition(const int size_name, + const int num_aqueous_constraints, + const int num_mineral_constraints, + AlquimiaGeochemicalCondition* condition); + void AllocateAlquimiaAqueousConstraintVector(const int num_constraints, + AlquimiaAqueousConstraintVector* constraint_list); + void AllocateAlquimiaAqueousConstraint(AlquimiaAqueousConstraint* constraint); + void AllocateAlquimiaMineralConstraintVector(const int num_constraints, + AlquimiaMineralConstraintVector* constraint_list); + void AllocateAlquimiaMineralConstraint(AlquimiaMineralConstraint* constraint); + + void FreeAlquimiaGeochemicalConditionVector(AlquimiaGeochemicalConditionVector* condition_list); + void FreeAlquimiaGeochemicalCondition(AlquimiaGeochemicalCondition* condition); + void FreeAlquimiaAqueousConstraintVector(AlquimiaAqueousConstraintVector* vector); + void FreeAlquimiaAqueousConstraint(AlquimiaAqueousConstraint* constraint); + void FreeAlquimiaMineralConstraintVector(AlquimiaMineralConstraintVector* vector); + void FreeAlquimiaMineralConstraint(AlquimiaMineralConstraint* constraint); + + // Data + void AllocateAlquimiaData(AlquimiaData* data); + void FreeAlquimiaData(AlquimiaData* data);*/ + +#ifdef __cplusplus +} +#endif /* __cplusplus */ + +#endif /* ALQUIMIA_C_MEMORY_H_ */ diff --git a/src/pks/ecosim/data/CMakeLists.txt b/src/pks/ecosim/data/CMakeLists.txt new file mode 100644 index 000000000..37d733eee --- /dev/null +++ b/src/pks/ecosim/data/CMakeLists.txt @@ -0,0 +1,100 @@ +# -*- mode: cmake -*- +# +# ATS +# Data management programs for EcoSIM_ATS +# + +get_property(AMANZI_TPLS_DIR GLOBAL PROPERTY AMANZI_TPLS_DIR) + +set(ECOSIM_INSTALL_PREFIX ${ECOSIM_DIR}/ecosim) +set(ECOSIM_LIB_LOCATION ${ECOSIM_DIR}/ecosim/local/lib) +set(ECOSIM_BUILD_PREFIX ${ECOSIM_DIR}/ecosim/build) +set(NETCDF_LIB ${ECOSIM_DIR}/lib) + +message("In ATS-EcoSIM data CMakeLists:") +message("ECOSIM_DIR:" ${ECOSIM_DIR}) +message("ECOSIM_INSTALL_PREFIX:" ${ECOSIM_INSTALL_PREFIX}) +message("ECOSIM_LIB_LOCATION:" ${ECOSIM_LIB_LOCATION}) +message("ECOSIM_BUILD_PREFIX:" ${ECOSIM_BUILD_PREFIX}) + +message("At include_directories in data") +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Utils/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Minimath/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Modelconfig/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Modelforc/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Mesh/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Modelpars/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Balances/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Ecosim_datatype/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/SoilPhys/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/SurfPhys/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/SnowPhys/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/HydroTherm/PhysData/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Ecosim_mods/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Prescribed_pheno/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/Plant_bgc/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/APIs/) +include_directories(${ECOSIM_BUILD_PREFIX}/f90src/APIData/) + +set(ats_ecosim_data_src_files + BGC_memory.cc + BGC_containers.cc + bgc_fortran_memory_mod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/BGC_containers.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/BGC_containers.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSCPLMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSEcoSIMInitMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/SharedDataMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/c_f_interface_module.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSEcoSIMAdvanceMod.F90 + ${ECOSIM_INSTALL_PREFIX}/f90src/ATSUtils/ATSUtilsMod.F90 +) + +set(ats_ecosim_data_inc_files + BGC_containers.hh + BGC_memory.hh +) + +file(GLOB ECOSIM_LIBRARIES + ${ECOSIM_LIB_LOCATION}/*.a +) + +set(ats_ecosim_data_link_libs + ${Teuchos_LIBRARIES} + ${Epetra_LIBRARIES} + ${ECOSIM_LIBRARIES} + error_handling + atk + mesh + data_structures + whetstone + operators + solvers + time_integration + state + pks + chemistry_pk + ats_pks + ats_eos + ats_operators + ${NETCDF_LIB}/libnetcdf.so + ${NETCDF_LIB}/libnetcdf.so.22 + ${NETCDF_LIB}/libnetcdff.so + ${NETCDF_LIB}/libnetcdff.so.7 + ${NETCDF_LIB}/libnetcdff.so.7.1.0 + gfortran +) + +message(STATUS "ats_ecosim_data_link_libs: ${ats_ecosim_data_link_libs}") + +# make the library +add_amanzi_library(ats_ecosim_data + SOURCE ${ats_ecosim_data_src_files} + HEADERS ${ats_ecosim_data_inc_files} + LINK_LIBS ${ats_ecosim_data_link_libs}) + +generate_evaluators_registration_header( + HEADERFILE ats_ecosim_data_registration.hh + LISTNAME ATS_ECOSIM_DATA_REG + INSTALL True + ) diff --git a/src/pks/ecosim/data/bgc_fortran_memory_mod.F90 b/src/pks/ecosim/data/bgc_fortran_memory_mod.F90 new file mode 100644 index 000000000..323de3bb0 --- /dev/null +++ b/src/pks/ecosim/data/bgc_fortran_memory_mod.F90 @@ -0,0 +1,225 @@ +!Memory and interface mangling functions for the ATS-EcoSIM coupler. + +module bgc_fortran_memory_mod + + + use BGCContainers_module, only : BGCSizes,BGCProperties,& + BGCState,BGCAuxiliaryData + use iso_c_binding, only : c_funptr + + implicit none + + type, bind(C) :: BGCInterface + + type(c_funptr) :: DataTest + type(c_funptr) :: Setup + type(c_funptr) :: Shutdown + type(c_funptr) :: Advance + + end type BGCInterface + + type, public :: BGCFortranInterface + type(BGCInterface) :: c_interface + + contains + procedure, public :: CreateInterface => Create_Fortran_BGC_Interface + procedure, public :: DataTest => BGC_Fortran_DataTest + procedure, public :: Setup => BGC_Fortran_Setup + procedure, public :: Shutdown => BGC_Fortran_Shutdown + procedure, public :: Advance => BGC_Fortran_Advance + end type BGCFortranInterface + + interface + subroutine CreateBGCInterface(engine_name, bgc_interface) bind(C, name='CreateBGCInterface') + use iso_c_binding, only: c_char + IMPORT + implicit none + character(kind=c_char) :: engine_name(*) + type(BGCInterface) :: bgc_interface + end subroutine + end interface + + ! Memory allocation subroutines + + interface + subroutine AllocateBGCState(sizes, state, ncells_per_col_, num_components, num_columns) bind(C, name='AllocateBGCState') + use BGCContainers_module, only : BGCSizes, BGCState + use, intrinsic :: iso_c_binding, only: c_int + implicit none + type(BGCSizes) :: sizes + type(BGCState) :: state + integer(c_int),VALUE :: ncells_per_col_ + integer(c_int),VALUE :: num_components + integer(c_int),VALUE :: num_columns + end subroutine + end interface + interface + subroutine FreeBGCState(state) bind(C, name='FreeBGCState') + use BGCContainers_module, only : BGCState + implicit none + type(BGCState) :: state + end subroutine + end interface + + interface + subroutine AllocateBGCProperties(sizes, properties, ncells_per_col_, num_columns) bind(C, name='AllocateBGCProperties') + use BGCContainers_module, only : BGCSizes, BGCProperties + use, intrinsic :: iso_c_binding, only: c_int + implicit none + type(BGCSizes) :: sizes + type(BGCProperties) :: properties + integer(c_int),VALUE :: ncells_per_col_ + integer(c_int),VALUE :: num_columns + end subroutine + end interface + interface + subroutine FreeBGCProperties(properties) bind(C, name='FreeBGCProperties') + use BGCContainers_module, only : BGCProperties + implicit none + type(BGCProperties) :: properties + end subroutine + end interface + + ! The following subroutines are methods of the engine itself + + interface + subroutine DataTest() bind(C) + use, intrinsic :: iso_c_binding + IMPORT + implicit none + end subroutine + end interface + + + interface + subroutine Setup(properties, state, sizes, num_iterations, num_columns, ncells_per_col_) bind(C) + + use, intrinsic :: iso_c_binding, only: c_char, c_bool, c_ptr, c_int + use BGCContainers_module, only : BGCSizes,BGCProperties,& + BGCState + IMPORT + implicit none + + integer(c_int),VALUE :: num_iterations + integer(c_int),VALUE :: num_columns + integer(c_int),VALUE :: ncells_per_col_ + + type(BGCProperties) :: properties + type(BGCState) :: state + type(BGCSizes) :: sizes + + end subroutine + end interface + + !gracefully shutdown the engine, cleanup memory + interface + subroutine Shutdown() bind(C) + use, intrinsic :: iso_c_binding, only : c_ptr + + implicit none + + end subroutine + end interface + + ! take one (or more?) reaction steps in operator split mode + interface + subroutine Advance(delta_t, properties, state, sizes, num_iterations, num_columns) bind(C) + use, intrinsic :: iso_c_binding, only : c_ptr, c_double, c_int + use BGCContainers_module, only : BGCSizes,BGCProperties,& + BGCState + implicit none + + real(c_double),VALUE :: delta_t + integer(c_int),VALUE :: num_iterations + integer(c_int),VALUE :: num_columns + + type(BGCProperties) :: properties + type(BGCState) :: state + type(BGCSizes) :: sizes + end subroutine + end interface + + contains + + subroutine BGC_Fortran_DataTest(this) + use, intrinsic :: iso_c_binding + + implicit none + class(BGCFortranInterface) :: this + procedure(DataTest), pointer :: engine_DataTest + + call c_f_procpointer(this%c_interface%Setup,engine_DataTest) + call engine_DataTest() + + end subroutine BGC_Fortran_DataTest + + subroutine BGC_Fortran_Setup(this, properties, state, sizes, num_iterations,& + num_columns, ncells_per_col_) + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_double, c_f_procpointer + use BGCContainers_module, only : BGCSizes, BGCProperties,& + BGCState + + implicit none + class(BGCFortranInterface) :: this + + real(c_double) :: delta_t + integer(c_int) :: num_columns + integer(c_int) :: num_iterations + integer(c_int) :: ncells_per_col_ + type(BGCProperties) :: properties + type(BGCState) :: state + type(BGCSizes) :: sizes + + procedure(Setup), pointer :: engine_Setup + + call c_f_procpointer(this%c_interface%Setup,engine_Setup) + call engine_Setup(properties, state, sizes, num_iterations, & + num_columns, ncells_per_col_) + + end subroutine BGC_Fortran_Setup + + subroutine BGC_Fortran_Shutdown(this) + use, intrinsic :: iso_c_binding, only : c_ptr,c_f_procpointer + + implicit none + class(BGCFortranInterface) :: this + procedure(Shutdown), pointer :: engine_Shutdown + + call c_f_procpointer(this%c_interface%Shutdown,engine_Shutdown) + call engine_Shutdown() + + end subroutine BGC_Fortran_Shutdown + + subroutine BGC_Fortran_Advance(this, delta_t, properties, state, sizes, num_iterations, num_columns) + use, intrinsic :: iso_c_binding, only : c_ptr, c_int, c_double, c_f_procpointer + use BGCContainers_module, only : BGCSizes, BGCProperties,& + BGCState + + implicit none + class(BGCFortranInterface) :: this + + real(c_double) :: delta_t + integer(c_int) :: num_columns + integer(c_int) :: num_iterations + type(BGCProperties) :: properties + type(BGCState) :: state + type(BGCSizes) :: sizes + + procedure(Advance), pointer :: engine_Advance + + call c_f_procpointer(this%c_interface%Advance,engine_Advance) + call engine_Advance(delta_t, properties, state, sizes, num_iterations, num_columns) + end subroutine BGC_Fortran_Advance + + subroutine Create_Fortran_BGC_Interface(this,engine_name) + use BGCContainers_module, only :kBGCMaxStringLength + use iso_c_binding, only: c_char,c_null_char + implicit none + class(BGCFortranInterface) :: this + character(kind=c_char,len=kBGCMaxStringLength) :: engine_name + + call CreateBGCInterface(trim(engine_name)//C_NULL_CHAR, this%c_interface) + + end subroutine + +end module bgc_fortran_memory_mod diff --git a/src/pks/ecosim/ecosim_interface.h b/src/pks/ecosim/ecosim_interface.h new file mode 100644 index 000000000..ac4327127 --- /dev/null +++ b/src/pks/ecosim/ecosim_interface.h @@ -0,0 +1,37 @@ + /***************************************************************************** + ** + ** C declarations of the ecosim interface + ** + ******************************************************************************/ + +#include "data/BGC_containers.hh" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +void ecosim_datatest(); + +void ecosim_setup( + BGCProperties* properties, + BGCState* state, + BGCSizes* sizes, + int num_iterations, + int num_columns, + int ncells_per_col_ +); + +void ecosim_shutdown(); + +void ecosim_advance( + double delta_t, + BGCProperties* properties, + BGCState* state, + BGCSizes* sizes, + int num_iterations, + int num_columns +); + +#ifdef __cplusplus +} +#endif /* __cplusplus */