diff --git a/host-configs/LLNL/dane-toss_4_x86_64_ib-gcc@12.1.1.cmake b/host-configs/LLNL/dane-toss_4_x86_64_ib-gcc@12.1.1.cmake index 6cd093dbe08..467b1bef578 100644 --- a/host-configs/LLNL/dane-toss_4_x86_64_ib-gcc@12.1.1.cmake +++ b/host-configs/LLNL/dane-toss_4_x86_64_ib-gcc@12.1.1.cmake @@ -15,6 +15,8 @@ set(CMAKE_C_COMPILER "/usr/tce/packages/gcc/gcc-12.1.1-magic/bin/gcc" CACHE PATH set(CMAKE_CXX_COMPILER "/usr/tce/packages/gcc/gcc-12.1.1-magic/bin/g++" CACHE PATH "") +set(CMAKE_Fortran_COMPILER "/usr/tce/packages/gcc/gcc-12.1.1-magic/bin/gfortran" CACHE PATH "") + set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG" CACHE STRING "") set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -g -DNDEBUG" CACHE STRING "") @@ -37,9 +39,10 @@ set(MPI_C_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-12.1.1-magic/b set(MPI_CXX_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-12.1.1-magic/bin/mpicxx" CACHE PATH "") -set(MPIEXEC "srun" CACHE PATH "") +set(MPI_Fortran_COMPILER "/usr/tce/packages/mvapich2/mvapich2-2.3.7-gcc-12.1.1-magic/bin/mpif90" CACHE PATH "") -set(MPIEXEC_NUMPROC_FLAG "-n" CACHE PATH "") +set(MPIEXEC_EXECUTABLE "/usr/bin/srun" CACHE PATH "") +set(MPIEXEC_NUMPROC_FLAG "-n" CACHE STRING "") #-------------------------------------------------------------------------------- # OpenMP diff --git a/inputFiles/Electrostatics/Electrostatics_alternative.xml b/inputFiles/Electrostatics/Electrostatics_alternative.xml new file mode 100644 index 00000000000..7f823fbd2fa --- /dev/null +++ b/inputFiles/Electrostatics/Electrostatics_alternative.xml @@ -0,0 +1,160 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/Electrostatics/Electrostatics_base.xml b/inputFiles/Electrostatics/Electrostatics_base.xml new file mode 100644 index 00000000000..8ab99d60083 --- /dev/null +++ b/inputFiles/Electrostatics/Electrostatics_base.xml @@ -0,0 +1,157 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/Electrostatics/Electrostatics_variable_reactivity.xml b/inputFiles/Electrostatics/Electrostatics_variable_reactivity.xml new file mode 100644 index 00000000000..ebc351fa868 --- /dev/null +++ b/inputFiles/Electrostatics/Electrostatics_variable_reactivity.xml @@ -0,0 +1,174 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/Electrostatics/LiSEDefect2D_butterfly.jou b/inputFiles/Electrostatics/LiSEDefect2D_butterfly.jou new file mode 100644 index 00000000000..c99895dd445 --- /dev/null +++ b/inputFiles/Electrostatics/LiSEDefect2D_butterfly.jou @@ -0,0 +1,186 @@ +# ============================================================================== +# Battery with Defect 2D +# ============================================================================== +reset + +# Domain width +{L = 6} + +# SE depth +{Ly1 = 1} + +# Cu depth +{Ly2 = 1.25} + +# Defect radius +{R = 1} + +# Defect depth +{h = 0.6} + +# Plane thickness +{thickness=0.01} + + +# ------------------------------------------------------------------------------ +# CREATE COPPER AND DEFECT +# ------------------------------------------------------------------------------ +create brick width {L} depth {Ly2} height {thickness} +volume 1 move X {L/2} Y {Ly1 + Ly2/2} Z {thickness/2} + +create cylinder radius {R} height {thickness} +volume 2 move X {L/2} Y {h} Z {thickness/2} + +# Ratio between the width of the inner square to the defect diameter +{factor = 0.5} + +create brick width {factor*R*2} depth {factor*R*2} height {thickness} +volume 3 move X {L/2} Y {h} Z {thickness/2} + +subtract volume 1 from volume 2 keep +subtract volume 4 from volume 2 + +subtract volume 2 from volume 1 keep +delete volume 1 + +subtract volume 2 from volume 3 keep +subtract volume 6 from volume 3 + +subtract volume 3 from volume 2 keep +delete volume 2 + +merge surface 14 44 +merge surface 24 50 +merge surface 39 43 +merge surface 42 45 + +create surface vertex 12 47 30 31 32 +sweep surface 51 vector 0 0 -1 distance {thickness} +subtract volume 5 from volume 8 keep +subtract volume 7 from volume 8 keep +delete volume 8 + +create surface vertex 13 49 28 34 33 +sweep surface 70 vector 0 0 -1 distance {thickness} +subtract volume 5 from volume 11 keep +subtract volume 7 from volume 11 keep +delete volume 11 + +subtract volume 9 12 from volume 7 keep +delete volume 7 + +subtract volume 10 13 from volume 5 keep +delete volume 5 + +merge surface 14 93 +merge surface 39 81 +merge surface 42 62 +merge surface 58 64 +merge surface 61 92 +merge surface 65 97 +merge surface 77 83 +merge surface 78 89 +merge surface 85 95 +merge surface 91 99 + +# ------------------------------------------------------------------------------ +# CREATE SE +# ------------------------------------------------------------------------------ +create vertex 0 0 {thickness} +create vertex {L} 0 {thickness} + +create surface vertex 114 102 49 47 75 87 132 131 +sweep surface 101 vector 0 0 -1 distance {thickness} + +merge surface 87 109 +merge surface 38 107 +merge surface 80 108 +merge surface 60 106 +merge surface 67 105 + +# ------------------------------------------------------------------------------ +# MESH +# ------------------------------------------------------------------------------ +# No. of intervals in left & right arcs +{Narc_lateral = 5} + +# No. of intervals in top arc +{Narc_top = 40} + +# No. of intervals in left & right radius +{Nrad_lateral = 5} + +# No. of intervals in outter radius +{Nrad_out= 20} + +# No. of intervals in SE depth +{Nse= 10} + +surface all scheme submap + +curve 111 152 interval {Narc_lateral} +curve 116 156 interval {Nrad_lateral} +mesh surface 59 79 + +curve 180 interval {Narc_top} +surface 94 scheme submap +mesh surface 94 + +curve 129 172 interval {Nrad_out} +mesh surface 66 88 + +mesh surface 96 + +mesh surface 40 + +curve 202 interval {Nse} +mesh surface 110 + +volume 14 scheme sweep +mesh volume 14 + +volume 3 9 10 12 13 15 scheme sweep +mesh volume 3 9 10 12 13 15 + +volume 16 scheme sweep +mesh volume 16 + +# ------------------------------------------------------------------------------ +# DEFINE BLOCK IDs +# ------------------------------------------------------------------------------ +block 1 volume 10 13 15 +block 1 name "COPPER" +block 2 volume 3 9 12 14 +block 2 name "DEFECT" +block 3 volume 16 +block 3 name "SE" + +# ------------------------------------------------------------------------------ +# SET BOUNDARY IDs +# ------------------------------------------------------------------------------ +nodeset 1 surface 86 102 +nodeset 1 name "xneg" + +nodeset 2 surface 68 104 +nodeset 2 name "xpos" + +nodeset 3 surface 103 +nodeset 3 name "yneg" + +nodeset 4 surface 100 +nodeset 4 name "ypos" + +nodeset 5 surface 41 63 69 82 84 90 98 101 +nodeset 5 name "zneg" + +nodeset 6 surface 40 59 66 79 88 94 96 110 +nodeset 6 name "zpos" + +nodeset 7 surface 38 60 80 +nodeset 7 name "LISEInterface" + +nodeset 8 surface 67 87 +nodeset 8 name "CUSEInterface" + +export Abaqus "LiSEDefect2D_butterfly.inp" overwrite + diff --git a/inputFiles/Electrostatics/LiSEdefect_Electrostatics_base.xml b/inputFiles/Electrostatics/LiSEdefect_Electrostatics_base.xml new file mode 100644 index 00000000000..c42f799f958 --- /dev/null +++ b/inputFiles/Electrostatics/LiSEdefect_Electrostatics_base.xml @@ -0,0 +1,112 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/Electrostatics/LiSEdefect_Electrostatics_main.xml b/inputFiles/Electrostatics/LiSEdefect_Electrostatics_main.xml new file mode 100644 index 00000000000..218b1b5543a --- /dev/null +++ b/inputFiles/Electrostatics/LiSEdefect_Electrostatics_main.xml @@ -0,0 +1,54 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/solidMechanics/elasticFiniteStrain_cornerPress.xml b/inputFiles/solidMechanics/elasticFiniteStrain_cornerPress.xml new file mode 100644 index 00000000000..edb837cdd4a --- /dev/null +++ b/inputFiles/solidMechanics/elasticFiniteStrain_cornerPress.xml @@ -0,0 +1,149 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/solidMechanics/elasticFiniteStrain_uniaxialPress.xml b/inputFiles/solidMechanics/elasticFiniteStrain_uniaxialPress.xml new file mode 100644 index 00000000000..2c78fb8bc57 --- /dev/null +++ b/inputFiles/solidMechanics/elasticFiniteStrain_uniaxialPress.xml @@ -0,0 +1,148 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/cmake/GeosxOptions.cmake b/src/cmake/GeosxOptions.cmake index 1e504a6ed90..206d0918ddd 100644 --- a/src/cmake/GeosxOptions.cmake +++ b/src/cmake/GeosxOptions.cmake @@ -130,6 +130,7 @@ option( GEOS_ENABLE_SIMPLEPDE "Enables simple PDE physics package" ON ) option( GEOS_ENABLE_SOLIDMECHANICS "Enables solid mechanics physics package" ON ) option( GEOS_ENABLE_SURFACEGENERATION "Enables surface generation physics package" ON ) option( GEOS_ENABLE_WAVEPROPAGATION "Enables wave propagation physics package" ON ) +option( GEOS_ENABLE_ELECTROSTATICS "Enables Butler-Volmer electrostatics physics package" ON ) #set(CMAKE_POSITION_INDEPENDENT_CODE ON CACHE BOOL "" FORCE) #blt_append_custom_compiler_flag(FLAGS_VAR CMAKE_CXX_FLAGS DEFAULT -rdynamic) diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 6af439fe23f..f928774ba20 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -187,6 +187,7 @@ set( constitutive_headers solid/ElasticIsotropicPressureDependent.hpp solid/ElasticTransverseIsotropic.hpp solid/ElasticOrthotropic.hpp + solid/ElasticIsotropicFiniteStrain.hpp solid/InvariantDecompositions.hpp solid/PorousDamageSolid.hpp solid/PerfectlyPlastic.hpp @@ -201,6 +202,7 @@ set( constitutive_headers solid/SolidModelDiscretizationOpsIsotropic.hpp solid/SolidModelDiscretizationOpsTransverseIsotropic.hpp solid/SolidModelDiscretizationOpsOrthotropic.hpp + solid/SolidModelDiscretizationOpsFullTensor.hpp solid/CeramicDamage.hpp solid/SolidFields.hpp solid/porosity/PorosityFields.hpp @@ -218,6 +220,8 @@ set( constitutive_headers thermalConductivity/SinglePhaseThermalConductivityBase.hpp thermalConductivity/SinglePhaseThermalConductivitySelector.hpp thermalConductivity/ThermalConductivityFields.hpp + electroChemistry/ElectroChemistryBase.hpp + electroChemistry/ButlerVolmerReaction.hpp ) # # Specify all sources @@ -332,6 +336,7 @@ set( constitutive_sources solid/ElasticIsotropicPressureDependent.cpp solid/ElasticTransverseIsotropic.cpp solid/ElasticOrthotropic.cpp + solid/ElasticIsotropicFiniteStrain.cpp solid/PorousDamageSolid.cpp solid/PerfectlyPlastic.cpp solid/PorousSolid.cpp @@ -349,6 +354,8 @@ set( constitutive_sources thermalConductivity/MultiPhaseVolumeWeightedThermalConductivity.cpp thermalConductivity/SinglePhaseThermalConductivity.cpp thermalConductivity/SinglePhaseThermalConductivityBase.cpp + electroChemistry/ElectroChemistryBase.cpp + electroChemistry/ButlerVolmerReaction.cpp ) set( dependencyList ${parallelDeps} functions denseLinearAlgebra ) diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp index 1c63ae20f7f..6802f40c0b8 100644 --- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp +++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp @@ -36,6 +36,7 @@ #include "solid/ElasticIsotropicPressureDependent.hpp" #include "solid/ElasticTransverseIsotropic.hpp" #include "solid/ElasticOrthotropic.hpp" +#include "solid/ElasticIsotropicFiniteStrain.hpp" #include "solid/PorousSolid.hpp" #include "solid/PorousDamageSolid.hpp" #include "solid/CompressibleSolid.hpp" @@ -55,6 +56,8 @@ #include "permeability/WillisRichardsPermeability.hpp" #include "contact/CoulombFriction.hpp" #include "contact/RateAndStateFriction.hpp" +#include "electroChemistry/ElectroChemistryBase.hpp" +#include "electroChemistry/ButlerVolmerReaction.hpp" namespace geos @@ -164,6 +167,7 @@ struct ConstitutivePassThru< SolidBase > ModifiedCamClay, DelftEgg, DruckerPrager, + ElasticIsotropicFiniteStrain, ElasticIsotropic, ElasticTransverseIsotropic, ElasticIsotropicPressureDependent, @@ -237,6 +241,7 @@ struct ConstitutivePassThruTriaxialDriver< SolidBase > ModifiedCamClay, DelftEgg, DruckerPrager, + ElasticIsotropicFiniteStrain, ElasticIsotropic, ElasticTransverseIsotropic, ElasticIsotropicPressureDependent, @@ -537,6 +542,36 @@ struct ConstitutivePassThru< CoupledSolidBase > } }; +/** + * Base material model for electrochemistry. + */ +template<> +struct ConstitutivePassThru +{ + template< typename LAMBDA > + static + void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler::execute(constitutiveRelation, + std::forward(lambda)); + } +}; + +/** + * @brief Material model for interface Butler-Volmer kinetics + */ +template<> +struct ConstitutivePassThru +{ + template< typename LAMBDA > + static + void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler::execute(constitutiveRelation, + std::forward(lambda)); + } +}; + } /* namespace constitutive */ } /* namespace geos */ diff --git a/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.cpp b/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.cpp new file mode 100644 index 00000000000..4630ff015f7 --- /dev/null +++ b/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.cpp @@ -0,0 +1,57 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ButlerVolmerReaction.cpp + */ + +#include "constitutive/electroChemistry/ButlerVolmerReaction.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +ButlerVolmerInterface::ButlerVolmerInterface(string const& name, Group * const parent) +: +ConstitutiveBase(name, parent), +m_krxn() +{ + registerWrapper(viewKeyStruct::defaultReactivityCoefficientString(), &m_defaultKrxn). + setApplyDefaultValue(1.0). + setInputFlag(InputFlags::REQUIRED). + setDescription("Default BV Reactivity Coefficient"); + + registerWrapper(viewKeyStruct::reactivityCoefficientString(), &m_krxn). + setApplyDefaultValue(-1.0). // will be overwritten + setPlotLevel(PlotLevel::LEVEL_0). + setDescription("Reactivity Coefficient Field"); +} + +ButlerVolmerInterface::~ButlerVolmerInterface() {} + +void ButlerVolmerInterface::postInputInitialization() +{ + this->getWrapper>(viewKeyStruct::reactivityCoefficientString()). + setApplyDefaultValue(m_defaultKrxn); +} + +REGISTER_CATALOG_ENTRY(ConstitutiveBase, ButlerVolmerInterface, string const&, Group * const) +} // namespace constitutive + +} // namespace geos \ No newline at end of file diff --git a/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.hpp b/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.hpp new file mode 100644 index 00000000000..4ea96047092 --- /dev/null +++ b/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.hpp @@ -0,0 +1,126 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ButlerVolmerReaction.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_BUTLERVOLMERREACTION_HPP +#define GEOS_CONSTITUTIVE_BUTLERVOLMERREACTION_HPP + +#include "common/DataLayouts.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "constitutive/ConstitutiveBase.hpp" + +namespace geos +{ + +namespace constitutive +{ + +/** + * @brief The abstract base class for reaction coefficient update + */ +class ButlerVolmerReactionUpdate +{ +public: + + /** + * @brief Constructor for the class performing the reactivity coefficient updates + * @param k_rxn the array of face-wise reactivity coefficient in the subregion + */ + ButlerVolmerReactionUpdate(arrayView1d const& k_rxn) + : + m_krxn(k_rxn) + {} + + /** + * @brief Number of elements storing reactivity coefficient data + * @return Number of elements + */ + localIndex numElem() const + { + return m_krxn.size(); + } + + /** + * @brief Get reactivity coefficient + * @param[in] k index of the element in the subRegion + * @return the reactivity coefficient of element k + */ + GEOS_HOST_DEVICE + inline + real64 getReactCoeff(localIndex const k) const + { + return m_krxn[k]; + } + +protected: + + /// View on the cell-wise conductivities + arrayView1d m_krxn; +}; + +/** + * @brief The class for interface Butler-Volmer kinetics + */ +class ButlerVolmerInterface : public ConstitutiveBase +{ +public: + + /** + * @brief Constructor for the abstract base class + * @param[in] name the name of the class + * @param[in] parent pointer to the parent Group + */ + ButlerVolmerInterface( string const & name, dataRepository::Group * const parent ); + + /** + * Destructor + */ + virtual ~ButlerVolmerInterface() override; + + static string catalogName() { return "ButlerVolmerInterface"; } + + virtual string getCatalogName() const override { return catalogName(); } + + /// Keys for data in this class + struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct + { + static constexpr char const* reactivityCoefficientString() {return "reactivityCoefficient";} + static constexpr char const* defaultReactivityCoefficientString() {return "defaultReactivityCoefficient";} + }; + + using KernelWrapper = ButlerVolmerReactionUpdate; + KernelWrapper createKernelUpdates() + { + return KernelWrapper(m_krxn); + } + +protected: + + /// Post-process XML input + virtual void postInputInitialization() override; + + array1d m_krxn; + + real64 m_defaultKrxn = 1.0; +}; + +} // namespace constitutive + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_BUTLERVOLMERREACTION_HPP \ No newline at end of file diff --git a/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.cpp b/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.cpp new file mode 100644 index 00000000000..a8f71a6e220 --- /dev/null +++ b/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.cpp @@ -0,0 +1,57 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ElectroChemistryBase.cpp + */ + +#include "constitutive/electroChemistry/ElectroChemistryBase.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace constitutive +{ + +ElectroChemistryBase::ElectroChemistryBase(string const& name, Group * const parent) +: +ConstitutiveBase(name, parent), +m_conductivity() +{ + registerWrapper(viewKeyStruct::defaultConductivityString(), &m_defaultConductivity). + setApplyDefaultValue(1.0). + setInputFlag(InputFlags::REQUIRED). + setDescription("Default Electro Conductivity"); + + registerWrapper(viewKeyStruct::conductivityString(), &m_conductivity). + setApplyDefaultValue(-1.0). // will be overwritten + setPlotLevel(PlotLevel::LEVEL_0). + setDescription("Electro Conductivity Field"); +} + +ElectroChemistryBase::~ElectroChemistryBase() {} + +void ElectroChemistryBase::postInputInitialization() +{ + this->getWrapper>(viewKeyStruct::conductivityString()). + setApplyDefaultValue(m_defaultConductivity); +} + +REGISTER_CATALOG_ENTRY(ConstitutiveBase, ElectroChemistryBase, string const&, Group * const) +} // namespace constitutive + +} // namespace geos \ No newline at end of file diff --git a/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.hpp b/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.hpp new file mode 100644 index 00000000000..602ecdef3f6 --- /dev/null +++ b/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.hpp @@ -0,0 +1,128 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ElectroChemistryBase.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_ELECTROCHEMISTRYBASE_HPP +#define GEOS_CONSTITUTIVE_ELECTROCHEMISTRYBASE_HPP + +#include "common/DataLayouts.hpp" +#include "common/GEOS_RAJA_Interface.hpp" +#include "constitutive/ConstitutiveBase.hpp" + +namespace geos +{ + +namespace constitutive +{ + +/** + * @brief The abstract base class for electrochemistry material + */ +class ElectroChemistryBaseUpdate +{ +public: + + /** + * @brief Constructor for the class performing the electro conductivity updates + * @param conductivity the array of cell-wise conductivities in the subregion + */ + ElectroChemistryBaseUpdate(arrayView1d const& conductivity) + : + m_conductivity(conductivity) + {} + + /** + * @brief Number of elements storing conductivity data + * @return Number of elements + */ + localIndex numElem() const + { + return m_conductivity.size(); + } + + /** + * @brief Get conductivity + * @param[in] k index of the cell in the subRegion + * @return the conductivity of element k + */ + GEOS_HOST_DEVICE + inline + real64 getConductivity(localIndex const k) const + { + return m_conductivity[k]; + } + +protected: + + /// View on the cell-wise conductivities + arrayView1d m_conductivity; +}; + +/** + * @brief The abstract base class for electro conductivity + */ +class ElectroChemistryBase : public ConstitutiveBase +{ +public: + + /** + * @brief Constructor for the abstract base class + * @param[in] name the name of the class + * @param[in] parent pointer to the parent Group + */ + ElectroChemistryBase( string const & name, dataRepository::Group * const parent ); + + /** + * Destructor + */ + virtual ~ElectroChemistryBase() override; + + static string catalogName() { return "ElectroChemistryBase"; } + + virtual string getCatalogName() const override { return catalogName(); } + + /// Keys for data in this class + struct viewKeyStruct : public ConstitutiveBase::viewKeyStruct + { + static constexpr char const* conductivityString() {return "conductivity";} + static constexpr char const* defaultConductivityString() {return "defaultConductivity";} + }; + + // virtual void allocateConstitutiveData( dataRepository::Group & parent, + // localIndex const numConstitutivePointsPerParentIndex ) override; + + using KernelWrapper = ElectroChemistryBaseUpdate; + KernelWrapper createKernelUpdates() + { + return KernelWrapper(m_conductivity); + } + +protected: + + /// Post-process XML input + virtual void postInputInitialization() override; + + array1d m_conductivity; + + real64 m_defaultConductivity = 1.0; +}; +} // namespace constitutive + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_ELECTROCHEMISTRYBASE_HPP \ No newline at end of file diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.cpp b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.cpp new file mode 100644 index 00000000000..7b569498747 --- /dev/null +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.cpp @@ -0,0 +1,43 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ElasticIsotropicFiniteStrain.cpp + */ + +#include "ElasticIsotropicFiniteStrain.hpp" + +namespace geos +{ +using namespace dataRepository; + +namespace constitutive +{ + +ElasticIsotropicFiniteStrain::ElasticIsotropicFiniteStrain(string const & name, Group * const parent) + : ElasticIsotropic( name, parent ) {} + +ElasticIsotropicFiniteStrain::~ElasticIsotropicFiniteStrain() +{} + +void ElasticIsotropicFiniteStrain::postInputInitialization() +{ + ElasticIsotropic::postInputInitialization(); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, ElasticIsotropicFiniteStrain, std::string const &, Group * const ) + +} // namespace constitutive +} // namespace geos \ No newline at end of file diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.hpp new file mode 100644 index 00000000000..1439464400a --- /dev/null +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.hpp @@ -0,0 +1,426 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ElasticIsotropicFiniteStrain.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPICFINITESTRAIN_HPP_ +#define GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPICFINITESTRAIN_HPP_ + +#include "ElasticIsotropic.hpp" +#include "PropertyConversions.hpp" +#include "SolidModelDiscretizationOpsFullTensor.hpp" +#include "LvArray/src/tensorOps.hpp" + +namespace geos{ + +namespace constitutive { + +/** + * @class ElasticIsotropicFiniteStrainUpdates + * + * Class to provide material updates that may be + * called from a kernel function. + */ +class ElasticIsotropicFiniteStrainUpdates : public ElasticIsotropicUpdates +{ +public: + ElasticIsotropicFiniteStrainUpdates( + arrayView1d< real64 const > const & bulkModulus, + arrayView1d< real64 const > const & shearModulus, + arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView3d< real64, solid::STRESS_USD > const & newStress, + arrayView3d< real64, solid::STRESS_USD > const & oldStress, + bool const & disableInelasticity ) + : ElasticIsotropicUpdates(bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity) + {} + + /// Default copy constructor + ElasticIsotropicFiniteStrainUpdates( ElasticIsotropicFiniteStrainUpdates const & ) = default; + + /// Default move constructor + ElasticIsotropicFiniteStrainUpdates( ElasticIsotropicFiniteStrainUpdates && ) = default; + + /// Deleted default constructor + ElasticIsotropicFiniteStrainUpdates() = delete; + + /// Deleted copy assignment operator + ElasticIsotropicFiniteStrainUpdates & operator=( ElasticIsotropicFiniteStrainUpdates const & ) = delete; + + /// Deleted move assignment operator + ElasticIsotropicFiniteStrainUpdates & operator=( ElasticIsotropicFiniteStrainUpdates && ) = delete; + + /// Use the uncompressed version of the stiffness bilinear form + using DiscretizationOps = SolidModelDiscretizationOpsFullTensor; + + // returns first Piola-Kirchhoff stress which is asymmetric + GEOS_HOST_DEVICE + void finiteStrainNoStateUpdate_StressOnly(localIndex const k, localIndex const q, + real64 const (&totalElasticStrain)[6], + real64 const (&fInv)[3][3], + real64 (&firstPiolaStress)[3][3], + real64 (&kirchhoffStress)[6] + ) const; + + GEOS_HOST_DEVICE + void finiteStrainUpdate_StressOnly(localIndex const k, localIndex const q, + real64 const (&totalElasticStrain)[6], + real64 const (&fInv)[3][3], + real64 (&firstPiolaStress)[3][3] + ) const; + + GEOS_HOST_DEVICE + void finiteStrainNoStateUpdate(localIndex const k, localIndex const q, + real64 const (&elasticDeformGrad)[3][3], + real64 (&firstPiolaStress)[3][3], + real64 (&kirchhoffStress)[6], + real64 (&stiffness)[9][9]) const; + + GEOS_HOST_DEVICE + void finiteStrainUpdate(localIndex const k, localIndex const q, + real64 const (&elasticDeformGrad)[3][3], + real64 (&firstPiolaStress)[3][3], + real64 (&stiffness)[9][9]) const; + + GEOS_HOST_DEVICE + void finiteStrainUpdate(localIndex const k, localIndex const q, + real64 const (&elasticDeformGrad)[3][3], + real64 (&firstPiolaStress)[3][3], + DiscretizationOps & stiffness) const; + + GEOS_HOST_DEVICE + void computeLogElasticStrain(real64 const (&elasticDeformGrad)[3][3], + real64 (&elasticStrain)[6], + real64 (&eigenValues)[3], + real64 (&eigenVectors)[3][3], + real64 (&eigenVectorsT)[3][3]) const; + +private: + GEOS_HOST_DEVICE + void computeMaterialTangentColumn(real64 const (&deltaCe)[6], real64 const (&M_hat)[6], + real64 const (&Q)[3][3], real64 const (&Q_T)[3][3], + real64 const (&dtau_dEe)[6][6], real64 const (&fInvT)[3][3], + real64 (&deltaP_mat)[3][3]) const; +}; + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::computeLogElasticStrain( + real64 const (&elasticDeformGrad)[3][3], real64 (&elasticStrain)[6], + real64 (&eigenValues)[3], real64 (&eigenVectors)[3][3], real64 (&eigenVectorsT)[3][3]) const +{ + real64 C_e[3][3] = {}; // Right Cauchy tensor + LvArray::tensorOps::Rij_eq_AkiAkj<3, 3>(C_e, elasticDeformGrad); + + real64 Ce_symmetric[6] = {}; + LvArray::tensorOps::denseToSymmetric<3>(Ce_symmetric, C_e); + + // initial eigen vector are stored in rows + LvArray::tensorOps::symEigenvectors<3>(eigenValues, eigenVectorsT, Ce_symmetric); + LvArray::tensorOps::transpose<3, 3>(eigenVectors, eigenVectorsT); + + // matrix logarithm + real64 logLambda[6] = { {0} }; + for( int i = 0; i < 3; i++ ) + { + logLambda[i] = 0.5 * log(eigenValues[i]); + } + + // elastic strain + // E_e = 0.5 * Q * ln(Lam) * Q^T + LvArray::tensorOps::Rij_eq_AikSymBklAjl<3>(elasticStrain, eigenVectors, logLambda); + + // scale shear components by 2.0 to use in small strain material model where stiffness is expressed in Voigt notation + elasticStrain[3] *= 2.0; + elasticStrain[4] *= 2.0; + elasticStrain[5] *= 2.0; +} + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::computeMaterialTangentColumn( + real64 const (&deltaCe)[6], real64 const (&M_hat)[6], real64 const (&Q)[3][3], + real64 const (&Q_T)[3][3], real64 const (&dtau_dEe)[6][6], real64 const (&fInvT)[3][3], + real64 (&deltaP_mat)[3][3]) const +{ + // Q^T * deltaCe * Q + real64 deltaCe_hat[6] = {}; + LvArray::tensorOps::Rij_eq_AikSymBklAjl<3>(deltaCe_hat, Q_T, deltaCe); + + // M \circ Q^T * deltaCe * Q + real64 deltaLnCe_hat[6] = {}; + LvArray::tensorOps::hadamardProduct<6>(deltaLnCe_hat, M_hat, deltaCe_hat); + + // derivative of matrix logarithm deltaEe = dEe / dFe = 1/2 * Q * (M \circ Q^T * deltaCe * Q) * Q^T + real64 deltaLnCe[6] = {}; + real64 deltaEe[6] = {}; + LvArray::tensorOps::Rij_eq_AikSymBklAjl<3>(deltaLnCe, Q, deltaLnCe_hat); + LvArray::tensorOps::copy<6>(deltaEe, deltaLnCe); + LvArray::tensorOps::scale<6>(deltaEe, 0.5); + + // scale shear components for stiffness in voigt notation + deltaEe[3] *= 2.0; + deltaEe[4] *= 2.0; + deltaEe[5] *= 2.0; + + // double contraction dtau / dEe * deltaEe + real64 dtau_voigt[6] = {}; + real64 deltaTau[3][3] = {}; + LvArray::tensorOps::Ri_eq_AijBj<6, 6>(dtau_voigt, dtau_dEe, deltaEe); + LvArray::tensorOps::symmetricToDense<3>(deltaTau, dtau_voigt); + + // deltaTau * F^-T + LvArray::tensorOps::Rij_eq_AikBkj<3, 3, 3>(deltaP_mat, deltaTau, fInvT); +} + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::finiteStrainNoStateUpdate_StressOnly(localIndex const k, localIndex const q, + real64 const (&totalElasticStrain)[6], + real64 const (&fInv)[3][3], + real64 (&firstPiolaStress)[3][3], + real64 (&kirchhoffStress)[6] + ) const +{ + this->smallStrainNoStateUpdate_StressOnly(k, q, totalElasticStrain, kirchhoffStress); + + LvArray::tensorOps::Rij_eq_symAikBjk<3>(firstPiolaStress, kirchhoffStress, fInv); +} + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::finiteStrainUpdate_StressOnly(localIndex const k, localIndex const q, + real64 const (&totalElasticStrain)[6], + real64 const (&fInv)[3][3], + real64 (&firstPiolaStress)[3][3] + ) const +{ + real64 kirchhoffStress[6] = {}; + finiteStrainNoStateUpdate_StressOnly(k, q, totalElasticStrain, fInv, firstPiolaStress, kirchhoffStress); + + // Can only save the symmetric kirchhoff stress right now + LvArray::tensorOps::copy<6>(m_oldStress[k][q], m_newStress[k][q]); + LvArray::tensorOps::copy<6>(m_newStress[k][q], kirchhoffStress); +} + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::finiteStrainNoStateUpdate(localIndex const k, localIndex const q, + real64 const (&elasticDeformGrad)[3][3], + real64 (&firstPiolaStress)[3][3], + real64 (&kirchhoffStress)[6], + real64 (&stiffness)[9][9]) const +{ + real64 fInv[3][3] = {}; + real64 fInvT[3][3] = {}; + + real64 totalElasticStrain[6] = {}; + + real64 eigenValues[3] = {}; + real64 eigenVectors[3][3] = {}; + real64 eigenVectorsT[3][3] = {}; + + LvArray::tensorOps::invert<3>(fInv, elasticDeformGrad); + LvArray::tensorOps::transpose<3, 3>(fInvT, fInv); + + computeLogElasticStrain(elasticDeformGrad, totalElasticStrain, eigenValues, eigenVectors, eigenVectorsT); + + finiteStrainNoStateUpdate_StressOnly(k, q, totalElasticStrain, fInv, firstPiolaStress, kirchhoffStress); + + real64 smallStrainStiffness[6][6] = {}; + this->getElasticStiffness(k, q, smallStrainStiffness); + + // derivative of log matrix: \partial ln(C_e) / \partial C_e, C_e = F^T F + // build symmetric spectral of log matrix derivative + real64 M_hat[6] = {}; + M_hat[0] = 1.0 / eigenValues[0]; + M_hat[1] = 1.0 / eigenValues[1]; + M_hat[2] = 1.0 / eigenValues[2]; + + if (abs(eigenValues[1] - eigenValues[2]) <= 1e-12) { + M_hat[3] = 1.0 / eigenValues[1]; + } else { + M_hat[3] = log(eigenValues[1] / eigenValues[2]) / (eigenValues[1] - eigenValues[2]); + } + + if (abs(eigenValues[0] - eigenValues[2]) <= 1e-12) { + M_hat[4] = 1.0 / eigenValues[0]; + } else { + M_hat[4] = log(eigenValues[0] / eigenValues[2]) / (eigenValues[0] - eigenValues[2]); + } + + if (abs(eigenValues[0] - eigenValues[1]) <= 1e-12) { + M_hat[5] = 1.0 / eigenValues[0]; + } else { + M_hat[5] = log(eigenValues[0] / eigenValues[1]) / (eigenValues[0] - eigenValues[1]); + } + + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + real64 deltaFe[3][3] = {0}; + deltaFe[i][j] = 1.0; + + // delCe = delFe^T * Fe + Fe^T * delFe + real64 dCe[3][3] = {0}; + for (int l = 0; l < 3; ++l) { + dCe[j][l] += deltaFe[i][j] * elasticDeformGrad[i][l]; + dCe[l][j] += elasticDeformGrad[i][l] * deltaFe[i][j]; + } + + real64 deltaCe[6] = {}; + LvArray::tensorOps::denseToSymmetric<3>(deltaCe, dCe); + + // material tangent + real64 deltaP_mat[3][3] = {}; + computeMaterialTangentColumn(deltaCe, M_hat, eigenVectors, eigenVectorsT, smallStrainStiffness, fInvT, deltaP_mat); + + // geometric tangent + real64 deltaP_geo[3][3] = {}; + real64 deltaFeT_FeInvT[3][3] = {0}; + for (int l = 0; l < 3; ++l) { + deltaFeT_FeInvT[j][l] += deltaFe[i][j] * fInvT[i][l]; + } + LvArray::tensorOps::Rij_eq_AikBkj<3, 3, 3>(deltaP_geo, firstPiolaStress, deltaFeT_FeInvT); + + real64 deltaP_total[3][3] = {}; + LvArray::tensorOps::copy<3, 3>(deltaP_total, deltaP_mat); + LvArray::tensorOps::scaledAdd<3, 3>(deltaP_total, deltaP_geo, -1.0); + + // input results into D operator + int col = 3 * i + j; + for (int m = 0; m < 3; ++m) { + for (int n = 0; n < 3; ++n) { + // row major flattening + int row = 3 * m + n; + stiffness[row][col] = deltaP_total[m][n]; + } + } + + } + } + +} + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::finiteStrainUpdate(localIndex const k, localIndex const q, + real64 const (&elasticDeformGrad)[3][3], + real64 (&firstPiolaStress)[3][3], + real64 (&stiffness)[9][9]) const +{ + real64 kirchhoffStress[6] = {}; + finiteStrainNoStateUpdate(k, q, elasticDeformGrad, firstPiolaStress, kirchhoffStress, stiffness); + + // Can only save the symmetric kirchhoff stress right now + LvArray::tensorOps::copy<6>(m_oldStress[k][q], m_newStress[k][q]); + LvArray::tensorOps::copy<6>(m_newStress[k][q], kirchhoffStress); +} + +GEOS_HOST_DEVICE +inline +void ElasticIsotropicFiniteStrainUpdates::finiteStrainUpdate(localIndex const k, localIndex const q, + real64 const (&elasticDeformGrad)[3][3], + real64 (&firstPiolaStress)[3][3], + DiscretizationOps & stiffness) const +{ + finiteStrainUpdate(k, q, elasticDeformGrad, firstPiolaStress, stiffness.m_c); +} + +/** + * @class ElasticIsotropicFiniteStrain + * + * Finite strain isotropic elastic material model. + */ +class ElasticIsotropicFiniteStrain : public ElasticIsotropic +{ +public: + /// @typedef Alias for ElasticIsotropicFiniteStrainUpdates + using KernelWrapper = ElasticIsotropicFiniteStrainUpdates; + + /** + * constructor + * @param[in] name name of the instance in the catalog + * @param[in] parent the group which contains this instance + */ + ElasticIsotropicFiniteStrain( string const & name, Group * const parent ); + + /** + * Default Destructor + */ + virtual ~ElasticIsotropicFiniteStrain() override; + + /// string name to use for this class in the catalog + static constexpr auto m_catalogNameString = "ElasticIsotropicFiniteStrain"; + + /** + * @brief Static catalog string + * @return A string that is used to register/lookup this class in the registry + */ + static std::string catalogName() { return m_catalogNameString; } + + /** + * @brief Get catalog name + * @return Name string + */ + virtual string getCatalogName() const override { return catalogName(); } + + /** + * @brief Create a instantiation of the ElasticFiniteStrainIsotropicUpdate class + * that refers to the data in this. + * @param includeState Flag whether to pass state arrays that may not be needed for "no-state" updates + * @return An instantiation of ElasticFiniteStrainIsotropicUpdate. + */ + KernelWrapper createKernelUpdates(bool const includeState = true) const + { + if (includeState) { + return KernelWrapper(m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_newStress, m_oldStress, m_disableInelasticity); + } else { + return KernelWrapper(m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + arrayView3d< real64, solid::STRESS_USD >(), + arrayView3d< real64, solid::STRESS_USD >(), + m_disableInelasticity); + } + } + + /** + * @brief Construct an update kernel for a derived type. + * @tparam UPDATE_KERNEL The type of update kernel from the derived type. + * @tparam PARAMS The parameter pack to hold the constructor parameters for the derived update kernel. + * @param constructorParams The constructor parameter for the derived type. + * @return An @p UPDATE_KERNEL object. + */ + template< typename UPDATE_KERNEL, typename ... PARAMS > + UPDATE_KERNEL createDerivedKernelUpdates( PARAMS && ... constructorParams ) + { + return UPDATE_KERNEL( std::forward< PARAMS >( constructorParams )..., + m_bulkModulus, + m_shearModulus, + m_thermalExpansionCoefficient, + m_newStress, + m_oldStress, + m_disableInelasticity ); + } + +protected: + virtual void postInputInitialization() override; +}; + +} // namespace constitutive + +} // namespace geos + +#endif /* GEOS_CONSTITUTIVE_SOLID_ELASTICISOTROPIC_HPP_ */ diff --git a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullTensor.hpp b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullTensor.hpp new file mode 100644 index 00000000000..4b9e5c21325 --- /dev/null +++ b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullTensor.hpp @@ -0,0 +1,81 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file SolidModelDiscretizationOpsFullTensor.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_SOLID_SOLIDMODELDISCRETIAZTIONOPSFULLTENSOR_HPP_ +#define GEOS_CONSTITUTIVE_SOLID_SOLIDMODELDISCRETIAZTIONOPSFULLTENSOR_HPP_ + +#include "SolidModelDiscretizationOps.hpp" + +namespace geos +{ +namespace constitutive +{ + +struct SolidModelDiscretizationOpsFullTensor : public SolidModelDiscretizationOps +{ + template< int NUM_SUPPORT_POINTS, + typename BASIS_GRADIENT > + GEOS_HOST_DEVICE + void BTDB( BASIS_GRADIENT const & gradN, + real64 const & detJxW, + real64 ( &elementStiffness )[NUM_SUPPORT_POINTS*3][NUM_SUPPORT_POINTS*3] ); + + GEOS_HOST_DEVICE + inline + void scaleParams( real64 const scale ) + { + LvArray::tensorOps::scale< 9, 9 >( m_c, scale ); + } + + real64 m_c[9][9]; +}; + +template< int NUM_SUPPORT_POINTS, + typename BASIS_GRADIENT > +GEOS_HOST_DEVICE +inline +void SolidModelDiscretizationOpsFullTensor::BTDB( BASIS_GRADIENT const & gradN, + real64 const & detJxW, + real64 (& elementStiffness)[NUM_SUPPORT_POINTS * 3][NUM_SUPPORT_POINTS * 3] ) +{ + real64 B[9][NUM_SUPPORT_POINTS * 3] = { {0} }; + + for (int a = 0; a < NUM_SUPPORT_POINTS; ++a) { + for (int i = 0; i < 3; ++i) { + int col = 3 * a + i; + for (int j = 0; j < 3; ++j) { + int row = 3 * i + j; + B[row][col] = gradN[a][j]; + } + } + } + + real64 DB[9][NUM_SUPPORT_POINTS * 3]; + real64 elmat[NUM_SUPPORT_POINTS * 3][NUM_SUPPORT_POINTS * 3]; + + LvArray::tensorOps::Rij_eq_AikBkj<9, NUM_SUPPORT_POINTS * 3, 9>(DB, m_c, B); + LvArray::tensorOps::Rij_eq_AkiBkj(elmat, B, DB); + LvArray::tensorOps::scale(elmat, detJxW); + LvArray::tensorOps::add(elementStiffness, elmat); +} + +} // namespace constitutive +} // namespace geos + +#endif /* GEOS_CONSTITUTIVE_SOLID_SOLIDMODELDISCRETIAZTIONOPSFULLTENSOR_HPP_ */ diff --git a/src/coreComponents/constitutive/solid/SolidUtilities.hpp b/src/coreComponents/constitutive/solid/SolidUtilities.hpp index 337b7f30509..9f2dd4601a2 100644 --- a/src/coreComponents/constitutive/solid/SolidUtilities.hpp +++ b/src/coreComponents/constitutive/solid/SolidUtilities.hpp @@ -157,6 +157,143 @@ struct SolidUtilities } } + /** + * @brief Perform a finite-difference check of the stiffness computation + * + * This method uses several stress evaluations and finite differencing to + * approximate the 9x9 stiffness matrix, and then computes an error between + * the coded stiffness method and the finite difference version. + * + * @note This method only works for models providing the finiteStrainUpdate + * method returning a 9x9 stiffness. + * + * @param solid the solid kernel wrapper + * @param k the element number + * @param q the quadrature index + * @param elasticDeformGrad elastic deformation gradient (on top of which a FD perturbation will be added) + * @param print flag to decide if debug output is printed or not + */ + template< typename SOLID_TYPE > + GEOS_HOST_DEVICE + static bool + checkFiniteStrainStiffness( SOLID_TYPE const & solid, + localIndex k, + localIndex q, + real64 const ( &elasticDeformGrad )[3][3], + bool print = false ) + { + real64 firstPiolaStress[3][3] = {}; + real64 stiffness[9][9] = {}; + solid.finiteStrainUpdate(0, 0, elasticDeformGrad, firstPiolaStress, stiffness); + + real64 stiffnessFD[9][9] = {}; + SolidUtilities::computeFiniteStrainFiniteDifferenceStiffness(solid, k, q, elasticDeformGrad, stiffnessFD); + + real64 error = 0; + real64 norm = 0; + real64 rerr = 0; + + for( localIndex i = 0; i < 9; ++i ) + { + for( localIndex j = 0; j < 9; ++j ) + { + error += pow( stiffnessFD[i][j] - stiffness[i][j], 2.0 ); + norm += pow( stiffness[i][j], 2.0 ); + } + } + error = sqrt(error); + norm = sqrt(norm); + rerr = error / norm; + + if( print ) + { + for( localIndex i = 0; i < 9; ++i ) + { + for( localIndex j = 0; j < 9; ++j ) + { + // printf( "[%12.5e vs %12.5e] ", stiffnessFD[i][j], stiffness[i][j] ); + printf( "%12.5e ", fabs(stiffnessFD[i][j] - stiffness[i][j]) ); + } + printf( "\n" ); + } + + printf("Abs err = %12.5e, Rel err = %12.5e\n", error, rerr); + } + + return (rerr < 1e-3); + } + + /** + * @brief Perform a finite-difference stiffness computation for finite strain material model + * + * This method uses stress evaluations and finite differencing to + * approximate the 9x9 stiffness matrix. + * + * @note This method only works for models providing the finiteStrainUpdate + * method returning a 9x9 stiffness, as it will primarily be used to check + * the hand coded tangent against a finite difference reference. + * + * @param solid the solid kernel wrapper + * @param k the element number + * @param q the quadrature index + * @param elasticDeformGrad elastic deformation gradient (on top of which a FD perturbation will be added) + * @param stiffnessFD finite different stiffness approximation + */ + template< typename SOLID_TYPE > + GEOS_HOST_DEVICE + static void + computeFiniteStrainFiniteDifferenceStiffness( SOLID_TYPE const & solid, + localIndex k, + localIndex q, + real64 const ( &elasticDeformGrad )[3][3], + real64 ( & stiffnessFD )[9][9] ) + { + real64 eps = 1e-7; + real64 F[3][3] = {}; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + F[i][j] = elasticDeformGrad[i][j]; + } + } + + real64 fInv[3][3] = {}; + real64 totalElasticStrain[6] = {}; + real64 eigenValues[3] = {}; + real64 eigenVectors[3][3] = {}; + real64 eigenVectorsT[3][3] = {}; + + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + // plus + real64 P_plus[3][3] = {}; + F[i][j] += eps; + LvArray::tensorOps::invert<3>(fInv, F); + + solid.computeLogElasticStrain(F, totalElasticStrain, eigenValues, eigenVectors, eigenVectorsT); + solid.finiteStrainUpdate_StressOnly(k, q, totalElasticStrain, fInv, P_plus); + F[i][j] -= eps; + + // minus + real64 P_minus[3][3] = {}; + F[i][j] -= eps; + LvArray::tensorOps::invert<3>(fInv, F); + + solid.computeLogElasticStrain(F, totalElasticStrain, eigenValues, eigenVectors, eigenVectorsT); + solid.finiteStrainUpdate_StressOnly(k, q, totalElasticStrain, fInv, P_minus); + F[i][j] += eps; + + int col = 3 * i + j; + for (int m = 0; m < 3; ++m) { + for (int n = 0; n < 3; ++n) { + // row major flattening + int row = 3 * m + n; + stiffnessFD[row][col] = (P_plus[m][n] - P_minus[m][n]) / (2 * eps); + } + } + } + } + } + /** * @brief Hypo update (small strain, large rotation). * diff --git a/src/coreComponents/constitutive/solid/unitTests/CMakeLists.txt b/src/coreComponents/constitutive/solid/unitTests/CMakeLists.txt index 36b456cb4b5..999f31cad17 100644 --- a/src/coreComponents/constitutive/solid/unitTests/CMakeLists.txt +++ b/src/coreComponents/constitutive/solid/unitTests/CMakeLists.txt @@ -1,9 +1,10 @@ # Unit tests for constitutive solid models (DamageSpectral etc.) set( gtest_geosx_tests - testDamage.cpp ) + testDamage.cpp + testFiniteDifferenceElasticFiniteStrain.cpp ) -set( dependencyList gtest blas lapack constitutive ${parallelDeps} ) +set( dependencyList gtest mainInterface blas lapack constitutive ${parallelDeps} ) if( ENABLE_CUDA AND ENABLE_CUDA_NVTOOLSEXT ) list( APPEND dependencyList CUDA::nvToolsExt ) diff --git a/src/coreComponents/constitutive/solid/unitTests/testFiniteDifferenceElasticFiniteStrain.cpp b/src/coreComponents/constitutive/solid/unitTests/testFiniteDifferenceElasticFiniteStrain.cpp new file mode 100644 index 00000000000..8bed576896b --- /dev/null +++ b/src/coreComponents/constitutive/solid/unitTests/testFiniteDifferenceElasticFiniteStrain.cpp @@ -0,0 +1,109 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include + +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutive/solid/ElasticIsotropicFiniteStrain.hpp" +#include "constitutive/solid/SolidUtilities.hpp" +#include "dataRepository/xmlWrapper.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" + +using namespace geos; +using namespace ::geos::constitutive; + +TEST( ElasticFiniteStrainTests, testMaterialTangentFiniteDifference ) +{ + // create a model and test xml input + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + ConstitutiveManager constitutiveManager( "constitutive", &rootGroup ); + + string const inputStream = + "" + " " + ""; + + xmlWrapper::xmlDocument xmlDocument; + xmlWrapper::xmlResult xmlResult = xmlDocument.loadString( inputStream ); + if( !xmlResult ) + { + GEOS_LOG_RANK_0( "XML parsed with errors!" ); + GEOS_LOG_RANK_0( "Error description: " << xmlResult.description()); + GEOS_LOG_RANK_0( "Error offset: " << xmlResult.offset ); + } + + xmlWrapper::xmlNode xmlConstitutiveNode = xmlDocument.getChild( "Constitutive" ); + constitutiveManager.processInputFileRecursive( xmlDocument, xmlConstitutiveNode ); + constitutiveManager.postInputInitializationRecursive(); + + localIndex constexpr numElem = 1; + localIndex constexpr numQuad = 1; + + dataRepository::Group disc( "discretization", &rootGroup ); + disc.resize( numElem ); + + ElasticIsotropicFiniteStrain& constitutive_model = constitutiveManager.getConstitutiveRelation< ElasticIsotropicFiniteStrain >( "lithium" ); + constitutive_model.allocateConstitutiveData( disc, numQuad ); + + // confirm allocation sizes + EXPECT_EQ( constitutive_model.size(), numElem ); + EXPECT_EQ( constitutive_model.numQuadraturePoints(), numQuad ); + + ElasticIsotropicFiniteStrain::KernelWrapper cmw = constitutive_model.createKernelUpdates(); + real64 deformationGradient[3][3] = {}; + deformationGradient[0][0] = 1.05; + deformationGradient[0][1] = 0.02; + deformationGradient[0][2] = 0.01; + deformationGradient[1][0] = 0.0; + deformationGradient[1][1] = 1.02; + deformationGradient[1][2] = 0.0; + deformationGradient[2][0] = 0.0; + deformationGradient[2][1] = 0.0; + deformationGradient[2][2] = 1.01; + + EXPECT_TRUE( SolidUtilities::checkFiniteStrainStiffness( cmw, 0, 0, deformationGradient, true ) ); + + deformationGradient[0][0] = 1.5; + deformationGradient[0][1] = 2.0; + deformationGradient[0][2] = 3.0; + deformationGradient[1][0] = 1.2; + deformationGradient[1][1] = 1.3; + deformationGradient[1][2] = 4.0; + deformationGradient[2][0] = 2.5; + deformationGradient[2][1] = 1.2; + deformationGradient[2][2] = 1.4; + + EXPECT_TRUE( SolidUtilities::checkFiniteStrainStiffness( cmw, 0, 0, deformationGradient, true ) ); +} + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::GeosxState state( geos::basicSetup( argc, argv ) ); + + int const result = RUN_ALL_TESTS(); + + geos::basicCleanup(); + + return result; +} \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index 246fc46be5f..f22c3af7639 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -118,6 +118,11 @@ if( GEOS_ENABLE_MULTIPHYSICS ) list( APPEND dependencyList multiPhysicsSolvers ) endif() +if( GEOS_ENABLE_ELECTROSTATICS ) + add_subdirectory( electrostatics ) + list( APPEND dependencyList electrostaticsSolvers ) +endif() + geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) diff --git a/src/coreComponents/physicsSolvers/electrostatics/CMakeLists.txt b/src/coreComponents/physicsSolvers/electrostatics/CMakeLists.txt new file mode 100644 index 00000000000..8f458c002df --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/CMakeLists.txt @@ -0,0 +1,45 @@ +# SPDX-License-Identifier: LGPL-2.1-only +# +# Copyright (c) 2016-2024 Lawrence Livermore National Security LLC +# Copyright (c) 2018-2024 TotalEnergies +# Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University +# Copyright (c) 2023-2024 Chevron +# Copyright (c) 2019- GEOS/GEOSX Contributors +# All rights reserved +# +# See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. +# +#-------------------------------------------------------------------------------------------------- + +#[[ +Package: Electrostatics +Contains: +#]] + +# Specify solver headers +set(electrostaticsSolvers_headers + Electrostatics.hpp + ElectrostaticsKernels.hpp) + +set(electrostaticsSolvers_sources + Electrostatics.cpp) + +set( dependencyList ${parallelDeps} physicsSolversBase ) + +geos_decorate_link_dependencies( LIST decoratedDependencies DEPENDENCIES ${dependencyList} ) + +blt_add_library( NAME electrostaticsSolvers + SOURCES ${electrostaticsSolvers_sources} + HEADERS ${electrostaticsSolvers_headers} + DEPENDS_ON ${decoratedDependencies} ${externalComponentDeps} + OBJECT ${GEOS_BUILD_OBJ_LIBS} + SHARED ${GEOS_BUILD_SHARED_LIBS} + ) + +target_include_directories( electrostaticsSolvers PUBLIC ${CMAKE_SOURCE_DIR}/coreComponents ) + +install( TARGETS electrostaticsSolvers LIBRARY DESTINATION ${CMAKE_INSTALL_PREFIX}/lib ) + +if( externalComponentDeps ) + target_include_directories( electrostaticsSolvers PUBLIC ${CMAKE_SOURCE_DIR}/externalComponents ) +endif() \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.cpp b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.cpp new file mode 100644 index 00000000000..f6dff6ae335 --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.cpp @@ -0,0 +1,382 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file Electrostatics.cpp + */ + +#include "Electrostatics.hpp" +#include "ElectrostaticsKernels.hpp" + +#include "constitutive/electroChemistry/ElectroChemistryBase.hpp" +#include "constitutive/electroChemistry/ButlerVolmerReaction.hpp" + +#include "fieldSpecification/FieldSpecificationManager.hpp" +#include "fieldSpecification/TractionBoundaryCondition.hpp" + +#include "mesh/DomainPartition.hpp" +#include "mesh/FaceElementSubRegion.hpp" +#include "mesh/CellElementSubRegion.hpp" +#include "mesh/mpiCommunications/NeighborCommunicator.hpp" + +namespace geos +{ + +using namespace dataRepository; +using namespace constitutive; + +Electrostatics::Electrostatics(const string& name, Group* const parent) +: +PhysicsSolverBase(name, parent), +m_fieldName("primaryField"), +m_timeIntegrationOption(TimeIntegrationOption::QuasiStatic) +{ + registerWrapper(viewKeyStruct::timeIntegrationOption(), &m_timeIntegrationOption). + setApplyDefaultValue(m_timeIntegrationOption). + setInputFlag(InputFlags:: OPTIONAL). + setDescription("Time integration method. Options are:\n* " + EnumStrings< TimeIntegrationOption >::concat( "\n*" )); + + registerWrapper(viewKeyStruct::fieldVarName(), &m_fieldName). + setRTTypeName(rtTypes::CustomTypes::groupNameRef). + setInputFlag(InputFlags::REQUIRED). + setDescription("Name of field variable"); + + registerWrapper(viewKeyStruct::surfaceGeneratorNameString(), &m_surfaceGeneratorName). + setInputFlag(InputFlags::OPTIONAL). + setDescription("Name of the surface generator to use"); +} + +Electrostatics::~Electrostatics() {} + +void Electrostatics::postInputInitialization() +{ + PhysicsSolverBase::postInputInitialization(); + + m_surfaceGenerator = this->getParent().getGroupPointer(m_surfaceGeneratorName); +} + +void Electrostatics::registerDataOnMesh(Group& meshBodies) +{ + forDiscretizationOnMeshTargets(meshBodies, + [&](string const&, MeshLevel& mesh, string_array const& regionNames) { + NodeManager& nodes = mesh.getNodeManager(); + + nodes.registerWrapper(m_fieldName) + .setApplyDefaultValue(0.0) + .setPlotLevel(PlotLevel::LEVEL_0) + .setDescription("Primary field variable"); + + ElementRegionManager& elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< CellElementSubRegion >(regionNames, + [&]( localIndex const, CellElementSubRegion& subRegion) { + subRegion.registerWrapper(viewKeyStruct::electroMaterialNamesString()). + setPlotLevel(PlotLevel::NOPLOT). + setRestartFlags(RestartFlags::NO_WRITE). + setSizedFromParent(0); + + string& electroMaterialName = subRegion.getReference(viewKeyStruct::electroMaterialNamesString()); + electroMaterialName = PhysicsSolverBase::getConstitutiveName(subRegion); + GEOS_ERROR_IF( electroMaterialName.empty(), GEOS_FMT("{}: ElectroChemistryBase model not found on subregion {}", + getDataContext(), subRegion.getName())); + }); + + elemManager.forElementSubRegions< FaceElementSubRegion >(regionNames, + [&]( localIndex const, FaceElementSubRegion& subRegion) { + subRegion.registerWrapper(viewKeyStruct::reactiveMaterialNamesString()). + setPlotLevel(PlotLevel::NOPLOT). + setRestartFlags(RestartFlags::NO_WRITE). + setSizedFromParent(0); + + string& reactiveMaterialName = subRegion.getReference(viewKeyStruct::reactiveMaterialNamesString()); + reactiveMaterialName = PhysicsSolverBase::getConstitutiveName(subRegion); + GEOS_ERROR_IF( reactiveMaterialName.empty(), GEOS_FMT("{}: ButlerVolmerInterface model not found on subregion {}", + getDataContext(), subRegion.getName())); + }); + }); +} + +real64 Electrostatics::solverStep(real64 const& time_n, real64 const& dt, + const int cycleNumber, DomainPartition& domain) +{ + GEOS_MARK_FUNCTION; + real64 dtReturn = dt; + + // if (m_surfaceGenerator != nullptr && cycleNumber == 0) + // { + // m_surfaceGenerator->solverStep(time_n, dt, cycleNumber, domain); + // } + + // setupSystem(domain, m_dofManager, m_localMatrix, m_rhs, m_solution, false); + if (cycleNumber == 0) { + FieldSpecificationManager& fieldSpecificationManager = FieldSpecificationManager::getInstance(); + forDiscretizationOnMeshTargets(domain.getMeshBodies(), + [&](string const&, MeshLevel& mesh, string_array const&) { + fieldSpecificationManager.applyInitialConditions(mesh); + }); + } + + implicitStepSetup(time_n, dt, domain); + + dtReturn = linearImplicitStep(time_n, dt, cycleNumber, domain); + + return dtReturn; +} + +void Electrostatics::implicitStepSetup(real64 const& GEOS_UNUSED_PARAM(time_n), + real64 const& GEOS_UNUSED_PARAM(dt), + DomainPartition& domain) +{ + Timestamp const meshModificationTimestamp = getMeshModificationTimestamp(domain); + + if (meshModificationTimestamp > getSystemSetupTimestamp()) + { + setupSystem(domain, m_dofManager, m_localMatrix, m_rhs, m_solution); + setSystemSetupTimestamp(meshModificationTimestamp); + } +} + +void Electrostatics::setupDofs(DomainPartition const& GEOS_UNUSED_PARAM(domain), + DofManager& dofManager) const +{ + GEOS_MARK_FUNCTION; + dofManager.addField(m_fieldName, FieldLocation::Node, 1, getMeshTargets()); + dofManager.addCoupling(m_fieldName, m_fieldName, DofManager::Connector::Elem); +} + +void Electrostatics::setupSystem(DomainPartition& domain, DofManager& dofManager, + CRSMatrix< real64, globalIndex >& localMatrix, + ParallelVector& rhs, ParallelVector& solution, + bool const GEOS_UNUSED_PARAM(setSparsity)) +{ + GEOS_LOG("Electrostatics::setupSystem"); + + GEOS_MARK_FUNCTION; + PhysicsSolverBase::setupSystem(domain, dofManager, localMatrix, rhs, solution, true); +} + +void Electrostatics::assembleSystem(real64 const GEOS_UNUSED_PARAM(time_n), real64 const dt, + DomainPartition& domain, DofManager const& dofManager, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs) +{ + GEOS_MARK_FUNCTION; + + localMatrix.zero(); + localRhs.zero(); + + forDiscretizationOnMeshTargets(domain.getMeshBodies(), + [&](string const&, MeshLevel& mesh, string_array const& regionNames) { + NodeManager& nodeManager = mesh.getNodeManager(); + string const dofKey = dofManager.getKey(m_fieldName); + arrayView1d const& dofIndex = nodeManager.getReference>(dofKey); + + ElectrostaticsKernelFactory kernelFactory(dofIndex, dofManager.rankOffset(), localMatrix, localRhs, dt, m_fieldName); + + finiteElement::regionBasedKernelApplication, constitutive::ElectroChemistryBase, CellElementSubRegion>( + mesh, regionNames, this->getDiscretizationName(), viewKeyStruct::electroMaterialNamesString(), kernelFactory); + }); + + applyButlerVolmerCurrent(dofManager, domain, localMatrix, localRhs); +} + +void Electrostatics::applyBoundaryConditions(real64 const time_n, real64 const dt, + DomainPartition& domain, DofManager const& dofManager, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs) +{ + applyPotentialBC(time_n + dt, dofManager, domain, localMatrix, localRhs); + + applyCurrentBC(time_n + dt, dofManager, domain, localRhs); +} + +// real64 Electrostatics::calculateResidualNorm(real64 const& GEOS_UNUSED_PARAM(time_n), real64 const& GEOS_UNUSED_PARAM(dt), +// DomainPartition const& domain, DofManager const& dofManager, +// arrayView1d const& localRhs) +// { +// GEOS_MARK_FUNCTION; +// +// real64 totalResidualNorm = 0.0; +// +// forDiscretizationOnMeshTargets(domain.getMeshBodies(), +// [&](string const&, MeshLevel const& mesh, string_array const&){ +// NodeManager const& nodeManager = mesh.getNodeManager(); +// string const dofKey = dofManager.getKey(m_fieldName); +// arrayView1d const dofNumber = nodeManager.getReference>(dofKey); +// +// globalIndex const rankOffset = dofManager.rankOffset(); +// arrayView1d const ghostRank = nodeManager.ghostRank(); +// +// RAJA::ReduceSum localSum(0.0); +// +// SortedArrayView const& targetNodes = nodeManager.sets() +// }); +// } + +void Electrostatics::applySystemSolution(DofManager const& dofManager, arrayView1d const& localSolution, + real64 const scalingFactor, real64 const dt, DomainPartition& domain) +{ + GEOS_UNUSED_VAR(dt); + dofManager.addVectorToField(localSolution, m_fieldName, m_fieldName, scalingFactor); + + forDiscretizationOnMeshTargets(domain.getMeshBodies(), + [&](string const&, MeshLevel& mesh, string_array const&) { + FieldIdentifiers fieldsToBeSync; + fieldsToBeSync.addFields(FieldLocation::Node, {m_fieldName}); + + CommunicationTools::getInstance().synchronizeFields(fieldsToBeSync, mesh, domain.getNeighbors(), true); + }); +} + +void Electrostatics::applyPotentialBC(real64 const time, DofManager const& dofManager, DomainPartition& domain, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs) +{ + FieldSpecificationManager const& fsManager = FieldSpecificationManager::getInstance(); + + forDiscretizationOnMeshTargets(domain.getMeshBodies(), + [&](string const&, MeshLevel& mesh, string_array const&) + { + fsManager.apply(time, mesh, m_fieldName, + [&](FieldSpecificationBase const& bc, string const&, SortedArrayView const& targetSet, + NodeManager& targetGroup, string const& GEOS_UNUSED_PARAM(fieldName)) + { + bc.applyBoundaryConditionToSystem>( + targetSet, time, targetGroup, m_fieldName, dofManager.getKey(m_fieldName), + dofManager.rankOffset(), localMatrix, localRhs); + }); + }); +} + +void Electrostatics::applyCurrentBC(real64 const time, DofManager const& dofManager, + DomainPartition& domain, arrayView1d const& localRhs) +{ + FieldSpecificationManager& fsManager = FieldSpecificationManager::getInstance(); + + forDiscretizationOnMeshTargets(domain.getMeshBodies(), + [&](string const&, MeshLevel& mesh, string_array const&) + { + FaceManager const& faceManager = mesh.getFaceManager(); + NodeManager const& nodeManager = mesh.getNodeManager(); + + string const dofKey = dofManager.getKey(m_fieldName); + arrayView1d const blockLocalDofNumber = nodeManager.getReference(dofKey); + globalIndex const dofRankOffset = dofManager.rankOffset(); + + fsManager.template apply(time, mesh, TractionBoundaryCondition::catalogName(), + [&](TractionBoundaryCondition const& bc, string const&, SortedArrayView const& targetSet, + Group&, string const&) + { + bc.launch(time, blockLocalDofNumber, dofRankOffset, faceManager, targetSet, localRhs); + }); + }); +} + +void Electrostatics::applyButlerVolmerCurrent(DofManager const& dofManager, DomainPartition& domain, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs) +{ + GEOS_MARK_FUNCTION; + + forDiscretizationOnMeshTargets(domain.getMeshBodies(), + [&](string const&, MeshLevel& mesh, string_array const&) + { + FaceManager const& faceManager = mesh.getFaceManager(); + NodeManager& nodeManager = mesh.getNodeManager(); + ElementRegionManager& elemManager = mesh.getElemManager(); + + arrayView1d const phi = nodeManager.getReference>(m_fieldName).toViewConst(); + + arrayView2d const faceNormal = faceManager.faceNormal(); + ArrayOfArraysView const facesToNodes = faceManager.nodeList().toViewConst(); + + string const dofKey = dofManager.getKey(m_fieldName); + arrayView1d const nodeDofNumber = nodeManager.getReference(dofKey); + globalIndex const rankOffset = dofManager.rankOffset(); + + constexpr localIndex maxNodesPerFace = 4; + constexpr localIndex maxDofPerElem = maxNodesPerFace * 2; + + elemManager.forElementSubRegions([&](FaceElementSubRegion& subRegion) { + string const & constitutiveName = subRegion.getReference(viewKeyStruct::reactiveMaterialNamesString()); + constitutive::ConstitutiveBase& constitutiveModel = subRegion.getConstitutiveModel( constitutiveName ); + constitutive::ButlerVolmerInterface& castedConstitutiveModel = dynamic_cast(constitutiveModel); + constitutive::ButlerVolmerInterface::KernelWrapper const m_constitutiveUpdate(castedConstitutiveModel.createKernelUpdates()); + + arrayView1d const area = subRegion.getElementArea(); + arrayView2d const elemsToFaces = subRegion.faceList().toViewConst(); + + forAll>(subRegion.size(), [=](localIndex const kfe) { + real64 const k_rxn = m_constitutiveUpdate.getReactCoeff(kfe); + + localIndex const kf0 = elemsToFaces[kfe][0], kf1 = elemsToFaces[kfe][1]; + + localIndex const numNodesPerFace = facesToNodes.sizeOfArray(kf0); + real64 const Ja = area[kfe] / numNodesPerFace; + + stackArray1d< globalIndex, maxDofPerElem > rowDof(numNodesPerFace*2); + stackArray1d< real64, maxDofPerElem > nodeRHS(numNodesPerFace*2); + stackArray2d< real64, maxDofPerElem *maxDofPerElem > dRdPhi(numNodesPerFace*2, numNodesPerFace*2); + + for (localIndex a = 0; a < numNodesPerFace; ++a) + { + localIndex const node0 = facesToNodes[kf0][a]; + localIndex const node1 = facesToNodes[kf1][a == 0 ? a : numNodesPerFace - a]; + + real64 phi_jump = phi[node0] - phi[node1]; + + rowDof[a] = nodeDofNumber[node0]; + rowDof[numNodesPerFace + a] = nodeDofNumber[node1]; + + // The factor 2.0 comes from linearizing BV with Taylor expansion + nodeRHS[a] += 2.0 * k_rxn * Ja * phi_jump / thermodynamicPotential; + nodeRHS[numNodesPerFace + a] -= 2.0 * k_rxn * Ja * phi_jump / thermodynamicPotential; + + // initial implementation with mass lumping + dRdPhi(a, a) += 2.0 * k_rxn * Ja / thermodynamicPotential; + dRdPhi(a, numNodesPerFace + a) -= 2.0 * k_rxn * Ja / thermodynamicPotential; + dRdPhi(numNodesPerFace + a, numNodesPerFace + a) += 2.0 * k_rxn * Ja / thermodynamicPotential; + dRdPhi(numNodesPerFace + a, a) -= 2.0 * k_rxn * Ja / thermodynamicPotential; + } + + for (localIndex idof = 0; idof < numNodesPerFace * 2; ++idof) + { + localIndex const localRow = LvArray::integerConversion(rowDof[idof] - rankOffset); + if (localRow >= 0 && localRow < localMatrix.numRows()) + { + localMatrix.addToRowBinarySearchUnsorted( + localRow, rowDof.data(), dRdPhi[idof].dataIfContiguous(), numNodesPerFace*2); + RAJA::atomicAdd(&localRhs[localRow], nodeRHS[idof]); + } + } + }); + }); + }); +} + +void Electrostatics::updateState(DomainPartition& domain) +{ + GEOS_UNUSED_VAR(domain); +} + +void Electrostatics::resetStateToBeginningOfStep(DomainPartition& GEOS_UNUSED_PARAM(domain)) {} + +void Electrostatics::implicitStepComplete(real64 const& GEOS_UNUSED_PARAM(time_n), + real64 const& GEOS_UNUSED_PARAM(dt), + DomainPartition& GEOS_UNUSED_PARAM(domain)) +{} + +REGISTER_CATALOG_ENTRY(PhysicsSolverBase, Electrostatics, string const&, dataRepository::Group* const) +} \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.hpp b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.hpp new file mode 100644 index 00000000000..320e037f7cd --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.hpp @@ -0,0 +1,148 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file Electrostatics.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_ELECTROSTATICS_HPP_ +#define GEOS_PHYSICSSOLVERS_ELECTROSTATICS_HPP_ + + +#include "common/format/EnumStrings.hpp" +#include "common/TimingMacros.hpp" +#include "mesh/mpiCommunications/CommunicationTools.hpp" +#include "mesh/mpiCommunications/MPI_iCommData.hpp" +#include "physicsSolvers/PhysicsSolverBase.hpp" + +namespace geos +{ + +/** + * @class Electrostatics + * + * This class implements a finite element solution to the charge balance equation with interface reaction. + */ +class Electrostatics : public PhysicsSolverBase +{ +public: + /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the + /// tree structure of classes) + Electrostatics(const string& name, Group* const parent); + + /// Destructor + virtual ~Electrostatics() override; + + /// "CatalogName()" return the string used as XML tag in the input file. It ties the XML tag with + /// this C++ classes. This is important. + static string catalogName() { return "Electrostatics"; } + + /** + * @copydoc PhysicsSolverBase::getCatalogName() + */ + string getCatalogName() const override { return catalogName(); } + + // virtual void initializePreSubGroups() override; + + /// This method ties properties with their supporting mesh + virtual void registerDataOnMesh(Group & meshBodies) override final; + + virtual real64 solverStep(real64 const& time_n, real64 const& dt, + integer const cycleNumber, + DomainPartition& domain) override; + + virtual void implicitStepSetup(real64 const& time_n, real64 const& dt, + DomainPartition& domain) override; + + virtual void setupDofs(DomainPartition const& domain, DofManager& dofManager) const override; + + virtual void setupSystem(DomainPartition& domain, DofManager& dofManager, + CRSMatrix& localMatrix, + ParallelVector& rhs, ParallelVector& solution, + bool const setSparsity = false) override; + + virtual void assembleSystem(real64 const time, real64 const dt, + DomainPartition& domain, DofManager const& dofManager, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs) override; + + virtual void applyBoundaryConditions(real64 const time, real64 const dt, + DomainPartition& domain, DofManager const& dofManager, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs) override; + + // virtual real64 calculateResidualNorm(real64 const& time_n, real64 const& dt, + // DomainPartition const& domain, DofManager const& dofManager, + // arrayView1d const& localRhs) override; + + virtual void applySystemSolution(DofManager const& dofManager, + arrayView1d const& localSolution, + real64 const scalingFactor, real64 const dt, + DomainPartition& domain) override; + + virtual void updateState(DomainPartition& domain) override final; + + virtual void resetStateToBeginningOfStep(DomainPartition& GEOS_UNUSED_PARAM(domain)) override; + + virtual void implicitStepComplete(real64 const& time, real64 const& dt, + DomainPartition& domain) override; + + void applyPotentialBC(real64 const time, DofManager const& dofManager, DomainPartition& domain, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs); + + void applyCurrentBC(real64 const time, DofManager const& dofManager, + DomainPartition& domain, arrayView1d const& localRhs); + + void applyButlerVolmerCurrent(DofManager const& dofManager, DomainPartition& domain, + CRSMatrixView const& localMatrix, + arrayView1d const& localRhs); + + enum class TimeIntegrationOption : integer + { + QuasiStatic, + ImplicitTransient + }; + + struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct + { + static constexpr char const* timeIntegrationOption() { return "timeIntegrationOption"; } + static constexpr char const* fieldVarName() { return "fieldName"; } + static constexpr char const* electroMaterialNamesString() {return "electroMaterialNames";} + static constexpr char const* reactiveMaterialNamesString() {return "reactiveMaterialNames";} + static constexpr char const* surfaceGeneratorNameString() {return "surfaceGeneratorName";} + }; + +protected: + virtual void postInputInitialization() override; + + string m_fieldName; + TimeIntegrationOption m_timeIntegrationOption; + + static constexpr double F = 96485.33212; // C/mol, Faraday's constant + static constexpr double T = 298.15; // K, Temperature + static constexpr double R = 8.314462618; // J/mol/K, ideal gas constant + real64 thermodynamicPotential = R * T / F; // Volt + +private: + PhysicsSolverBase *m_surfaceGenerator; + string m_surfaceGeneratorName; +}; + +ENUM_STRINGS(Electrostatics::TimeIntegrationOption, + "QuasiStatic", "ImplicitTransient"); +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_ELECTROSTATICS_HPP_ */ \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/electrostatics/ElectrostaticsKernels.hpp b/src/coreComponents/physicsSolvers/electrostatics/ElectrostaticsKernels.hpp new file mode 100644 index 00000000000..e378dc8955c --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/ElectrostaticsKernels.hpp @@ -0,0 +1,209 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + + +/** + * @file ElectrostaticsKernels.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_ELECTROSTATICS_ELECTROSTATICSKERNELS_HPP_ +#define GEOS_PHYSICSSOLVERS_ELECTROSTATICS_ELECTROSTATICSKERNELS_HPP_ + +#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" +#include "finiteElement/kernelInterface/SparsityKernelBase.hpp" + +namespace geos +{ +//***************************************************************************** +/** + * @brief Implements kernels for solving electrostatics equation. + * @copydoc geos::finiteElement::KernelBase + * @tparam NUM_NODES_PER_ELEM The number of nodes per element for the + * @p SUBREGION_TYPE. + * @tparam UNUSED An unused parameter since we are assuming that the test and + * trial space have the same number of support points. + * + * ### ElectrostaticsKernel Description + * Implements the KernelBase interface functions required for solving electrostatics + * equation using one of the finite element kernel application functions such as + * geos::finiteElement::RegionBasedKernelApplication. + * + * In this implementation, the template parameter @p NUM_NODES_PER_ELEM is used + * in place of both @p NUM_TEST_SUPPORT_POINTS_PER_ELEM and + * @p NUM_TRIAL_SUPPORT_POINTS_PER_ELEM, which are assumed to be equal. This + * results in the @p UNUSED template parameter as only the NUM_NODES_PER_ELEM + * is passed to the ImplicitKernelBase template to form the base class. + * + * Additionally, the number of degrees of freedom per support point for both + * the test and trial spaces are specified as `1` when specifying the base + * class. + */ +template +class ElectrostaticsKernel : public finiteElement::ImplicitKernelBase +{ +public: + using Base = finiteElement::ImplicitKernelBase; + + using Base::numDofPerTestSupportPoint; + using Base::numDofPerTrialSupportPoint; + using Base::maxNumTestSupportPointsPerElem; + using Base::m_dofNumber; + using Base::m_dofRankOffset; + using Base::m_matrix; + using Base::m_rhs; + using Base::m_elemsToNodes; + using Base::m_constitutiveUpdate; + using Base::m_finiteElementSpace; + using Base::m_meshData; + using Base::m_dt; + + /** + * @brief Constructor + * @copydoc geos::finiteElement::ImplicitKernelBase::ImplicitKernelBase + * @param fieldName The name of the primary field + */ + ElectrostaticsKernel(NodeManager const& nodeManager, EdgeManager const& edgeManager, + FaceManager const& faceManager, localIndex const targetRegionIndex, + SUBREGION_TYPE const& elementSubRegion, FE_TYPE const& finiteElementSpace, + CONSTITUTIVE_TYPE& inputConstitutiveType, arrayView1d const inputDofNumber, + globalIndex const rankOffset, CRSMatrixView const inputMatrix, + arrayView1d const inputRhs, real64 const inputDt, string const fieldName) + : + Base(nodeManager, edgeManager, faceManager, targetRegionIndex, elementSubRegion, finiteElementSpace, + inputConstitutiveType, inputDofNumber, rankOffset, inputMatrix, inputRhs, inputDt), + m_X(nodeManager.referencePosition()), + m_potential(nodeManager.template getReference>(fieldName)) + {} + + /** + * @class StackVariables + * @copydoc geos::finiteElement::ImplicitKernelBase::StackVariables + * + * Adds a stack array for the primary field. + */ + struct StackVariables : Base::StackVariables + { + public: + GEOS_HOST_DEVICE + StackVariables() + : Base::StackVariables(), xLocal(), primaryField_local{0.0} {} + +#if !defined(CALC_FEM_SHAPE_IN_KERNEL) + int xLocal; +#else + real64 xLocal[maxNumTestSupportPointsPerElem][3]; +#endif + + real64 primaryField_local[maxNumTestSupportPointsPerElem]; + }; + + /** + * @brief Copy global values from primary field to a local stack array. + * @copydoc geos::finiteElement::ImplicitKernelBase::setup + * + * For the ElectrostaticsKernel implementation, global values from the + * primaryField, and degree of freedom numbers are placed into element local + * stack storage. + */ + GEOS_HOST_DEVICE + inline void setup(localIndex const k, StackVariables& stack) const + { + m_finiteElementSpace.template setup(k, m_meshData, stack.feStack); + stack.numRows = m_finiteElementSpace.template numSupportPoints(stack.feStack); + stack.numCols = stack.numRows; + for (localIndex a = 0; a < stack.numRows; ++a) + { + localIndex const localNodeIndex = m_elemsToNodes(k, a); + +#if defined(CALC_FEM_SHAPE_IN_KERNEL) + for (int i = 0; i < 3; ++i) + { + stack.xLocal[a][i] = m_X[localNodeIndex][i]; + } +#endif + + stack.primaryField_local[a] = m_potential[localNodeIndex]; + stack.localRowDofIndex[a] = m_dofNumber[localNodeIndex]; + stack.localColDofIndex[a] = m_dofNumber[localNodeIndex]; + } + m_finiteElementSpace.template addGradGradStabilizationMatrix( + stack.feStack, stack.localJacobian); + } + + /** + * @copydoc geos::finiteElement::ImplicitKernelBase::quadraturePointKernel + */ + GEOS_HOST_DEVICE + inline void quadraturePointKernel(localIndex const k, localIndex const q, StackVariables& stack) const + { + real64 const conductivity = m_constitutiveUpdate.getConductivity(k); + real64 dNdX[maxNumTestSupportPointsPerElem][3]; + real64 const detJxW = m_finiteElementSpace.template getGradN(k, q, stack.xLocal, stack.feStack, dNdX); + + for (localIndex a = 0; a < stack.numRows; ++a) + { + for (localIndex b = 0; b < stack.numCols; ++b) + { + stack.localJacobian[a][b] += conductivity * LvArray::tensorOps::AiBi<3>(dNdX[a], dNdX[b]) * detJxW; + } + } + } + + /** + * @copydoc geos::finiteElement::ImplicitKernelBase::complete + * + * Form element residual from the fully formed element Jacobian dotted with + * the primary field and map the element local Jacobian/Residual to the + * global matrix/vector. + */ + GEOS_HOST_DEVICE + inline real64 complete(localIndex const k, StackVariables & stack) const + { + GEOS_UNUSED_VAR(k); + real64 maxForce = 0; + + for (localIndex a = 0; a < stack.numRows; ++a) + { + for (localIndex b = 0; b < stack.numRows; ++b) + { + stack.localResidual[a] += stack.localJacobian[a][b] * stack.primaryField_local[b]; + } + } + + for (localIndex a = 0; a < stack.numRows; ++a) + { + localIndex const dof = LvArray::integerConversion(stack.localRowDofIndex[a] - m_dofRankOffset); + if (dof < 0 || dof >= m_matrix.numRows()) continue; + m_matrix.template addToRowBinarySearchUnsorted( + dof, stack.localColDofIndex, stack.localJacobian[a], stack.numCols); + + RAJA::atomicAdd(&m_rhs[dof], stack.localResidual[a]); + maxForce = fmax(maxForce, fabs(stack.localResidual[a])); + } + + return maxForce; + } + +protected: + arrayView2d const m_X; + arrayView1d const m_potential; +}; + +using ElectrostaticsKernelFactory = finiteElement::KernelFactory const, + globalIndex const, CRSMatrixView const, + arrayView1d const, real64 const, string const>; +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_ELECTROSTATICS_ELECTROSTATICSKERNELS_HPP_ \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt index ced15c47ba5..4a8d25eeb37 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/solidMechanics/CMakeLists.txt @@ -36,6 +36,8 @@ set( solidMechanicsSolvers_headers kernels/ImplicitSmallStrainNewmark_impl.hpp kernels/ImplicitSmallStrainQuasiStatic.hpp kernels/ImplicitSmallStrainQuasiStatic_impl.hpp + kernels/ImplicitFiniteStrainQuasiStatic.hpp + kernels/ImplicitFiniteStrainQuasiStatic_impl.hpp SolidMechanicsStateReset.hpp SolidMechanicsStatistics.hpp contact/ContactSolverBase.hpp @@ -82,6 +84,7 @@ set( kernelTemplateFileList "" ) list( APPEND kernelTemplateFileList kernels/SolidMechanicsKernels.cpp.template + kernels/SolidMechanicsImplicitFiniteStrainKernels.cpp.template kernels/SolidMechanicsFixedStressThermoPoromechanicsKernels.cpp.template ) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 18e7714e97d..1d1d4b9e5cf 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -22,6 +22,7 @@ #include "SolidMechanicsLagrangianFEM.hpp" #include "kernels/ImplicitSmallStrainNewmark.hpp" #include "kernels/ImplicitSmallStrainQuasiStatic.hpp" +#include "kernels/ImplicitFiniteStrainQuasiStatic.hpp" #include "kernels/ExplicitSmallStrain.hpp" #include "kernels/ExplicitFiniteStrain.hpp" #include "kernels/FixedStressThermoPoromechanics.hpp" @@ -47,6 +48,7 @@ #include "physicsSolvers/LogLevelsInfo.hpp" #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsKernelsDispatchTypeList.hpp" +#include "physicsSolvers/solidMechanics/kernels/SolidMechanicsImplicitFiniteStrainKernelsDispatchTypeList.hpp" #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsFixedStressThermoPoromechanicsKernelsDispatchTypeList.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" @@ -1012,14 +1014,25 @@ void SolidMechanicsLagrangianFEM::setSparsityPattern( DomainPartition & domain, arrayView1d< globalIndex const > const dofNumber = nodeManager.getReference< globalIndex_array >( dofManager.getKey( solidMechanics::totalDisplacement::key() ) ); - finiteElement:: - fillSparsity< CellElementSubRegion, - solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, - regionNames, - getDiscretizationName(), - dofNumber, - dofManager.rankOffset(), - pattern ); + if (m_strainTheory == 0) { + finiteElement:: + fillSparsity< CellElementSubRegion, + solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh, + regionNames, + getDiscretizationName(), + dofNumber, + dofManager.rankOffset(), + pattern ); + } else if (m_strainTheory == 1) { + finiteElement:: + fillSparsity< CellElementSubRegion, + solidMechanicsLagrangianFEMKernels::ImplicitFiniteStrainQuasiStatic >( mesh, + regionNames, + getDiscretizationName(), + dofNumber, + dofManager.rankOffset(), + pattern ); + } } ); @@ -1140,14 +1153,26 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM { if( m_timeIntegrationOption == TimeIntegrationOption::QuasiStatic ) { - m_maxForce = assemblyLaunch< SolidMechanicsKernelsDispatchTypeList, - solidMechanicsLagrangianFEMKernels::QuasiStaticFactory >( mesh, - dofManager, - regionNames, - viewKeyStruct::solidMaterialNamesString(), - localMatrix, - localRhs, - dt ); + if (m_strainTheory == 0) { + m_maxForce = assemblyLaunch< SolidMechanicsKernelsDispatchTypeList, + solidMechanicsLagrangianFEMKernels::QuasiStaticFactory >( mesh, + dofManager, + regionNames, + viewKeyStruct::solidMaterialNamesString(), + localMatrix, + localRhs, + dt ); + } else if (m_strainTheory == 1) { + m_maxForce = + assemblyLaunch< SolidMechanicsImplicitFiniteStrainKernelsDispatchTypeList, + solidMechanicsLagrangianFEMKernels::ImplicitFiniteStrainQuasiStaticFactory >( mesh, + dofManager, + regionNames, + viewKeyStruct::solidMaterialNamesString(), + localMatrix, + localRhs, + dt ); + } } else if( m_timeIntegrationOption == TimeIntegrationOption::ImplicitDynamic ) { diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json index 9eb41224111..13edef0b6a8 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernelSpecs.json @@ -74,6 +74,30 @@ "CellElementSubRegion#ElasticIsotropic#H1_Prism11_VEM_Gauss1" ] }, + "SolidMechanicsImplicitFiniteStrainKernels": { + "vars": [ "SUBREGION_TYPE", "CONSTITUTIVE_TYPE", "FE_TYPE" ], + "constants": [ + [ "ImplicitFiniteStrainQuasiStaticPolicy", "geos::parallelDevicePolicy< GEOS_BLOCK_SIZE >" ] + ], + "combinations": { + "SUBREGION_TYPE": [ + "CellElementSubRegion" + ], + "CONSTITUTIVE_TYPE": [ + "ElasticIsotropicFiniteStrain" + ], + "FE_TYPE": [ + "H1_Hexahedron_Lagrange1_GaussLegendre2", + "H1_Wedge_Lagrange1_Gauss6", + "H1_Tetrahedron_Lagrange1_Gauss1", + "H1_Tetrahedron_Lagrange1_Gauss5", + "H1_Tetrahedron_Lagrange1_Gauss14", + "H1_Pyramid_Lagrange1_Gauss5" + ] + }, + "explicit": [] + }, + "SolidMechanicsFixedStressThermoPoromechanicsKernels": { "vars": [ "SUBREGION_TYPE", diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic.hpp new file mode 100644 index 00000000000..71b89af9ada --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic.hpp @@ -0,0 +1,239 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ImplicitFiniteStrainQuasistatic.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITFINITESTRAINQUASISTATIC_HPP_ +#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITFINITESTRAINQUASISTATIC_HPP_ + +#include "finiteElement/kernelInterface/ImplicitKernelBase.hpp" +#include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" + +namespace geos +{ + +namespace solidMechanicsLagrangianFEMKernels +{ + +/** + * @brief Implements kernels for solving finite strain quasi-static equilibrium. + * @copydoc geos::finiteElement::ImplicitKernelBase + * @tparam NUM_NODES_PER_ELEM The number of nodes per element for the + * @p SUBREGION_TYPE. + * @tparam UNUSED An unused parameter since we are assuming that the test and + * trial space have the same number of support points. + * + * ### QuasiStatic Description + * Implements the KernelBase interface functions required for solving the + * quasi-static equilibrium equations using one of the + * "finite element kernel application" functions such as + * geos::finiteElement::RegionBasedKernelApplication. + * + * In this implementation, the template parameter @p NUM_NODES_PER_ELEM is used + * in place of both @p NUM_TEST_SUPPORT_POINTS_PER_ELEM and + * @p NUM_TRIAL_SUPPORT_POINTS_PER_ELEM, which are assumed to be equal. This + * results in the @p UNUSED template parameter as only the NUM_NODES_PER_ELEM + * is passed to the ImplicitKernelBase template to form the base class. + * + * Additionally, the number of degrees of freedom per support point for both + * the test and trial spaces are specified as `3` when specifying the base + * class. + */ +template +class ImplicitFiniteStrainQuasiStatic : + public finiteElement::ImplicitKernelBase +{ +public: + /// Alias for the base class; + using Base = finiteElement::ImplicitKernelBase; + + /// Maximum number of nodes per element, which is equal to the maxNumTestSupportPointPerElem and + /// maxNumTrialSupportPointPerElem by definition. When the FE_TYPE is not a Virtual Element, this + /// will be the actual number of nodes per element. + static constexpr int numNodesPerElem = Base::maxNumTestSupportPointsPerElem; + using Base::numDofPerTestSupportPoint; + using Base::numDofPerTrialSupportPoint; + using Base::m_dofNumber; + using Base::m_dofRankOffset; + using Base::m_matrix; + using Base::m_rhs; + using Base::m_elemsToNodes; + using Base::m_constitutiveUpdate; + using Base::m_finiteElementSpace; + using Base::m_meshData; + using Base::m_dt; + + /** + * @brief Constructor + * @copydoc geos::finiteElement::ImplicitKernelBase::ImplicitKernelBase + * @param inputGravityVector The gravity vector. + */ + ImplicitFiniteStrainQuasiStatic(NodeManager const & nodeManager, EdgeManager const & edgeManager, + FaceManager const & faceManager, localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, arrayView1d const inputDofNumber, + globalIndex const rankOffset, CRSMatrixView const inputMatrix, + arrayView1d const inputRhs, real64 const inputDt, real64 const (&inputGravityVector)[3]); + + //***************************************************************************** + /** + * @class StackVariables + * @copydoc geos::finiteElement::ImplicitKernelBase::StackVariables + * + * Adds a stack array for the displacement, incremental displacement, and the + * constitutive stiffness. + */ + struct StackVariables : public Base::StackVariables + { + public: + + /// Constructor. + GEOS_HOST_DEVICE + StackVariables() + : Base::StackVariables(), + xLocal(), u_local(), uhat_local(), constitutiveStiffness() + {} + +#if !defined(CALC_FEM_SHAPE_IN_KERNEL) + /// Dummy + int xLocal; +#else + /// C-array stack storage for element local the nodal positions. + real64 xLocal[ numNodesPerElem ][ 3 ]; +#endif + + /// Stack storage for the element local nodal displacement + real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint]; + + /// Stack storage for the element local nodal incremental displacement + real64 uhat_local[numNodesPerElem][numDofPerTrialSupportPoint]; + + /// Stack storage for the constitutive stiffness at a quadrature point. + real64 constitutiveStiffness[ 9 ][ 9 ]; + }; + //***************************************************************************** + + /** + * @brief Copy global values from primary field to a local stack array. + * @copydoc ::geos::finiteElement::ImplicitKernelBase::setup + * + * For the QuasiStatic implementation, global values from the displacement, + * incremental displacement, and degree of freedom numbers are placed into + * element local stack storage. + */ + GEOS_HOST_DEVICE + void setup(localIndex const k, StackVariables & stack) const; + + /** + * @brief Internal struct to provide no-op defaults used in the inclusion + * of lambda functions into kernel component functions. + * @struct NoOpFunctors + */ + struct NoOpFunctors + { + /** + * @brief operator() no-op used for adding an additional dynamics term + * inside the jacobian assembly loop. + * @param a Node index for the row. + * @param b Node index for the col. + */ + GEOS_HOST_DEVICE inline constexpr + void operator() ( localIndex const a, localIndex const b ) + { + GEOS_UNUSED_VAR( a ); + GEOS_UNUSED_VAR( b ); + } + + /** + * @brief operator() no-op used for modifying the stress tensor prior to + * integrating the divergence to produce nodal forces. + * @param stress The stress array. + */ + GEOS_HOST_DEVICE inline constexpr + void operator() ( real64 (& stress)[3][3] ) + { + GEOS_UNUSED_VAR( stress ); + } + }; + + /** + * @copydoc geos::finiteElement::KernelBase::quadraturePointKernel + * @tparam STRESS_MODIFIER Type of optional functor to allow for the + * modification of stress prior to integration. + * @param stressModifier An optional functor to allow for the modification + * of stress prior to integration. + * For solid mechanics kernels, the strain increment is calculated, and the + * constitutive update is called. In addition, the constitutive stiffness + * stack variable is filled by the constitutive model. + */ + template< typename STRESS_MODIFIER = NoOpFunctors > + GEOS_HOST_DEVICE + void quadraturePointKernel(localIndex const k, localIndex const q, StackVariables & stack, + STRESS_MODIFIER && stressModifier = NoOpFunctors{}) const; + + /** + * @copydoc geos::finiteElement::ImplicitKernelBase::complete + */ + GEOS_HOST_DEVICE + inline + real64 complete( localIndex const k, + StackVariables & stack ) const; + + /** + * @copydoc geos::finiteElement::KernelBase::kernelLaunch + */ + template< typename POLICY, + typename KERNEL_TYPE > + static real64 + kernelLaunch( localIndex const numElems, + KERNEL_TYPE const & kernelComponent ); + +protected: + /// The array containing the nodal position array. + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X; + + /// The rank-global displacement array. + arrayView2d< real64 const, nodes::TOTAL_DISPLACEMENT_USD > const m_disp; + + /// The rank-global incremental displacement array. + arrayView2d< real64 const, nodes::INCR_DISPLACEMENT_USD > const m_uhat; + + /// The gravity vector. + real64 const m_gravityVector[3]; + + /// The rank global density + arrayView2d< real64 const > const m_density; + +}; + +/// The factory used to construct a ImplicitFiniteStrainQuasiStatic kernel. +using ImplicitFiniteStrainQuasiStaticFactory = + finiteElement::KernelFactory< ImplicitFiniteStrainQuasiStatic, + arrayView1d< globalIndex const > const, + globalIndex, + CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, + real64 const, + real64 const (&)[3] >; + +} // namespace solidMechanicsLagrangianFEMKernels + +} // namespace geos + +#include "finiteElement/kernelInterface/SparsityKernelBase.hpp" + +#endif // GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITFINITESTRAINQUASISTATIC_HPP_ \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic_impl.hpp new file mode 100644 index 00000000000..43f5dc66779 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic_impl.hpp @@ -0,0 +1,176 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file ImplicitFinitStrainQuasiStatic_impl.hpp + */ + +#ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITFINITESTRAINQUASISTATIC_IMPL_HPP_ +#define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITFINITESTRAINQUASISTATIC_IMPL_HPP_ + +#include "ImplicitFiniteStrainQuasiStatic.hpp" +namespace geos +{ + +namespace solidMechanicsLagrangianFEMKernels +{ + +template +ImplicitFiniteStrainQuasiStatic :: ImplicitFiniteStrainQuasiStatic( + NodeManager const & nodeManager, EdgeManager const & edgeManager, + FaceManager const & faceManager, localIndex const targetRegionIndex, + SUBREGION_TYPE const & elementSubRegion, FE_TYPE const & finiteElementSpace, + CONSTITUTIVE_TYPE & inputConstitutiveType, arrayView1d const inputDofNumber, + globalIndex const rankOffset, CRSMatrixView const inputMatrix, + arrayView1d const inputRhs, real64 const inputDt, real64 const (&inputGravityVector)[3]) + : Base(nodeManager, edgeManager, faceManager, targetRegionIndex, elementSubRegion, finiteElementSpace, + inputConstitutiveType, inputDofNumber, rankOffset, inputMatrix, inputRhs, inputDt), + m_X(nodeManager.referencePosition()), + m_disp(nodeManager.getField()), + m_uhat(nodeManager.getField()), + m_gravityVector{inputGravityVector[0], inputGravityVector[1], inputGravityVector[2]}, + m_density(inputConstitutiveType.getDensity()) + {} + +template +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ImplicitFiniteStrainQuasiStatic :: setup( + localIndex const k, StackVariables& stack) const +{ + m_finiteElementSpace.template setup(k, m_meshData, stack.feStack); + + localIndex const numSupportPoints = m_finiteElementSpace.template numSupportPoints(stack.feStack); + + stack.numRows = 3 * numSupportPoints; + stack.numCols = stack.numRows; + + for (localIndex a = 0; a < numSupportPoints; ++a) + { + localIndex const localNodeIndex = m_elemsToNodes(k, a); + + for (int i = 0; i < numDofPerTestSupportPoint; ++i) + { +#if defined(CALC_FEM_SHAPE_IN_KERNEL) + stack.xLocal[a][i] = m_X[localNodeIndex][i]; +#endif + stack.u_local[a][i] = m_disp[localNodeIndex][i]; + stack.uhat_local[a][i] = m_uhat[localNodeIndex][i]; + stack.localRowDofIndex[a * 3 + i] = m_dofNumber[localNodeIndex] + i; + stack.localColDofIndex[a * 3 + i] = m_dofNumber[localNodeIndex] + i; + } + } + + // Hanyu: not adding stabilization to the local jacobian since this is for FEM +} + +template +template +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ImplicitFiniteStrainQuasiStatic::quadraturePointKernel( + localIndex const k, localIndex const q, StackVariables & stack, STRESS_MODIFIER && stressModifier) const +{ + real64 dNdX[numNodesPerElem][3]; + real64 const detJxW = m_finiteElementSpace.template getGradN(k, q, stack.xLocal, stack.feStack, dNdX); + + real64 stress[3][3] = { {0} }; + typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness; + + real64 dUhatdX[3][3] = { {0} }; + real64 dUdX[3][3] = { {0} }; + real64 F[3][3] = { {0} }; + + FE_TYPE::gradient(dNdX, stack.u_local, dUdX); + FE_TYPE::gradient(dNdX, stack.uhat_local, dUhatdX); + + // calculate deformation gradient + LvArray::tensorOps::copy<3, 3>(F, dUdX); + LvArray::tensorOps::add<3, 3>(F, dUhatdX); + LvArray::tensorOps::addIdentity< 3 >(F, 1.0); + + m_constitutiveUpdate.finiteStrainUpdate(k, q, F, stress, stiffness); + + stressModifier(stress); + + for (localIndex i = 0; i < 3; ++i) { + for (localIndex j = 0; j < 3; ++j) { + stress[i][j] *= -detJxW; + } + } + + real64 const gravityForce[3] = { m_gravityVector[0] * m_density( k, q )* detJxW, + m_gravityVector[1] * m_density( k, q )* detJxW, + m_gravityVector[2] * m_density( k, q )* detJxW }; + + real64 N[numNodesPerElem]; + FE_TYPE::calcN( q, stack.feStack, N ); + FE_TYPE::plusGradNajAijPlusNaFi( dNdX, + stress, + N, + gravityForce, + reinterpret_cast< real64 (&)[numNodesPerElem][3] >(stack.localResidual) ); + + stiffness.template BTDB(dNdX, -detJxW, stack.localJacobian); +} + +template +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 ImplicitFiniteStrainQuasiStatic::complete( + localIndex const k, StackVariables & stack) const +{ + GEOS_UNUSED_VAR( k ); + real64 maxForce = 0; + + localIndex const numSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack ); + + // #pragma unroll + for( int localNode = 0; localNode < numSupportPoints; ++localNode ) + { + // #pragma unroll + for( int dim = 0; dim < numDofPerTestSupportPoint; ++dim ) + { + localIndex const dof = LvArray::integerConversion< localIndex >( stack.localRowDofIndex[ numDofPerTestSupportPoint * localNode + dim ] - m_dofRankOffset ); + if( dof < 0 || dof >= m_matrix.numRows() ) + continue; + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof, + stack.localRowDofIndex, + stack.localJacobian[ numDofPerTestSupportPoint * localNode + dim ], + stack.numRows ); + + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[ dof ], stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ); + maxForce = fmax( maxForce, fabs( stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ) ); + } + } + + return maxForce; +} + +template +template +GEOS_FORCE_INLINE +real64 +ImplicitFiniteStrainQuasiStatic::kernelLaunch(localIndex const numElems, + KERNEL_TYPE const & kernelComponent) +{ + return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent ); +} + +} // namespace solidMechanicsLagrangianFEMKernels + +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITFINITESTRAINQUASISTATIC_IMPL_HPP_ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsImplicitFiniteStrainKernels.cpp.template b/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsImplicitFiniteStrainKernels.cpp.template new file mode 100644 index 00000000000..231e328c807 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/SolidMechanicsImplicitFiniteStrainKernels.cpp.template @@ -0,0 +1,38 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic_impl.hpp" + +using ImplicitFiniteStrainQuasiStaticPolicy = @ImplicitFiniteStrainQuasiStaticPolicy@; + +#define INSTANTIATION( NAME )\ +template class NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ >; \ +template real64 NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ >::kernelLaunch< NAME##Policy, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ > > \ + ( localIndex const, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@ > const & ); \ + + +namespace geos +{ +using namespace constitutive; +using namespace finiteElement; +namespace solidMechanicsLagrangianFEMKernels +{ + INSTANTIATION( ImplicitFiniteStrainQuasiStatic ) +} +} + +#undef INSTANTIATION diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in b/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in index c14654cf0c7..cf8eb3fba33 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/policies.hpp.in @@ -21,6 +21,7 @@ using ExplicitFiniteStrainPolicy = @ExplicitFiniteStrainPolicy@; using FixedStressThermoPoromechanicsPolicy = @FixedStressThermoPoromechanicsPolicy@; using ImplicitSmallStrainNewmarkPolicy = @ImplicitSmallStrainNewmarkPolicy@; using ImplicitSmallStrainQuasiStaticPolicy = @ImplicitSmallStrainQuasiStaticPolicy@; +using ImplicitFiniteStrainQuasiStaticPolicy = @ImplicitFiniteStrainQuasiStaticPolicy@; #endif /* GEOS_CORECOMPONENTS_PHYSICSSOLVERSE_SOLIDMECHANICS_KERNELS_CONFIG_HPP */ diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 9c5a0188b78..9551fdaa471 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -371,6 +371,10 @@ + + + + @@ -655,6 +659,10 @@ + + + + @@ -783,6 +791,10 @@ + + + + @@ -795,6 +807,10 @@ + + + + @@ -2751,6 +2767,7 @@ Information output from lower logLevels is added with the desired log level + @@ -3908,6 +3925,57 @@ When set to `all` output both convergence & iteration information to a csv.--> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -6485,6 +6553,7 @@ Information output from lower logLevels is added with the desired log level + @@ -6517,9 +6586,11 @@ Information output from lower logLevels is added with the desired log level + + @@ -6706,6 +6777,12 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + + + + + @@ -7514,6 +7591,22 @@ For instance, if "oil" is before "gas" in "phaseNames", the table order should b + + + + + + + + + + + + + + + + @@ -7602,6 +7695,12 @@ For instance, if "oil" is before "gas" in "phaseNames", the table order should b + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 4a2f9c949c6..0b5824df857 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -572,6 +572,7 @@ A field can represent a physical variable. (pressure, temperature, global compos + @@ -983,6 +984,15 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + + + + + + @@ -1575,6 +1585,7 @@ A field can represent a physical variable. (pressure, temperature, global compos + @@ -1607,9 +1618,11 @@ A field can represent a physical variable. (pressure, temperature, global compos + + @@ -1825,6 +1838,10 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + @@ -2676,6 +2693,20 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + + + + + + + + + + + @@ -2742,6 +2773,10 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + +