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 90ae2cad970..045e1481a68 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_base.xml b/inputFiles/Electrostatics/Electrostatics_base.xml new file mode 100644 index 00000000000..e98455fb2b4 --- /dev/null +++ b/inputFiles/Electrostatics/Electrostatics_base.xml @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/Electrostatics/Electrostatics_variable_reactivity.xml b/inputFiles/Electrostatics/Electrostatics_variable_reactivity.xml new file mode 100644 index 00000000000..86c36193678 --- /dev/null +++ b/inputFiles/Electrostatics/Electrostatics_variable_reactivity.xml @@ -0,0 +1,173 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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..bcf13c576af --- /dev/null +++ b/inputFiles/Electrostatics/LiSEdefect_Electrostatics_main.xml @@ -0,0 +1,53 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 75075f8aae7..867dd2f1d3c 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< ElectroChemistryBase > +{ + template< typename LAMBDA > + static + void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler< ElectroChemistryBase >::execute( constitutiveRelation, + std::forward< LAMBDA >( lambda )); + } +}; + +/** + * @brief Material model for interface Butler-Volmer kinetics + */ +template<> +struct ConstitutivePassThru< ButlerVolmerInterface > +{ + template< typename LAMBDA > + static + void execute( ConstitutiveBase & constitutiveRelation, LAMBDA && lambda ) + { + ConstitutivePassThruHandler< ButlerVolmerInterface >::execute( constitutiveRelation, + std::forward< LAMBDA >( 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..d6c4142943e --- /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_reactivityCoeff() +{ + registerWrapper( viewKeyStruct::defaultReactivityCoefficientString(), &m_defaultKrxn ). + setApplyDefaultValue( 1.0 ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Default BV Reactivity Coefficient" ); + + registerWrapper( viewKeyStruct::reactivityCoefficientString(), &m_reactivityCoeff ). + setApplyDefaultValue( -1.0 ). // will be overwritten + setPlotLevel( PlotLevel::LEVEL_0 ). + setDescription( "Reactivity Coefficient Field" ); +} + +ButlerVolmerInterface::~ButlerVolmerInterface() {} + +void ButlerVolmerInterface::postInputInitialization() +{ + this->getWrapper< array1d< real64 > >( viewKeyStruct::reactivityCoefficientString()). + setApplyDefaultValue( m_defaultKrxn ); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, ButlerVolmerInterface, string const &, Group * const ) +} // namespace constitutive + +} // namespace geos diff --git a/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.hpp b/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.hpp new file mode 100644 index 00000000000..a3e890bbe17 --- /dev/null +++ b/src/coreComponents/constitutive/electroChemistry/ButlerVolmerReaction.hpp @@ -0,0 +1,127 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 reactivityCoeff the array of face-wise reactivity coefficient in the subregion + */ + ButlerVolmerReactionUpdate( arrayView1d< real64 > const & reactivityCoeff ) + : + m_reactivityCoeff( reactivityCoeff ) + {} + + /** + * @brief Number of elements storing reactivity coefficient data + * @return Number of elements + */ + localIndex numElem() const + { + return m_reactivityCoeff.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_reactivityCoeff[k]; + } + +protected: + + /// View on the face-wise reactivity coefficient + /// Also commonly referred to as k_rxn in electrochemistry + arrayView1d< real64 > m_reactivityCoeff; +}; + +/** + * @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_reactivityCoeff ); + } + +protected: + + /// Post-process XML input + virtual void postInputInitialization() override; + + array1d< real64 > m_reactivityCoeff; + + real64 m_defaultKrxn = 1.0; +}; + +} // namespace constitutive + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_BUTLERVOLMERREACTION_HPP diff --git a/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.cpp b/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.cpp new file mode 100644 index 00000000000..eb855d2fe4a --- /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< array1d< real64 > >( viewKeyStruct::conductivityString()). + setApplyDefaultValue( m_defaultConductivity ); +} + +REGISTER_CATALOG_ENTRY( ConstitutiveBase, ElectroChemistryBase, string const &, Group * const ) +} // namespace constitutive + +} // namespace geos diff --git a/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.hpp b/src/coreComponents/constitutive/electroChemistry/ElectroChemistryBase.hpp new file mode 100644 index 00000000000..fa7c058a902 --- /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< real64 > 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< real64 > 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< real64 > m_conductivity; + + real64 m_defaultConductivity = 1.0; +}; +} // namespace constitutive + +} // namespace geos + +#endif //GEOS_CONSTITUTIVE_ELECTROCHEMISTRYBASE_HPP diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.cpp b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.cpp new file mode 100644 index 00000000000..c494bafd0b6 --- /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 diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.hpp new file mode 100644 index 00000000000..b4674599629 --- /dev/null +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicFiniteStrain.hpp @@ -0,0 +1,448 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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; + + /// Compute the Hencky strain, namely ln(F), where F = I + grad(u) is the deformation gradient + 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: + /// Compute the forth-order material tangent tensor D_m, where \delta P = D_m * \delta F + 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..4401a7fc3fa --- /dev/null +++ b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullTensor.hpp @@ -0,0 +1,84 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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< NUM_SUPPORT_POINTS * 3, NUM_SUPPORT_POINTS * 3, 9 >( elmat, B, DB ); + LvArray::tensorOps::scale< NUM_SUPPORT_POINTS * 3, NUM_SUPPORT_POINTS * 3 >( elmat, detJxW ); + LvArray::tensorOps::add< NUM_SUPPORT_POINTS * 3, NUM_SUPPORT_POINTS * 3 >( 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..ef6a992f0b4 100644 --- a/src/coreComponents/constitutive/solid/SolidUtilities.hpp +++ b/src/coreComponents/constitutive/solid/SolidUtilities.hpp @@ -157,6 +157,149 @@ 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..bb05b3daa25 --- /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; +} 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..734ecc0e0fd --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.cpp @@ -0,0 +1,371 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" ) +{ + 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< PhysicsSolverBase >( m_surfaceGeneratorName ); +} + +void Electrostatics::registerDataOnMesh( Group & meshBodies ) +{ + forDiscretizationOnMeshTargets( meshBodies, + [&]( string const &, MeshLevel & mesh, string_array const & regionNames ) { + NodeManager & nodes = mesh.getNodeManager(); + + nodes.registerWrapper< real64_array >( 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< string >( viewKeyStruct::electroMaterialNamesString()). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setSizedFromParent( 0 ); + + string & electroMaterialName = subRegion.getReference< string >( viewKeyStruct::electroMaterialNamesString()); + electroMaterialName = PhysicsSolverBase::getConstitutiveName< ElectroChemistryBase >( 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< string >( viewKeyStruct::reactiveMaterialNamesString()). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setSizedFromParent( 0 ); + + string & reactiveMaterialName = subRegion.getReference< string >( viewKeyStruct::reactiveMaterialNamesString()); + reactiveMaterialName = PhysicsSolverBase::getConstitutiveName< ButlerVolmerInterface >( 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( 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< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > 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< globalIndex const > const & dofIndex = nodeManager.getReference< array1d< globalIndex > >( dofKey ); + + ElectrostaticsKernelFactory kernelFactory( dofIndex, dofManager.rankOffset(), localMatrix, localRhs, dt, m_fieldName ); + + finiteElement::regionBasedKernelApplication< parallelDevicePolicy<>, 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< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > 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< real64 const > 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< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) +{ + FieldSpecificationManager const & fsManager = FieldSpecificationManager::getInstance(); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), + [&]( string const &, MeshLevel & mesh, string_array const & ) + { + fsManager.apply< NodeManager >( time, mesh, m_fieldName, + [&]( FieldSpecificationBase const & bc, string const &, SortedArrayView< localIndex const > const & targetSet, + NodeManager & targetGroup, string const & GEOS_UNUSED_PARAM( fieldName )) + { + bc.applyBoundaryConditionToSystem< FieldSpecificationEqual, parallelDevicePolicy<> >( + 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< real64 > 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< globalIndex const > const blockLocalDofNumber = nodeManager.getReference< globalIndex_array >( dofKey ); + globalIndex const dofRankOffset = dofManager.rankOffset(); + + fsManager.template apply< FaceManager, TractionBoundaryCondition >( time, mesh, TractionBoundaryCondition::catalogName(), + [&]( TractionBoundaryCondition const & bc, string const &, SortedArrayView< localIndex const > const & targetSet, + Group &, string const & ) + { + bc.launch( time, blockLocalDofNumber, dofRankOffset, faceManager, nodeManager.referencePosition(), targetSet, localRhs ); + } ); + } ); +} + +void Electrostatics::applyButlerVolmerCurrent( DofManager const & dofManager, DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > 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< real64 const > const phi = nodeManager.getReference< array1d< real64 > >( m_fieldName ).toViewConst(); + + arrayView2d< real64 const > const faceNormal = faceManager.faceNormal(); + ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst(); + + string const dofKey = dofManager.getKey( m_fieldName ); + arrayView1d< globalIndex > const nodeDofNumber = nodeManager.getReference< globalIndex_array >( dofKey ); + globalIndex const rankOffset = dofManager.rankOffset(); + + constexpr localIndex maxNodesPerFace = 4; + constexpr localIndex maxDofPerElem = maxNodesPerFace * 2; + + elemManager.forElementSubRegions< FaceElementSubRegion >( [&]( FaceElementSubRegion & subRegion ) { + string const & constitutiveName = subRegion.getReference< string >( viewKeyStruct::reactiveMaterialNamesString()); + constitutive::ConstitutiveBase & constitutiveModel = subRegion.getConstitutiveModel( constitutiveName ); + constitutive::ButlerVolmerInterface & castedConstitutiveModel = dynamic_cast< constitutive::ButlerVolmerInterface & >(constitutiveModel); + constitutive::ButlerVolmerInterface::KernelWrapper const m_constitutiveUpdate( castedConstitutiveModel.createKernelUpdates()); + + arrayView1d< real64 > const area = subRegion.getElementArea(); + arrayView2d< localIndex const > const elemsToFaces = subRegion.faceList().toViewConst(); + + forAll< parallelDevicePolicy<> >( 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< localIndex >( rowDof[idof] - rankOffset ); + if( localRow >= 0 && localRow < localMatrix.numRows()) + { + localMatrix.addToRowBinarySearchUnsorted< parallelDeviceAtomic >( + localRow, rowDof.data(), dRdPhi[idof].dataIfContiguous(), numNodesPerFace*2 ); + RAJA::atomicAdd< parallelDeviceAtomic >( &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 ) +} diff --git a/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.hpp b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.hpp new file mode 100644 index 00000000000..097375f5304 --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/Electrostatics.hpp @@ -0,0 +1,138 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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< real64, globalIndex > & 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< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ) override; + + virtual void applyBoundaryConditions( real64 const time, real64 const dt, + DomainPartition & domain, DofManager const & dofManager, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > 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< real64 const > 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< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + void applyCurrentBC( real64 const time, DofManager const & dofManager, + DomainPartition & domain, arrayView1d< real64 > const & localRhs ); + + void applyButlerVolmerCurrent( DofManager const & dofManager, DomainPartition & domain, + CRSMatrixView< real64, globalIndex const > const & localMatrix, + arrayView1d< real64 > const & localRhs ); + + struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct + { + 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; + + 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; +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_ELECTROSTATICS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/electrostatics/ElectrostaticsKernels.hpp b/src/coreComponents/physicsSolvers/electrostatics/ElectrostaticsKernels.hpp new file mode 100644 index 00000000000..4f1dc26305a --- /dev/null +++ b/src/coreComponents/physicsSolvers/electrostatics/ElectrostaticsKernels.hpp @@ -0,0 +1,205 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +class ElectrostaticsKernel : public finiteElement::ImplicitKernelBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE, 1, 1 > +{ +public: + using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE, 1, 1 >; + + 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< globalIndex const > const inputDofNumber, + globalIndex const rankOffset, CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > 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< array1d< real64 > >( 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 } {} + + real64 xLocal[maxNumTestSupportPointsPerElem][3]; + + 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< FE_TYPE >( k, m_meshData, stack.feStack ); + stack.numRows = m_finiteElementSpace.getNumSupportPoints( 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< FE_TYPE, numDofPerTrialSupportPoint >( + 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 = FE_TYPE::calcGradN( 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< localIndex >( stack.localRowDofIndex[a] - m_dofRankOffset ); + if( dof < 0 || dof >= m_matrix.numRows()) continue; + m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( + dof, stack.localColDofIndex, stack.localJacobian[a], stack.numCols ); + + RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[dof], stack.localResidual[a] ); + maxForce = fmax( maxForce, fabs( stack.localResidual[a] )); + } + + return maxForce; + } + +protected: + arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const m_X; + arrayView1d< real64 const > const m_potential; +}; + +using ElectrostaticsKernelFactory = finiteElement::KernelFactory< ElectrostaticsKernel, arrayView1d< globalIndex const > const, + globalIndex const, CRSMatrixView< real64, globalIndex const > const, + arrayView1d< real64 > const, real64 const, string const >; +} // namespace geos + +#endif // GEOS_PHYSICSSOLVERS_ELECTROSTATICS_ELECTROSTATICSKERNELS_HPP_ 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 525be734d47..c051a7889af 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" @@ -1063,14 +1065,28 @@ 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 ); + } } ); @@ -1191,14 +1207,29 @@ 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..01ef9d79721 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic.hpp @@ -0,0 +1,234 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +class ImplicitFiniteStrainQuasiStatic : + public finiteElement::ImplicitKernelBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE, 3, 3 > +{ +public: + /// Alias for the base class; + using Base = finiteElement::ImplicitKernelBase< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE, 3, 3 >; + + /// 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< globalIndex const > const inputDofNumber, + globalIndex const rankOffset, CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > 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() + {} + + /// C-array stack storage for element local the nodal positions. + real64 xLocal[ numNodesPerElem ][ 3 ]; + + /// 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_ 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..4e8080e86d1 --- /dev/null +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitFiniteStrainQuasiStatic_impl.hpp @@ -0,0 +1,180 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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" +#include "finiteElement/elementFormulations/FiniteElementOperators.hpp" + +namespace geos +{ + +namespace solidMechanicsLagrangianFEMKernels +{ + +template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +ImplicitFiniteStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE > :: 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< globalIndex const > const inputDofNumber, + globalIndex const rankOffset, CRSMatrixView< real64, globalIndex const > const inputMatrix, + arrayView1d< real64 > 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< fields::solidMechanics::totalDisplacement >()), + m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >()), + m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] }, + m_density( inputConstitutiveType.getDensity()) +{} + +template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ImplicitFiniteStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE > :: setup( + localIndex const k, StackVariables & stack ) const +{ + m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack ); + + localIndex const numSupportPoints = m_finiteElementSpace.getNumSupportPoints( 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< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +template< typename STRESS_MODIFIER > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ImplicitFiniteStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::quadraturePointKernel( + localIndex const k, localIndex const q, StackVariables & stack, STRESS_MODIFIER && stressModifier ) const +{ + real64 dNdX[numNodesPerElem][3]; + real64 const detJxW = FE_TYPE::calcGradN( 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} }; + + finiteElement::feOps::gradient( dNdX, stack.u_local, dUdX ); + finiteElement::feOps::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 ); + finiteElement::feOps::plusGradNajAijPlusNaFi( dNdX, + stress, + N, + gravityForce, + reinterpret_cast< real64 (&)[numNodesPerElem][3] >(stack.localResidual) ); + + stiffness.template BTDB< numNodesPerElem >( dNdX, -detJxW, stack.localJacobian ); +} + +template< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +real64 ImplicitFiniteStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::complete( + localIndex const k, StackVariables & stack ) const +{ + GEOS_UNUSED_VAR( k ); + real64 maxForce = 0; + + localIndex const numSupportPoints = m_finiteElementSpace.getNumSupportPoints( 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< typename SUBREGION_TYPE, typename CONSTITUTIVE_TYPE, typename FE_TYPE > +template< typename POLICY, typename KERNEL_TYPE > +GEOS_FORCE_INLINE +real64 +ImplicitFiniteStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >::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..2c4c6979a46 --- /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@_impl >; \ +template real64 NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl >::kernelLaunch< NAME##Policy, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl > > \ + ( localIndex const, \ + NAME < @SUBREGION_TYPE@, @CONSTITUTIVE_TYPE@, @FE_TYPE@_impl> 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 728d7ec8eff..79c4d90c457 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -375,6 +375,10 @@ + + + + @@ -671,6 +675,10 @@ + + + + @@ -799,6 +807,10 @@ + + + + @@ -811,6 +823,10 @@ + + + + @@ -2768,6 +2784,7 @@ Information output from lower logLevels is added with the desired log level + @@ -3973,6 +3990,48 @@ Information output from lower logLevels is added with the desired log level + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -6635,6 +6694,7 @@ Information output from lower logLevels is added with the desired log level + @@ -6667,9 +6727,11 @@ Information output from lower logLevels is added with the desired log level + + @@ -6856,6 +6918,12 @@ The expected format is "{ waterMax, oilMax }", in that order--> + + + + + + @@ -7664,6 +7732,22 @@ For instance, if "oil" is before "gas" in "phaseNames", the table order should b + + + + + + + + + + + + + + + + @@ -7752,6 +7836,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 fc732d195f6..f16aad57777 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -573,6 +573,7 @@ A field can represent a physical variable. (pressure, temperature, global compos + @@ -996,6 +997,15 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + + + + + + @@ -1603,6 +1613,7 @@ A field can represent a physical variable. (pressure, temperature, global compos + @@ -1635,9 +1646,11 @@ A field can represent a physical variable. (pressure, temperature, global compos + + @@ -1853,6 +1866,10 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + @@ -2704,6 +2721,20 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + + + + + + + + + + + + @@ -2770,6 +2801,10 @@ A field can represent a physical variable. (pressure, temperature, global compos + + + +