diff --git a/framework/doc/content/source/kokkos/bcs/KokkosADCoupledVarNeumannBC.md b/framework/doc/content/source/kokkos/bcs/KokkosADCoupledVarNeumannBC.md new file mode 100644 index 000000000000..2090e43710f8 --- /dev/null +++ b/framework/doc/content/source/kokkos/bcs/KokkosADCoupledVarNeumannBC.md @@ -0,0 +1,20 @@ +# KokkosADCoupledVarNeumannBC + +!if! function=hasCapability('kokkos') + +This is the Kokkos version of [ADCoupledVarNeumannBC](CoupledVarNeumannBC.md). See the original document for details. + +## Example Input Syntax + +!listing test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann.i start=[right] end=[] include-end=true + +!syntax parameters /BCs/KokkosADCoupledVarNeumannBC + +!syntax inputs /BCs/KokkosADCoupledVarNeumannBC + +!syntax children /BCs/KokkosADCoupledVarNeumannBC + +!if-end! + +!else +!include kokkos/kokkos_warning.md diff --git a/framework/doc/content/source/kokkos/bcs/KokkosADDirichletBC.md b/framework/doc/content/source/kokkos/bcs/KokkosADDirichletBC.md new file mode 100644 index 000000000000..26949349b010 --- /dev/null +++ b/framework/doc/content/source/kokkos/bcs/KokkosADDirichletBC.md @@ -0,0 +1,20 @@ +# KokkosADDirichletBC + +!if! function=hasCapability('kokkos') + +This is the Kokkos version of [ADDirichletBC](ADDirichletBC.md). See the original document for details. + +## Example Input Syntax + +!listing test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_test.i start=[left] end=[] include-end=true + +!syntax parameters /BCs/KokkosADDirichletBC + +!syntax inputs /BCs/KokkosADDirichletBC + +!syntax children /BCs/KokkosADDirichletBC + +!if-end! + +!else +!include kokkos/kokkos_warning.md diff --git a/framework/doc/content/source/kokkos/bcs/KokkosADNeumannBC.md b/framework/doc/content/source/kokkos/bcs/KokkosADNeumannBC.md new file mode 100644 index 000000000000..78e63a312dbf --- /dev/null +++ b/framework/doc/content/source/kokkos/bcs/KokkosADNeumannBC.md @@ -0,0 +1,20 @@ +# KokkosADNeumannBC + +!if! function=hasCapability('kokkos') + +This is the Kokkos version of [ADNeumannBC](ADNeumannBC.md). See the original document for details. + +## Example Input Syntax + +!listing test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_neumannbc_test.i start=[right] end=[] include-end=true + +!syntax parameters /BCs/KokkosADNeumannBC + +!syntax inputs /BCs/KokkosADNeumannBC + +!syntax children /BCs/KokkosADNeumannBC + +!if-end! + +!else +!include kokkos/kokkos_warning.md diff --git a/framework/doc/content/source/kokkos/kernels/KokkosADCoupledTimeDerivative.md b/framework/doc/content/source/kokkos/kernels/KokkosADCoupledTimeDerivative.md new file mode 100644 index 000000000000..bbc3a84dbc41 --- /dev/null +++ b/framework/doc/content/source/kokkos/kernels/KokkosADCoupledTimeDerivative.md @@ -0,0 +1,20 @@ +# KokkosADCoupledTimeDerivative + +!if! function=hasCapability('kokkos') + +This is the Kokkos version of [ADCoupledTimeDerivative](ADCoupledTimeDerivative.md). See the original document for details. + +## Example Input Syntax + +!listing test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_coupled_time_derivative_test.i start=[time_v] end=[] include-end=true + +!syntax parameters /Kernels/KokkosADCoupledTimeDerivative + +!syntax inputs /Kernels/KokkosADCoupledTimeDerivative + +!syntax children /Kernels/KokkosADCoupledTimeDerivative + +!if-end! + +!else +!include kokkos/kokkos_warning.md diff --git a/framework/doc/content/source/kokkos/kernels/KokkosADDiffusion.md b/framework/doc/content/source/kokkos/kernels/KokkosADDiffusion.md new file mode 100644 index 000000000000..bd3d38e3dff2 --- /dev/null +++ b/framework/doc/content/source/kokkos/kernels/KokkosADDiffusion.md @@ -0,0 +1,20 @@ +# KokkosADDiffusion + +!if! function=hasCapability('kokkos') + +This is the Kokkos version of [ADDiffusion](ADDiffusion.md). See the original document for details. + +## Example Input Syntax + +!listing test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_test.i start=[diff] end=[] include-end=true + +!syntax parameters /Kernels/KokkosADDiffusion + +!syntax inputs /Kernels/KokkosADDiffusion + +!syntax children /Kernels/KokkosADDiffusion + +!if-end! + +!else +!include kokkos/kokkos_warning.md diff --git a/framework/doc/content/source/kokkos/kernels/KokkosADTimeDerivative.md b/framework/doc/content/source/kokkos/kernels/KokkosADTimeDerivative.md new file mode 100644 index 000000000000..90fddf7d4883 --- /dev/null +++ b/framework/doc/content/source/kokkos/kernels/KokkosADTimeDerivative.md @@ -0,0 +1,20 @@ +# KokkosADTimeDerivative + +!if! function=hasCapability('kokkos') + +This is the Kokkos version of [ADTimeDerivative](ADTimeDerivative.md). See the original document for details. + +## Example Input Syntax + +!listing test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_coupled_time_derivative_test.i start=[time_u] end=[] include-end=true + +!syntax parameters /Kernels/KokkosADTimeDerivative + +!syntax inputs /Kernels/KokkosADTimeDerivative + +!syntax children /Kernels/KokkosADTimeDerivative + +!if-end! + +!else +!include kokkos/kokkos_warning.md diff --git a/framework/doc/content/syntax/Kokkos/index.md b/framework/doc/content/syntax/Kokkos/index.md index aeb48cad6c6a..2f60c6d45834 100644 --- a/framework/doc/content/syntax/Kokkos/index.md +++ b/framework/doc/content/syntax/Kokkos/index.md @@ -317,7 +317,7 @@ For example, a postprocessor value which is always provided as read-only can be !alert note The data type should be copy-constructable. -### Static Polymorphism with Curiously Recurring Template Pattern (CRTP) id=kokkos_crtp +### Static Polymorphism id=kokkos_crtp The primary challenge in porting MOOSE to GPU lies in its heavy reliance on dynamic polymorphism using virtual functions. Polymorphism is the centerpiece of the +*template method pattern*+, which is a behavioral design pattern in object-oriented programming that establishes the skeleton of an algorithm in a base class while permitting derived classes to override certain steps without altering the algorithm’s overall structure, and it is the key design pattern of MOOSE. @@ -332,8 +332,9 @@ Furthermore, not every GPU backend supports virtual functions. Aside from portability, using virtual functions on GPU should be avoided if possible for performance, especially when the virtual functions are called in critical paths. The vtable lookup itself incurs overheads, and using function pointers prevents inlining. GPU compilers heavily rely on the inlining to generate an optimized code, and being unable to inline functions will likely lead to a performance hit. +Therefore, any polymorphism on GPU should be implemented statically. -Therefore, any polymorphism on GPU should be implemented statically, which can be achieved by the CRTP. +While not directly being used in Kokkos-MOOSE, the most famous static polymorphism design pattern is the Curiously Recurring Template Pattern (CRTP). The CRTP is a programming idiom that involves a class template inheriting from a template instantiation of itself, which is a technique used to achieve static (compile-time) polymorphism. The following pseudo-codes demonstrate a typical template method pattern implemented with the dynamic polymorphism and its equivalent implementation with the static polymorphism using the CRTP: @@ -405,9 +406,21 @@ Because the class type of the final derived class should be seen by the base cla If you accidentally derive a class from a non-template class, the base class will not be able to see the derived class. In this case, you are encouraged to prevent unintended inheritance by explicitly adding the `final` keyword to the last level derived class. -The Kokkos-MOOSE base classes are carefully designed to avoid the CRTP by leveraging a registry design pattern where external dispatchers are registered together with the objects. -Namely, the base classes themselves are not template classes, which alleviates the burden of users in dealing with class templates. -The hook methods provided by Kokkos-MOOSE are also template functions with respect to your object type, so you can directly cast `this` pointer in the hook methods without having to rely on class templates. +Kokkos-MOOSE is carefully designed to avoid the CRTP by leveraging a registry design pattern, where external dispatchers are registered together with the objects. +However, the basic principles of static polymorphism in the CRTP (templates, static casting of `this`, methods being public, function hiding) are still applicable and important to understand. +In Kokkos-MOOSE, the base classes are not template classes, which alleviates the burden of users in dealing with class templates. +Instead, the hook methods provided by Kokkos-MOOSE are template functions with respect to your object type, so you can directly cast `this` pointer in the hook methods without having to rely on class templates. +For example, the `computeQpResidual` hook method of Kokkos-MOOSE [Kernels](syntax/KokkosKernels/index.md) is defined as a template function: + +```cpp +template +KOKKOS_FUNCTION Real computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const; +``` + +The function template argument `Derived` replaces the class template argument in the CRTP and corresponds to your derived object type. +You can safely cast `this` pointer to the `Derived` type statically in your base class and directly call derived class methods, which are still required to be public methods. ### Separate Compilation diff --git a/framework/doc/content/syntax/KokkosBCs/index.md b/framework/doc/content/syntax/KokkosBCs/index.md index 1529b4f9955f..6e99d9c18b66 100644 --- a/framework/doc/content/syntax/KokkosBCs/index.md +++ b/framework/doc/content/syntax/KokkosBCs/index.md @@ -8,12 +8,11 @@ Before reading this documentation, consider reading the following materials firs - [Getting Started with Kokkos-MOOSE](syntax/Kokkos/index.md) to understand the programming practices for Kokkos-MOOSE, - [Kokkos Kernels System](syntax/KokkosKernels/index.md) to understand the common design pattern of objects in Kokkos-MOOSE. -!alert note -Kokkos-MOOSE boundary conditions do not support automatic differention yet. - The basic design pattern of Kokkos-MOOSE kernels described in [Kokkos Kernels System](syntax/Kokkos/index.md) applies to the boundary conditions as well. -You can create your own integrated and nodal boundary conditions by subclassing `Moose::Kokkos::IntegratedBC` and `Moose::Kokkos::NodalBC`, respectively, and following the same pattern with kernels including registering your boundary conditions with either `registerKokkosBoundaryCondition()` or `registerKokkosResidualObject()`. -Especially, integrated boundary conditions have identical interfaces with kernels, so they will not be explained here in detail. +You can create your own integrated and nodal boundary conditions by subclassing `Moose::Kokkos::IntegratedBC` and `Moose::Kokkos::NodalBC`, respectively, and following the same pattern with kernels including registering your boundary conditions with `registerKokkosResidualObject()`. +[Automatic differentiation (AD)](automatic_differentiation/index.md) versions of boundary conditions are also available and can be derived and registered in an analogous manner with the [AD kernels](syntax/KokkosKernels/index.md#kokkos_ad_kernel). + +Integrated boundary conditions have identical interfaces with kernels, so they will not be explained here in detail. See the following source codes of `KokkosCoupledVarNeumannBC` for an example of an integrated boundary condition: !listing framework/include/kokkos/bcs/KokkosCoupledVarNeumannBC.h id=kokkos-neumann-header @@ -42,7 +41,7 @@ To keep the consistency between interfaces, however, it is still passed as an ar `_current_node`, which is a pointer to the current libMesh node object, does not have a direct replacement. Instead, the node index can be queried by `datum.node()` and used to retrieve mesh data from the Kokkos mesh object. The node coordinate can also be obtained by `datum.q_point(qp)`. -As a result, the following residual function in `DirichletBCBase`: +As a result, the following residual function: ```cpp Real @@ -52,7 +51,7 @@ DirichletBCBase::computeQpResidual() } ``` -becomes the following in `Moose::Kokkos::DirichletBCBase`: +becomes the following: ```cpp template diff --git a/framework/doc/content/syntax/KokkosKernels/index.md b/framework/doc/content/syntax/KokkosKernels/index.md index 275d4ab8469c..9f8e9f4806f5 100644 --- a/framework/doc/content/syntax/KokkosKernels/index.md +++ b/framework/doc/content/syntax/KokkosKernels/index.md @@ -8,12 +8,12 @@ Before reading this documentation, consider reading the following materials firs - [Getting Started with Kokkos-MOOSE](syntax/Kokkos/index.md) to understand the programming practices for Kokkos-MOOSE. !alert note -Kokkos-MOOSE kernels do not support coupling with scalar variables and automatic differention yet. +Kokkos-MOOSE kernels do not support coupling with scalar variables yet. The Kokkos-MOOSE kernels are designed to resemble the original MOOSE kernels as much as possible for easier porting and adaptation. However, some differences still exist due to the fundamentally different programming paradigm between CPU and GPU. You can create your own kernel by subclassing `Moose::Kokkos::Kernel` as is done in the original MOOSE by inheriting `Kernel`. -However, your kernel should now be registered with either `registerKokkosKernel()` or `registerKokkosResidualObject()` instead of `registerMooseObject()`. +However, your kernel should now be registered with `registerKokkosResidualObject()` instead of `registerMooseObject()`. Also, the signatures of hook methods are different. In the original MOOSE, the following virtual functions should or optionally have been overriden: @@ -47,8 +47,7 @@ KOKKOS_FUNCTION Real computeQpOffDiagJacobian(const unsigned int i, AssemblyDatum & datum) const; ``` -The template argument `Derived` corresponds to your object type, and you can safely cast `this` pointer to the `Derived` type statically in your base class. -It can be useful for implementing polymorphism in a [CRTP-like](syntax/Kokkos/index.md#kokkos_crtp) fashion, as you can directly call the derived class methods using the cast pointer. +The template argument `Derived` can be used for implementing static polymorphism in a [CRTP-like](syntax/Kokkos/index.md#kokkos_crtp) fashion by statically casting `this` pointer to the derived type and directly calling the derived class methods using the cast pointer. Analogously to the original MOOSE, `computeQpResidual()` must be provided in the derived class, and the definition of `computeQpJacobian()` and `computeQpOffDiagJacobian()` are optional. The optional methods have default definitions in the base class, and redefining them in the derived class hides the base class definitions. @@ -99,13 +98,13 @@ KokkosDiffusion::computeQpJacobian(const unsigned int i, } ``` -See the following source codes of `KokkosCoupledForce` for another example of a kernel: +See the following source codes of `KokkosBodyForce` for another example of a kernel: -!listing framework/include/kokkos/kernels/KokkosCoupledForce.h id=kokkos-force-header - caption=The `KokkosCoupledForce` header file. +!listing framework/include/kokkos/kernels/KokkosBodyForce.h id=kokkos-force-header + caption=The `KokkosBodyForce` header file. -!listing framework/src/kokkos/kernels/KokkosCoupledForce.K id=kokkos-force-source language=cpp - caption=The `KokkosCoupledForce` source file. +!listing framework/src/kokkos/kernels/KokkosBodyForce.K id=kokkos-force-source language=cpp + caption=The `KokkosBodyForce` source file. !alert note [Every GPU function needs to be inlineable](syntax/Kokkos/index.md#kokkos_execution_space) and thus should be defined in headers. @@ -198,6 +197,21 @@ See the following source codes of `KokkosCoupledTimeDerivative` for an example o !listing framework/src/kokkos/kernels/KokkosCoupledTimeDerivative.K id=kokkos-time-derivative-source language=cpp caption=The `KokkosCoupledTimeDerivative` source file. +## Automatic Differentiation id=kokkos_ad_kernel + +Kokkos-MOOSE kernels also support [automatic differentiation (AD)](automatic_differentiation/index.md). +AD kernels can be derived from `Moose::Kokkos::ADKernel` and should be registered with `registerKokkosADResidualObject()`. +AD kernels require `computeQpResidual()` to be defined with the following signature, where everything remains the same with the ordinary kernels except the return type being `Moose::Kokkos::ADReal`: + +```cpp +template +KOKKOS_FUNCTION Moose::Kokkos::ADReal computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const; +``` + +`computeQpJacobian()` and `computeQpOffDiagJacobian()` are unused, as AD automatically assembles the Jacobian. + !syntax list /Kernels objects=True actions=False subsystems=False !if-end! diff --git a/framework/doc/content/syntax/KokkosNodalKernels/index.md b/framework/doc/content/syntax/KokkosNodalKernels/index.md index ff5cd3509072..45ec4ef48715 100644 --- a/framework/doc/content/syntax/KokkosNodalKernels/index.md +++ b/framework/doc/content/syntax/KokkosNodalKernels/index.md @@ -9,7 +9,7 @@ Before reading this documentation, consider reading the following materials firs - [Kokkos Kernels System](syntax/KokkosKernels/index.md) to understand the common design pattern of objects in Kokkos-MOOSE, - [Kokkos BCs System](syntax/KokkosBCs/index.md) to understand the design pattern of nodal boundary conditions in Kokkos-MOOSE. -You can create your own nodal kernels by inheriting `Moose::Kokkos::NodalKernel` and following the same pattern with kernels and boundary conditions including registering with either `registerKokkosNodalKernel()` or `registerKokkosResidualObject()`. +You can create your own nodal kernels by inheriting `Moose::Kokkos::NodalKernel` and following the same pattern with kernels and boundary conditions including registering with `registerKokkosResidualObject()`. The interfaces of nodal kernels are identical to the nodal boundary conditions described in [Kokkos BCs System](syntax/KokkosBCs/index.md), so they will not be explained here in detail. See the following source codes of `KokkosCoupledForceNodalKernel` for an example of a nodal kernel: diff --git a/framework/include/base/Assembly.h b/framework/include/base/Assembly.h index 0f49fce3a54a..c4a47bb8f7de 100644 --- a/framework/include/base/Assembly.h +++ b/framework/include/base/Assembly.h @@ -1825,11 +1825,7 @@ class Assembly void saveLocalADArray(std::vector & re, unsigned int i, unsigned int ntest, - const ADRealEigenVector & v) const - { - for (unsigned int j = 0; j < v.size(); ++j, i += ntest) - re[i] += v(j); - } + const ADRealEigenVector & v) const; /** * Helper function for assembling diagonal Jacobian contriubutions on local diff --git a/framework/include/interfaces/Coupleable.h b/framework/include/interfaces/Coupleable.h index 1bd15d7dc495..4498b9eaf068 100644 --- a/framework/include/interfaces/Coupleable.h +++ b/framework/include/interfaces/Coupleable.h @@ -1950,6 +1950,104 @@ class Coupleable unsigned int comp = 0) const; Moose::Kokkos::VariableValue kokkosCoupledNodalDots(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledVectorTagValueByName(const std::string & var_name, + const std::string & tag_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagValuesByName(const std::string & var_name, + const std::string & tag_name) const; + Moose::Kokkos::ADVariableGradient kokkosADCoupledVectorTagGradientByName( + const std::string & var_name, const std::string & tag_name, unsigned int comp = 0) const; + Moose::Kokkos::ADVariableGradient + kokkosADCoupledVectorTagGradientsByName(const std::string & var_name, + const std::string & tag_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledVectorTagNodalValueByName( + const std::string & var_name, const std::string & tag_name, unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagNodalValuesByName(const std::string & var_name, + const std::string & tag_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledVectorTagDofValueByName( + const std::string & var_name, const std::string & tag_name, unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagDofValuesByName(const std::string & var_name, + const std::string & tag_name) const; + + Moose::Kokkos::ADVariableValue kokkosADCoupledVectorTagValue(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagValues(const std::string & var_name, + const std::string & tag_param_name) const; + Moose::Kokkos::ADVariableGradient + kokkosADCoupledVectorTagGradient(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableGradient + kokkosADCoupledVectorTagGradients(const std::string & var_name, + const std::string & tag_param_name) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagNodalValue(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagNodalValues(const std::string & var_name, + const std::string & tag_param_name) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagDofValue(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledVectorTagDofValues(const std::string & var_name, + const std::string & tag_param_name) const; + + Moose::Kokkos::ADVariableValue kokkosADCoupledValue(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledValues(const std::string & var_name) const; + Moose::Kokkos::ADVariableGradient kokkosADCoupledGradient(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableGradient kokkosADCoupledGradients(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalValue(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalValues(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDofValue(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDofValues(const std::string & var_name) const; + + Moose::Kokkos::ADVariableValue kokkosADCoupledValueOld(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledValuesOld(const std::string & var_name) const; + Moose::Kokkos::ADVariableGradient kokkosADCoupledGradientOld(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableGradient kokkosADCoupledGradientsOld(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalValueOld(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalValuesOld(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDofValueOld(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDofValuesOld(const std::string & var_name) const; + + Moose::Kokkos::ADVariableValue kokkosADCoupledValueOlder(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledValuesOlder(const std::string & var_name) const; + Moose::Kokkos::ADVariableGradient kokkosADCoupledGradientOlder(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableGradient + kokkosADCoupledGradientsOlder(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalValueOlder(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue + kokkosADCoupledNodalValuesOlder(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDofValueOlder(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDofValuesOlder(const std::string & var_name) const; + + Moose::Kokkos::ADVariableValue kokkosADCoupledDot(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledDots(const std::string & var_name) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalDot(const std::string & var_name, + unsigned int comp = 0) const; + Moose::Kokkos::ADVariableValue kokkosADCoupledNodalDots(const std::string & var_name) const; + Moose::Kokkos::Scalar kokkosCoupledDotDu(const std::string & var_name, unsigned int comp = 0) const; diff --git a/framework/include/interfaces/TaggingInterface.h b/framework/include/interfaces/TaggingInterface.h index a7b4e88f446e..ff3df3d35d61 100644 --- a/framework/include/interfaces/TaggingInterface.h +++ b/framework/include/interfaces/TaggingInterface.h @@ -680,21 +680,3 @@ TaggingInterface::setResidual(SystemBase & sys, const SetResidualFunctor set_res if (sys.hasVector(tag_id)) set_residual_functor(sys.getVector(tag_id)); } - -inline void -TaggingInterface::addResiduals(Assembly & assembly, const ADResidualsPacket & packet) -{ - addResiduals(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor); -} - -inline void -TaggingInterface::addResidualsAndJacobian(Assembly & assembly, const ADResidualsPacket & packet) -{ - addResidualsAndJacobian(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor); -} - -inline void -TaggingInterface::addJacobian(Assembly & assembly, const ADResidualsPacket & packet) -{ - addJacobian(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor); -} diff --git a/framework/include/kokkos/base/KokkosDatum.h b/framework/include/kokkos/base/KokkosDatum.h index 6cb116d1d4ae..b4745c44b7ee 100644 --- a/framework/include/kokkos/base/KokkosDatum.h +++ b/framework/include/kokkos/base/KokkosDatum.h @@ -179,11 +179,6 @@ class Datum */ KOKKOS_FUNCTION Real3 normals(const unsigned int qp); - /** - * Deprecated, will be removed - */ - KOKKOS_FUNCTION void reinit() {} - /** * Set local parallelization option * @param local_thread_id The current local thread ID @@ -389,6 +384,7 @@ class AssemblyDatum : public Datum const unsigned int comp = 0) : Datum(elem, side, assembly, systems), _tag(ivar.tag()), + _sys(ivar.sys(comp)), _ivar(ivar.var(comp)), _jvar(jvar), _ife(systems[ivar.sys(comp)].getFETypeID(_ivar)), @@ -415,6 +411,7 @@ class AssemblyDatum : public Datum const unsigned int comp = 0) : Datum(node, assembly, systems), _tag(ivar.tag()), + _sys(ivar.sys(comp)), _ivar(ivar.var(comp)), _jvar(jvar), _ife(systems[ivar.sys(comp)].getFETypeID(_ivar)), @@ -437,6 +434,11 @@ class AssemblyDatum : public Datum * @returns The number of local DOFs */ KOKKOS_FUNCTION unsigned int n_jdofs() const { return _n_jdofs; } + /** + * Get the system number of variable + * @returns The system number of variable + */ + KOKKOS_FUNCTION unsigned int sys() const { return _sys; } /** * Get the variable number * @returns The variable number @@ -467,12 +469,26 @@ class AssemblyDatum : public Datum * @returns The variable FE type ID */ KOKKOS_FUNCTION unsigned int jfe() const { return _jfe; } + /** + * Set whether to compute derivatives for automatic differentiation (AD) + * @param flag Whether to compute derivatives + */ + KOKKOS_FUNCTION void do_derivatives(const bool flag) { _do_derivatives = flag; } + /** + * Get whether to compute derivatives for automatic differentiation (AD) + * @returns Whether to compute derivatives + */ + KOKKOS_FUNCTION bool do_derivatives() const { return _do_derivatives; } protected: /** * Solution tag ID */ const TagID _tag; + /** + * System number + */ + const unsigned int _sys; /** * Variable numbers */ @@ -485,6 +501,10 @@ class AssemblyDatum : public Datum * Number of local DOFs */ const unsigned int _n_idofs = 1, _n_jdofs = 1; + /** + * Whether to compute derivatives for automatic differentiation (AD) + */ + bool _do_derivatives = true; }; } // namespace Moose::Kokkos diff --git a/framework/include/kokkos/base/KokkosDispatcher.h b/framework/include/kokkos/base/KokkosDispatcher.h index a94db4bf2266..0ad71a1ed5d0 100644 --- a/framework/include/kokkos/base/KokkosDispatcher.h +++ b/framework/include/kokkos/base/KokkosDispatcher.h @@ -391,17 +391,28 @@ class DispatcherRegistry registerMooseObjectAliased(app, classname, alias); \ callRegisterKokkosResidualObjectFunction(classname, alias) -#define registerKokkosKernel(app, classname) registerKokkosResidualObject(app, classname) -#define registerKokkosKernelAliased(app, classname, alias) \ - registerKokkosResidualObjectAliased(app, classname, alias) +// AD Kernel, NodalKernel, BC -#define registerKokkosNodalKernel(app, classname) registerKokkosResidualObject(app, classname) -#define registerKokkosNodalKernelAliased(app, classname, alias) \ - registerKokkosResidualObjectAliased(app, classname, alias) +#define callRegisterKokkosADResidualObjectFunction(classname, objectname) \ + static char registerKokkosADResidualObject##classname() \ + { \ + using namespace Moose::Kokkos; \ + \ + DispatcherRegistry::addDispatcher(objectname); \ + \ + return 0; \ + } \ + \ + static char combineNames(kokkos_dispatcher_ad_residual_object_##classname, __COUNTER__) = \ + registerKokkosADResidualObject##classname() -#define registerKokkosBoundaryCondition(app, classname) registerKokkosResidualObject(app, classname) -#define registerKokkosBoundaryConditionAliased(app, classname, alias) \ - registerKokkosResidualObjectAliased(app, classname, alias) +#define registerKokkosADResidualObject(app, classname) \ + registerMooseObject(app, classname); \ + callRegisterKokkosADResidualObjectFunction(classname, #classname) + +#define registerKokkosADResidualObjectAliased(app, classname, alias) \ + registerMooseObjectAliased(app, classname, alias); \ + callRegisterKokkosADResidualObjectFunction(classname, alias) // Material diff --git a/framework/include/kokkos/base/KokkosJaggedArray.h b/framework/include/kokkos/base/KokkosJaggedArray.h index ba837117424e..48d8b30033d2 100644 --- a/framework/include/kokkos/base/KokkosJaggedArray.h +++ b/framework/include/kokkos/base/KokkosJaggedArray.h @@ -411,7 +411,10 @@ JaggedArray::finalize() _dims.copyToDevice(); _offsets.copyToDevice(); - _data.create(_offsets.last() + stride); + + // Pad an extra element at the end to avoid accessing the bound in the following operators when + // the last inner array has zero size + _data.create(_offsets.last() + stride + 1); _finalized = true; } diff --git a/framework/include/kokkos/base/KokkosResidualObject.h b/framework/include/kokkos/base/KokkosResidualObject.h index 26be9b30dae8..e90367032d67 100644 --- a/framework/include/kokkos/base/KokkosResidualObject.h +++ b/framework/include/kokkos/base/KokkosResidualObject.h @@ -64,7 +64,7 @@ class ResidualObject : public ::ResidualObject, { mooseError("computeOffDiagJacobian() is not used for Kokkos residual objects."); } - virtual void computeResidualAndJacobian() override final + virtual void computeResidualAndJacobian() override { computeResidual(); computeJacobian(); @@ -131,7 +131,7 @@ class ResidualObject : public ::ResidualObject, const unsigned int comp = 0) const; /** * Accumulate or set local nodal residual contribution to tagged vectors - * @param add The flag whether to add or set the local residual + * @param add Whether to add or set the local residual * @param local_re The local nodal residual contribution * @param node The contiguous node ID * @param comp The variable component @@ -155,9 +155,21 @@ class ResidualObject : public ::ResidualObject, const unsigned int j, const unsigned int jvar, const unsigned int comp = 0) const; + /** + * Accumulate local elemental Jacobian contribution to tagged matrices using automatic + * differentiation (AD) + * @param local_ke The local elemental Jacobian contribution + * @param datum The AssemblyDatum object of the current thread + * @param i The test function DOF index + * @param comp The variable component + */ + KOKKOS_FUNCTION void accumulateTaggedElementalMatrix(const DNDerivativeType & local_ke, + const AssemblyDatum & datum, + const unsigned int i, + const unsigned int comp = 0) const; /** * Accumulate or set local nodal Jacobian contribution to tagged matrices - * @param add The flag whether to add or set the local residual + * @param add Whether to add or set the local Jacobian * @param local_ke The local nodal Jacobian contribution * @param node The contiguous node ID * @param jvar The variable number for column @@ -168,6 +180,18 @@ class ResidualObject : public ::ResidualObject, const ContiguousNodeID node, const unsigned int jvar, const unsigned int comp = 0) const; + /** + * Accumulate or set local nodal Jacobian contribution to tagged matrices using automatic + * differentiation (AD) + * @param add Whether to add or set the local Jacobian + * @param local_ke The local elemental Jacobian contribution + * @param node The contiguous node ID + * @param comp The variable component + */ + KOKKOS_FUNCTION void accumulateTaggedNodalMatrix(const bool add, + const DNDerivativeType & local_ke, + const ContiguousNodeID node, + const unsigned int comp = 0) const; /** * The common loop structure template for computing elemental residual @@ -265,6 +289,29 @@ ResidualObject::accumulateTaggedElementalMatrix(const Real local_ke, } } +KOKKOS_FUNCTION inline void +ResidualObject::accumulateTaggedElementalMatrix(const DNDerivativeType & local_ke, + const AssemblyDatum & datum, + const unsigned int i, + const unsigned int comp) const +{ + auto & sys = kokkosSystem(_kokkos_var.sys(comp)); + auto row = sys.getElemLocalDofIndex(datum.elem().id, i, _kokkos_var.var(comp)); + + for (dof_id_type t = 0; t < _matrix_tags.size(); ++t) + { + auto tag = _matrix_tags[t]; + + if (sys.isMatrixTagActive(tag) && !sys.hasNodalBCMatrixTag(row, tag)) + for (unsigned int j = 0; j < local_ke.size(); ++j) + { + auto col = local_ke.raw_index(j); + + ::Kokkos::atomic_add(&sys.getMatrixValue(row, col, tag), local_ke.raw_at(j)); + } + } +} + KOKKOS_FUNCTION inline void ResidualObject::accumulateTaggedNodalMatrix(const bool add, const Real local_ke, @@ -298,6 +345,38 @@ ResidualObject::accumulateTaggedNodalMatrix(const bool add, } } +KOKKOS_FUNCTION inline void +ResidualObject::accumulateTaggedNodalMatrix(const bool add, + const DNDerivativeType & local_ke, + const ContiguousNodeID node, + const unsigned int comp) const +{ + auto & sys = kokkosSystem(_kokkos_var.sys(comp)); + auto row = sys.getNodeLocalDofIndex(node, 0, _kokkos_var.var(comp)); + + for (dof_id_type t = 0; t < _matrix_tags.size(); ++t) + { + auto tag = _matrix_tags[t]; + auto & matrix = sys.getMatrix(tag); + + if (sys.isMatrixTagActive(tag)) + { + if (!add) + matrix.zero(row); + + for (unsigned int j = 0; j < local_ke.size(); ++j) + { + auto col = local_ke.raw_index(j); + + if (add) + matrix(row, col) += local_ke.raw_at(j); + else + matrix(row, col) = local_ke.raw_at(j); + } + } + } +} + template KOKKOS_FUNCTION void ResidualObject::computeResidualInternal(AssemblyDatum & datum, function body) const diff --git a/framework/include/kokkos/base/KokkosTypes.h b/framework/include/kokkos/base/KokkosTypes.h index 0bd85ea266b8..2f28905028aa 100644 --- a/framework/include/kokkos/base/KokkosTypes.h +++ b/framework/include/kokkos/base/KokkosTypes.h @@ -13,6 +13,10 @@ #include "KokkosScalar.h" #include "KokkosJaggedArray.h" +#ifdef MOOSE_KOKKOS_SCOPE +#include "KokkosADReal.h" +#endif + #include "MooseError.h" #include "MooseUtils.h" @@ -21,32 +25,44 @@ namespace Moose::Kokkos { +template +struct Vector3; + +using Real3 = Vector3; +using ADReal3 = Vector3; + struct Real33; -struct Real3 +template +struct Vector3 { - Real v[3]; + T v[3]; #ifdef MOOSE_KOKKOS_SCOPE - Real3(const libMesh::TypeVector & vector); - KOKKOS_INLINE_FUNCTION Real3() { *this = 0; } - KOKKOS_INLINE_FUNCTION Real3(const Real scalar) { *this = scalar; } - KOKKOS_INLINE_FUNCTION Real3(const Real3 & vector) { *this = vector; } - KOKKOS_INLINE_FUNCTION Real3(const Real x, const Real y, const Real z); - - KOKKOS_INLINE_FUNCTION Real3 operator-() const; - KOKKOS_INLINE_FUNCTION Real & operator()(unsigned int i) { return v[i]; } - KOKKOS_INLINE_FUNCTION Real operator()(unsigned int i) const { return v[i]; } - - Real3 & operator=(const libMesh::TypeVector & vector); - KOKKOS_INLINE_FUNCTION Real3 & operator=(const Real3 & vector); - KOKKOS_INLINE_FUNCTION Real3 & operator=(const Real scalar); - - KOKKOS_INLINE_FUNCTION void operator+=(const Real scalar); - KOKKOS_INLINE_FUNCTION void operator+=(const Real3 vector); - KOKKOS_INLINE_FUNCTION void operator-=(const Real scalar); - KOKKOS_INLINE_FUNCTION void operator-=(const Real3 vector); - KOKKOS_INLINE_FUNCTION void operator*=(const Real scalar); + Vector3(const libMesh::TypeVector & vector); + KOKKOS_INLINE_FUNCTION Vector3() { *this = T{}; } + KOKKOS_INLINE_FUNCTION Vector3(const T & scalar) { *this = scalar; } + KOKKOS_INLINE_FUNCTION Vector3(const Vector3 & vector) { *this = vector; } + KOKKOS_INLINE_FUNCTION Vector3(const T & x, const T & y, const T & z); + + KOKKOS_INLINE_FUNCTION Vector3 operator-() const; + KOKKOS_INLINE_FUNCTION T & operator()(unsigned int i) { return v[i]; } + KOKKOS_INLINE_FUNCTION const T & operator()(unsigned int i) const { return v[i]; } + + Vector3 & operator=(const libMesh::TypeVector & vector); + + template + KOKKOS_INLINE_FUNCTION Vector3 & operator=(const Vector3 & vector); + KOKKOS_INLINE_FUNCTION Vector3 & operator=(const Vector3 & vector); + KOKKOS_INLINE_FUNCTION Vector3 & operator=(const T & scalar); + + template + KOKKOS_INLINE_FUNCTION void operator+=(const Vector3 & vector); + KOKKOS_INLINE_FUNCTION void operator+=(const T & scalar); + template + KOKKOS_INLINE_FUNCTION void operator-=(const Vector3 & vector); + KOKKOS_INLINE_FUNCTION void operator-=(const T & scalar); + KOKKOS_INLINE_FUNCTION void operator*=(const T & scalar); KOKKOS_INLINE_FUNCTION Real norm() const; KOKKOS_INLINE_FUNCTION Real dot_product(const Real3 vector) const; @@ -82,42 +98,59 @@ struct Real33 #ifdef MOOSE_KOKKOS_SCOPE -inline Real3::Real3(const libMesh::TypeVector & vector) +template +Vector3::Vector3(const libMesh::TypeVector & vector) { v[0] = vector(0); v[1] = vector(1); v[2] = vector(2); } +template KOKKOS_INLINE_FUNCTION -Real3::Real3(const Real x, const Real y, const Real z) +Vector3::Vector3(const T & x, const T & y, const T & z) { v[0] = x; v[1] = y; v[2] = z; } -KOKKOS_INLINE_FUNCTION Real3 -Real3::operator-() const +template +Vector3 & +Vector3::operator=(const libMesh::TypeVector & vector) { - Real3 vector(*this); + v[0] = vector(0); + v[1] = vector(1); + v[2] = vector(2); + + return *this; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +Vector3::operator-() const +{ + Vector3 vector(*this); vector *= -1; return vector; } -inline Real3 & -Real3::operator=(const libMesh::TypeVector & vector) +template +KOKKOS_INLINE_FUNCTION Vector3 & +Vector3::operator=(const Vector3 & vector) { - v[0] = vector(0); - v[1] = vector(1); - v[2] = vector(2); + v[0] = vector.v[0]; + v[1] = vector.v[1]; + v[2] = vector.v[2]; return *this; } -KOKKOS_INLINE_FUNCTION Real3 & -Real3::operator=(const Real3 & vector) +template +template +KOKKOS_INLINE_FUNCTION Vector3 & +Vector3::operator=(const Vector3 & vector) { v[0] = vector.v[0]; v[1] = vector.v[1]; @@ -126,8 +159,9 @@ Real3::operator=(const Real3 & vector) return *this; } -KOKKOS_INLINE_FUNCTION Real3 & -Real3::operator=(const Real scalar) +template +KOKKOS_INLINE_FUNCTION Vector3 & +Vector3::operator=(const T & scalar) { v[0] = scalar; v[1] = scalar; @@ -136,60 +170,133 @@ Real3::operator=(const Real scalar) return *this; } +template +template KOKKOS_INLINE_FUNCTION void -Real3::operator+=(const Real scalar) -{ - v[0] += scalar; - v[1] += scalar; - v[2] += scalar; -} - -KOKKOS_INLINE_FUNCTION void -Real3::operator+=(const Real3 vector) +Vector3::operator+=(const Vector3 & vector) { v[0] += vector.v[0]; v[1] += vector.v[1]; v[2] += vector.v[2]; } +template KOKKOS_INLINE_FUNCTION void -Real3::operator-=(const Real scalar) +Vector3::operator+=(const T & scalar) { - v[0] -= scalar; - v[1] -= scalar; - v[2] -= scalar; + v[0] += scalar; + v[1] += scalar; + v[2] += scalar; } +template +template KOKKOS_INLINE_FUNCTION void -Real3::operator-=(const Real3 vector) +Vector3::operator-=(const Vector3 & vector) { v[0] -= vector.v[0]; v[1] -= vector.v[1]; v[2] -= vector.v[2]; } +template KOKKOS_INLINE_FUNCTION void -Real3::operator*=(const Real scalar) +Vector3::operator-=(const T & scalar) +{ + v[0] -= scalar; + v[1] -= scalar; + v[2] -= scalar; +} + +template +KOKKOS_INLINE_FUNCTION void +Vector3::operator*=(const T & scalar) { v[0] *= scalar; v[1] *= scalar; v[2] *= scalar; } +template +KOKKOS_INLINE_FUNCTION Vector3 +operator+(const T & left, const Vector3 & right) +{ + return {left + right.v[0], left + right.v[1], left + right.v[2]}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator+(const Vector3 & left, const T & right) +{ + return {left.v[0] + right, left.v[1] + right, left.v[2] + right}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator+(const Vector3 & left, const Vector3 & right) +{ + return {left.v[0] + right.v[0], left.v[1] + right.v[1], left.v[2] + right.v[2]}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator-(const T & left, const Vector3 & right) +{ + return {left - right.v[0], left - right.v[1], left - right.v[2]}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator-(const Vector3 & left, const T & right) +{ + return {left.v[0] - right, left.v[1] - right, left.v[2] - right}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator-(const Vector3 & left, const Vector3 & right) +{ + return {left.v[0] - right.v[0], left.v[1] - right.v[1], left.v[2] - right.v[2]}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator*(const T & left, const Vector3 & right) +{ + return {left * right.v[0], left * right.v[1], left * right.v[2]}; +} + +template +KOKKOS_INLINE_FUNCTION Vector3 +operator*(const Vector3 & left, const T & right) +{ + return {left.v[0] * right, left.v[1] * right, left.v[2] * right}; +} + +template +KOKKOS_INLINE_FUNCTION T +operator*(const Vector3 & left, const Vector3 & right) +{ + return left.v[0] * right.v[0] + left.v[1] * right.v[1] + left.v[2] * right.v[2]; +} + +template <> KOKKOS_INLINE_FUNCTION Real -Real3::norm() const +Vector3::norm() const { return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); } +template <> KOKKOS_INLINE_FUNCTION Real -Real3::dot_product(const Real3 vector) const +Vector3::dot_product(const Real3 vector) const { return v[0] * vector.v[0] + v[1] * vector.v[1] + v[2] * vector.v[2]; } +template <> KOKKOS_INLINE_FUNCTION Real3 -Real3::cross_product(const Real3 vector) const +Vector3::cross_product(const Real3 vector) const { Real3 cross; @@ -200,8 +307,9 @@ Real3::cross_product(const Real3 vector) const return cross; } +template <> KOKKOS_INLINE_FUNCTION Real33 -Real3::cartesian_product(const Real3 vector) const +Vector3::cartesian_product(const Real3 vector) const { Real33 tensor; @@ -325,24 +433,6 @@ Real33::col(const unsigned int j) const return Real3(a[0][j], a[1][j], a[2][j]); } -KOKKOS_INLINE_FUNCTION Real3 -operator*(const Real left, const Real3 right) -{ - return {left * right.v[0], left * right.v[1], left * right.v[2]}; -} - -KOKKOS_INLINE_FUNCTION Real3 -operator*(const Real3 left, const Real right) -{ - return {left.v[0] * right, left.v[1] * right, left.v[2] * right}; -} - -KOKKOS_INLINE_FUNCTION Real -operator*(const Real3 left, const Real3 right) -{ - return left.v[0] * right.v[0] + left.v[1] * right.v[1] + left.v[2] * right.v[2]; -} - KOKKOS_INLINE_FUNCTION Real3 operator*(const Real33 left, const Real3 right) { @@ -376,12 +466,6 @@ operator+(const Real3 left, const Real right) return {left.v[0] + right, left.v[1] + right, left.v[2] + right}; } -KOKKOS_INLINE_FUNCTION Real3 -operator+(const Real3 left, const Real3 right) -{ - return {left.v[0] + right.v[0], left.v[1] + right.v[1], left.v[2] + right.v[2]}; -} - KOKKOS_INLINE_FUNCTION Real3 operator-(const Real left, const Real3 right) { @@ -395,9 +479,45 @@ operator-(const Real3 left, const Real right) } KOKKOS_INLINE_FUNCTION Real3 -operator-(const Real3 left, const Real3 right) +operator*(const Real left, const Real3 right) { - return {left.v[0] - right.v[0], left.v[1] - right.v[1], left.v[2] - right.v[2]}; + return {left * right.v[0], left * right.v[1], left * right.v[2]}; +} + +KOKKOS_INLINE_FUNCTION Real3 +operator*(const Real3 left, const Real right) +{ + return {left.v[0] * right, left.v[1] * right, left.v[2] * right}; +} + +template ::type, ADReal>::value>::type> +KOKKOS_INLINE_FUNCTION ADReal3 +operator*(const Real3 left, const T & right) +{ + return {left(0) * right, left(1) * right, left(2) * right}; +} + +template ::type, ADReal>::value>::type> +KOKKOS_INLINE_FUNCTION ADReal3 +operator*(const T & left, const Real3 right) +{ + return {left * right(0), left * right(1), left * right(2)}; +} + +KOKKOS_INLINE_FUNCTION ADReal +operator*(const Real3 left, const ADReal3 & right) +{ + return left(0) * right(0) + left(1) * right(1) + left(2) * right(2); +} + +KOKKOS_INLINE_FUNCTION ADReal +operator*(const ADReal3 & left, const Real3 right) +{ + return left(0) * right(0) + left(1) * right(1) + left(2) * right(2); } #endif diff --git a/framework/include/kokkos/base/KokkosVariable.h b/framework/include/kokkos/base/KokkosVariable.h index eddf6b23a278..f630e2d90699 100644 --- a/framework/include/kokkos/base/KokkosVariable.h +++ b/framework/include/kokkos/base/KokkosVariable.h @@ -38,42 +38,100 @@ class Variable * @param variable The MOOSE variable * @param tag The vector tag ID */ - Variable(const MooseVariableBase & variable, const TagID tag) { init(variable, tag); } + Variable(const MooseVariableFieldBase & variable, const TagID tag) { init(variable, tag); } /** * Constructor * Initialize the variable with a MOOSE variable and vector tag name * @param variable The MOOSE variable * @param tag_name The vector tag name */ - Variable(const MooseVariableBase & variable, const TagName & tag_name = Moose::SOLUTION_TAG) + Variable(const MooseVariableFieldBase & variable, const TagName & tag_name = Moose::SOLUTION_TAG) { init(variable, tag_name); } + /** + * Constructor + * Initialize the variable with multiple MOOSE variables and vector tag ID + * @param variables The MOOSE variables + * @param tag The vector tag ID + */ + ///@{ + Variable(const std::vector & variables, const TagID tag) + { + init(variables, tag); + } + Variable(const std::vector & variables, const TagID tag) + { + init(variables, tag); + } + ///@} + /** + * Constructor + * Initialize the variable with multiple MOOSE variables and vector tag name + * @param variables The MOOSE variables + * @param tag The vector tag ID + */ + ///@{ + Variable(const std::vector & variables, + const TagName & tag_name = Moose::SOLUTION_TAG) + { + init(variables, tag_name); + } + Variable(const std::vector & variables, + const TagName & tag_name = Moose::SOLUTION_TAG) + { + init(variables, tag_name); + } + ///@} /** * Initialize the variable with a MOOSE variable and vector tag ID * @param variable The MOOSE variable * @param tag The vector tag ID */ - void init(const MooseVariableBase & variable, const TagID tag); + void init(const MooseVariableFieldBase & variable, const TagID tag); /** * Initialize the variable with a MOOSE variable and vector tag name * @param variable The MOOSE variable * @param tag_name The vector tag name */ - void init(const MooseVariableBase & variable, const TagName & tag_name = Moose::SOLUTION_TAG); + void init(const MooseVariableFieldBase & variable, + const TagName & tag_name = Moose::SOLUTION_TAG); /** - * Initialize the variable with coupled MOOSE variables - * @param variables The coupled MOOSE variables + * Initialize the variable with multiple MOOSE variables and vector tag ID + * @param variables The MOOSE variables * @param tag The vector tag ID */ - void - init(const std::vector & variables, const TagID tag, CoupleableKey); + ///@{ + void init(const std::vector & variables, const TagID tag); + void init(const std::vector & variables, const TagID tag); + ///@} + /** + * Initialize the variable with multiple MOOSE variables and vector tag name + * @param variables The MOOSE variables + * @param tag_name The vector tag name + */ + ///@{ + void init(const std::vector & variables, + const TagName & tag_name = Moose::SOLUTION_TAG); + void init(const std::vector & variables, + const TagName & tag_name = Moose::SOLUTION_TAG); + ///@} /** * Initialize the variable with coupled default values * @param values The default coupled values */ void init(const std::vector & values, CoupleableKey); + /** + * Get the MOOSE variable of a component + * @param comp The variable component + * @returns The MOOSE variable + */ + const MooseVariableFieldBase * mooseVar(unsigned int comp = 0) + { + return _moose_var.size() ? _moose_var[comp] : nullptr; + } + /** * Get whether the variable is initialized * @returns Whether the variable is initialized @@ -89,6 +147,16 @@ class Variable * @returns Whether the variable is nodal */ KOKKOS_FUNCTION bool nodal() const { return _nodal; } + /** + * Get whether the tag is time derivative + * @returns Whether the tag is time derivative + */ + KOKKOS_FUNCTION bool dot() const { return _dot; } + /** + * Get whether the tag is old/older value + * @returns Whether the tag is old/older value + */ + KOKKOS_FUNCTION bool old() const { return _old; } /** * Get the number of components * @returns The number of components @@ -136,6 +204,14 @@ class Variable * Whether the variable is nodal */ bool _nodal = false; + /** + * Whether the tag is time derivative + */ + bool _dot = false; + /** + * Whether the tag is old/older value + */ + bool _old = false; /** * Number of components */ @@ -144,6 +220,10 @@ class Variable * Vector tag ID */ TagID _tag = Moose::INVALID_TAG_ID; + /** + * MOOSE variable of each component + */ + Array _moose_var; /** * Variable number of each component */ diff --git a/framework/include/kokkos/base/KokkosVariableValue.h b/framework/include/kokkos/base/KokkosVariableValue.h index 0cadb3e76c0d..7b6bf288d8bd 100644 --- a/framework/include/kokkos/base/KokkosVariableValue.h +++ b/framework/include/kokkos/base/KokkosVariableValue.h @@ -11,7 +11,7 @@ #include "KokkosDatum.h" -#include "MooseVariableBase.h" +#include "MooseVariableFieldBase.h" namespace Moose::Kokkos { @@ -109,37 +109,75 @@ class VariableTestGradient : datum.assembly().getGradPhiFace(elem.subdomain, elem.type, fe)(side)(i, qp)); } }; + +using ADVariablePhiValue = VariablePhiValue; +using ADVariablePhiGradient = VariablePhiGradient; +using ADVariableTestValue = VariableTestValue; +using ADVariableTestGradient = VariableTestGradient; + ///@} /** * The Kokkos wrapper classes for MOOSE-like variable value access */ ///@{ -class VariableValue +template +class VariableValueTempl { + using real_type = std::conditional_t; + public: /** * Default constructor */ - VariableValue() = default; + VariableValueTempl() = default; /** * Constructor * @param var The Kokkos variable * @param dof Whether to get DOF values */ - VariableValue(Variable var, bool dof = false) : _var(var), _dof(dof) {} + VariableValueTempl(Variable var, bool dof = false) : _var(var), _dof(dof) {} /** * Constructor * @param var The MOOSE variable * @param tag The vector tag name * @param dof Whether to get DOF values */ - VariableValue(const MooseVariableBase & var, - const TagName & tag = Moose::SOLUTION_TAG, - bool dof = false) + VariableValueTempl(const MooseVariableFieldBase & var, + const TagName & tag = Moose::SOLUTION_TAG, + bool dof = false) : _var(var, tag), _dof(dof) { } + /** + * Constructor + * @param vars The MOOSE variables + * @param tag The vector tag name + * @param dof Whether to get DOF values + */ + ///@{ + VariableValueTempl(const std::vector & vars, + const TagName & tag = Moose::SOLUTION_TAG, + bool dof = false) + : _var(vars, tag), _dof(dof) + { + } + VariableValueTempl(const std::vector & vars, + const TagName & tag = Moose::SOLUTION_TAG, + bool dof = false) + : _var(vars, tag), _dof(dof) + { + } + ///@} + + /** + * Copy constructor for parallel dispatch + */ + VariableValueTempl(const VariableValueTempl & object); + /** + * Copy assignment operator + */ + VariableValueTempl & operator=(const VariableValueTempl & object); /** * Get whether the variable was coupled @@ -150,11 +188,24 @@ class VariableValue /** * Get the current variable value * @param datum The Datum object of the current thread - * @param qp The local quadrature point index + * @param idx The local quadrature point or DOF index * @param comp The variable component * @returns The variable value */ - KOKKOS_FUNCTION Real operator()(Datum & datum, unsigned int qp, unsigned int comp = 0) const; + KOKKOS_FUNCTION auto operator()(Datum & datum, unsigned int idx, unsigned int comp = 0) const + { + return get(datum, idx, comp); + } + + /** + * Get the current variable value + * @param datum The AssemblyDatum object of the current thread + * @param idx The local quadrature point or DOF index + * @param comp The variable component + * @returns The variable value + */ + KOKKOS_FUNCTION auto + operator()(AssemblyDatum & datum, unsigned int idx, unsigned int comp = 0) const; /** * Get the Kokkos variable @@ -163,37 +214,185 @@ class VariableValue KOKKOS_FUNCTION const Variable & variable() const { return _var; } private: + /** + * Get the current variable value + * @param datum The Datum object of the current thread + * @param idx The local quadrature point or DOF index + * @param comp The variable component + * @param seed The derivative seed (only meaningful for AD) + * @returns The variable value + */ + KOKKOS_FUNCTION auto + get(Datum & datum, unsigned int idx, unsigned int comp = 0, Real seed = 0) const; + /** * Coupled Kokkos variable */ Variable _var; + /** + * Derivative seed of each component for AD + */ + Array _seed; /** * Flag whether DOF values are requested */ bool _dof = false; }; -class VariableGradient +template +VariableValueTempl::VariableValueTempl(const VariableValueTempl & object) + : _var(object._var), _seed(object._seed), _dof(object._dof) +{ + if constexpr (is_ad) + if (_var.coupled()) + { + if (!_seed.isAlloc()) + _seed.create(_var.components()); + + for (unsigned int comp = 0; comp < _var.components(); ++comp) + _seed[comp] = + _var.dot() ? _var.mooseVar(comp)->sys().duDotDu(_var.var(comp)) : (_var.old() ? 0 : 1); + + _seed.copyToDevice(); + } +} + +template +VariableValueTempl & +VariableValueTempl::operator=(const VariableValueTempl & object) +{ + _var = object._var; + _dof = object._dof; + + return *this; +} + +template +KOKKOS_FUNCTION auto +VariableValueTempl::operator()(AssemblyDatum & datum, + unsigned int idx, + unsigned int comp) const +{ + if constexpr (is_ad) + { + Real seed = + datum.do_derivatives() && _var.coupled() && _var.sys(comp) == datum.sys() ? _seed[comp] : 0; + + return get(datum, idx, comp, seed); + } + else + return get(datum, idx, comp); +} + +template +KOKKOS_FUNCTION auto +VariableValueTempl::get(Datum & datum, + unsigned int idx, + unsigned int comp, + [[maybe_unused]] Real seed) const { + KOKKOS_ASSERT(_var.initialized()); + + real_type value; + + if (_var.coupled()) + { + auto & sys = datum.system(_var.sys(comp)); + auto var = _var.var(comp); + auto tag = _var.tag(); + + if (_dof) + { + unsigned int dof; + + if (datum.isNodal()) + { + auto node = datum.node(); + dof = sys.getNodeLocalDofIndex(node, 0, var); + } + else + { + auto elem = datum.elem().id; + dof = sys.getElemLocalDofIndex(elem, idx, var); + } + + if constexpr (is_ad) + value = sys.getVectorDofADValue(dof, tag, seed); + else + value = sys.getVectorDofValue(dof, tag); + } + else + { + auto & elem = datum.elem(); + auto side = datum.side(); + + if constexpr (is_ad) + value = side == libMesh::invalid_uint + ? sys.getVectorQpADValue(elem, datum.qpOffset(), idx, var, tag, seed) + : sys.getVectorQpADValueFace(elem, side, idx, var, tag, seed); + else + value = side == libMesh::invalid_uint + ? sys.getVectorQpValue(elem, datum.qpOffset() + idx, var, tag) + : sys.getVectorQpValueFace(elem, side, idx, var, tag); + } + } + else + value = _var.value(comp); + + return value; +} + +template +class VariableGradientTempl +{ + using real3_type = std::conditional_t; + public: /** * Default constructor */ - VariableGradient() = default; + VariableGradientTempl() = default; /** * Constructor * @param var The Kokkos variable */ - VariableGradient(Variable var) : _var(var) {} + VariableGradientTempl(Variable var) : _var(var) {} /** * Constructor * @param var The MOOSE variable * @param tag The vector tag name */ - VariableGradient(const MooseVariableBase & var, const TagName & tag = Moose::SOLUTION_TAG) + VariableGradientTempl(const MooseVariableFieldBase & var, + const TagName & tag = Moose::SOLUTION_TAG) : _var(var, tag) { } + /** + * Constructor + * @param vars The MOOSE variables + * @param tag The vector tag name + */ + ///@{ + VariableGradientTempl(const std::vector & vars, + const TagName & tag = Moose::SOLUTION_TAG) + : _var(vars, tag) + { + } + VariableGradientTempl(const std::vector & vars, + const TagName & tag = Moose::SOLUTION_TAG) + : _var(vars, tag) + { + } + ///@} + + /** + * Copy constructor for parallel dispatch + */ + VariableGradientTempl(const VariableGradientTempl & object); + /** + * Copy assignment operator + */ + VariableGradientTempl & operator=(const VariableGradientTempl & object); /** * Get whether the variable was coupled @@ -204,11 +403,24 @@ class VariableGradient /** * Get the current variable gradient * @param datum The Datum object of the current thread - * @param idx The local quadrature point or DOF index + * @param qp The local quadrature point index * @param comp The variable component * @returns The variable gradient */ - KOKKOS_FUNCTION Real3 operator()(Datum & datum, unsigned int idx, unsigned int comp = 0) const; + KOKKOS_FUNCTION auto operator()(Datum & datum, unsigned int qp, unsigned int comp = 0) const + { + return get(datum, qp, comp); + } + + /** + * Get the current variable gradient + * @param datum The AssemblyDatum object of the current thread + * @param qp The local quadrature point index + * @param comp The variable component + * @returns The variable gradient + */ + KOKKOS_FUNCTION auto + operator()(AssemblyDatum & datum, unsigned int qp, unsigned int comp = 0) const; /** * Get the Kokkos variable @@ -217,76 +429,126 @@ class VariableGradient KOKKOS_FUNCTION const Variable & variable() const { return _var; } private: + /** + * Get the current variable gradient + * @param datum The Datum object of the current thread + * @param qp The local quadrature point index + * @param comp The variable component + * @param seed The derivative seed (only meaningful for AD) + * @returns The variable gradient + */ + KOKKOS_FUNCTION auto + get(Datum & datum, unsigned int qp, unsigned int comp = 0, Real seed = 0) const; + /** * Coupled Kokkos variable */ Variable _var; + /** + * Derivative seed of each component for AD + */ + Array _seed; }; -///@} -KOKKOS_FUNCTION inline Real -VariableValue::operator()(Datum & datum, unsigned int idx, unsigned int comp) const +template +VariableGradientTempl::VariableGradientTempl(const VariableGradientTempl & object) + : _var(object._var), _seed(object._seed) { - KOKKOS_ASSERT(_var.initialized()); - - if (_var.coupled()) - { - auto & sys = datum.system(_var.sys(comp)); - auto var = _var.var(comp); - auto tag = _var.tag(); - - if (_dof) + if constexpr (is_ad) + if (_var.coupled()) { - unsigned int dof; + if (!_seed.isAlloc()) + _seed.create(_var.components()); - if (datum.isNodal()) - { - auto node = datum.node(); - dof = sys.getNodeLocalDofIndex(node, 0, var); - } - else - { - auto elem = datum.elem().id; - dof = sys.getElemLocalDofIndex(elem, idx, var); - } + for (unsigned int comp = 0; comp < _var.components(); ++comp) + _seed[comp] = + _var.dot() ? _var.mooseVar(comp)->sys().duDotDu(_var.var(comp)) : (_var.old() ? 0 : 1); - return sys.getVectorDofValue(dof, tag); + _seed.copyToDevice(); } - else - { - auto & elem = datum.elem(); - auto side = datum.side(); - auto offset = datum.qpOffset(); +} - return side == libMesh::invalid_uint ? sys.getVectorQpValue(elem, offset + idx, var, tag) - : sys.getVectorQpValueFace(elem, side, idx, var, tag); - } +template +VariableGradientTempl & +VariableGradientTempl::operator=(const VariableGradientTempl & object) +{ + _var = object._var; + + return *this; +} + +template +KOKKOS_FUNCTION auto +VariableGradientTempl::operator()(AssemblyDatum & datum, + unsigned int qp, + unsigned int comp) const +{ + if constexpr (is_ad) + { + Real seed = + datum.do_derivatives() && _var.coupled() && _var.sys(comp) == datum.sys() ? _seed[comp] : 0; + + return get(datum, qp, comp, seed); } else - return _var.value(comp); + return get(datum, qp, comp); } -KOKKOS_FUNCTION inline Real3 -VariableGradient::operator()(Datum & datum, unsigned int qp, unsigned int comp) const +template +KOKKOS_FUNCTION auto +VariableGradientTempl::get(Datum & datum, + unsigned int qp, + unsigned int comp, + [[maybe_unused]] Real seed) const { KOKKOS_ASSERT(_var.initialized()); + real3_type grad; + if (_var.coupled()) { KOKKOS_ASSERT(!datum.isNodal()); auto & elem = datum.elem(); auto side = datum.side(); - auto offset = datum.qpOffset(); - return side == libMesh::invalid_uint - ? datum.system(_var.sys(comp)) - .getVectorQpGrad(elem, offset + qp, _var.var(comp), _var.tag()) - : datum.system(_var.sys(comp)) - .getVectorQpGradFace(elem, side, datum.J(qp), qp, _var.var(comp), _var.tag()); + if constexpr (is_ad) + grad = + side == libMesh::invalid_uint + ? datum.system(_var.sys(comp)) + .getVectorQpADGrad( + elem, datum.J(qp), datum.qpOffset(), qp, _var.var(comp), _var.tag(), seed) + : datum.system(_var.sys(comp)) + .getVectorQpADGradFace( + elem, side, datum.J(qp), qp, _var.var(comp), _var.tag(), seed); + else + grad = + side == libMesh::invalid_uint + ? datum.system(_var.sys(comp)) + .getVectorQpGrad(elem, datum.qpOffset() + qp, _var.var(comp), _var.tag()) + : datum.system(_var.sys(comp)) + .getVectorQpGradFace(elem, side, datum.J(qp), qp, _var.var(comp), _var.tag()); } - else - return Real3(0); + + return grad; } +using VariableValue = VariableValueTempl; +using ADVariableValue = VariableValueTempl; +using VariableGradient = VariableGradientTempl; +using ADVariableGradient = VariableGradientTempl; + +template <> +struct ArrayDeepCopy +{ + static constexpr bool value = true; +}; + +template <> +struct ArrayDeepCopy +{ + static constexpr bool value = true; +}; +///@} + } // namespace Moose::Kokkos diff --git a/framework/include/kokkos/bcs/KokkosADCoupledVarNeumannBC.h b/framework/include/kokkos/bcs/KokkosADCoupledVarNeumannBC.h new file mode 100644 index 000000000000..9e76cce88af1 --- /dev/null +++ b/framework/include/kokkos/bcs/KokkosADCoupledVarNeumannBC.h @@ -0,0 +1,47 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosADIntegratedBC.h" + +/** + * Implements a Neumann BC where grad(u)=_coupled_var on the boundary. + * Uses the term produced from integrating the diffusion operator by parts. + */ +class KokkosADCoupledVarNeumannBC : public Moose::Kokkos::ADIntegratedBC +{ +public: + static InputParameters validParams(); + + KokkosADCoupledVarNeumannBC(const InputParameters & parameters); + + template + KOKKOS_FUNCTION Moose::Kokkos::ADReal + computeQpResidual(const unsigned int i, const unsigned int qp, AssemblyDatum & datum) const; + +protected: + /// Variable providing the value of grad(u) on the boundary. + const Moose::Kokkos::ADVariableValue _coupled_var; + + /// A coefficient that is multiplied with the residual contribution + const Real _coef; + + /// Scale factor + const Moose::Kokkos::ADVariableValue _scale_factor; +}; + +template +KOKKOS_FUNCTION Moose::Kokkos::ADReal +KokkosADCoupledVarNeumannBC::computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const +{ + return -_test(datum, i, qp) * _scale_factor(datum, qp) * _coef * _coupled_var(datum, qp); +} diff --git a/framework/include/kokkos/bcs/KokkosADIntegratedBC.h b/framework/include/kokkos/bcs/KokkosADIntegratedBC.h new file mode 100644 index 000000000000..6d87a82fa89a --- /dev/null +++ b/framework/include/kokkos/bcs/KokkosADIntegratedBC.h @@ -0,0 +1,141 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosIntegratedBCBase.h" + +namespace Moose::Kokkos +{ + +/** + * The base class for a user to derive their own Kokkos integrated boundary conditions using + * automatic differentiation (AD). + * + * The user should define computeQpResidual() as inlined public methods in their derived class (not + * virtual override). The signature of computeQpResidual() expected to be defined in the derived + * class is as follows: + * + * @tparam Derived The object type + * @param i The element-local DOF index + * @param qp The local quadrature point index + * @param datum The AssemblyDatum object of the current thread + * @returns The residual contribution + * + * template + * KOKKOS_FUNCTION Moose::Kokkos::ADReal computeQpResidual(const unsigned int i, + * const unsigned int qp, + * AssemblyDatum & datum) const; + * + * Note that computeQpJacobian() and computeQpOffDiagJacobian() are unused for AD integrated + * boundary conditions. + */ +class ADIntegratedBC : public IntegratedBCBase +{ +public: + static InputParameters validParams(); + + /** + * Constructor + */ + ADIntegratedBC(const InputParameters & parameters); + + virtual void computeResidual() override; + virtual void computeJacobian() override; + virtual void computeResidualAndJacobian() override; + + /** + * The parallel computation entry function called by Kokkos + */ + template + KOKKOS_FUNCTION void operator()(ResidualLoop, const ThreadID tid, const Derived & bc) const; + + /** + * Compute residual + * @param bc The boundary condition object of the final derived type + * @param datum The AssemblyDatum object of the current thread + */ + template + KOKKOS_FUNCTION void computeResidualInternal(const Derived & bc, AssemblyDatum & datum) const; + +protected: + /** + * Dispatch parallel calculation + */ + virtual void dispatch(); + + /** + * Current test function + */ + const ADVariableTestValue _test; + /** + * Gradient of the current test function + */ + const ADVariableTestGradient _grad_test; + /** + * Current shape function + */ + const ADVariablePhiValue _phi; + /** + * Gradient of the current shape function + */ + const ADVariablePhiGradient _grad_phi; + /** + * Current solution at quadrature points + */ + const ADVariableValue _u; + /** + * Gradient of the current solution at quadrature points + */ + const ADVariableGradient _grad_u; + /** + * Whether computing residual + */ + bool _computing_residual = false; + /** + * Whether computing Jacobian + */ + bool _computing_jacobian = false; +}; + +template +KOKKOS_FUNCTION void +ADIntegratedBC::operator()(ResidualLoop, const ThreadID tid, const Derived & bc) const +{ + auto [elem, side] = kokkosBoundaryElementSideID(_thread(tid, 1)); + + AssemblyDatum datum( + elem, side, kokkosAssembly(), kokkosSystems(), _kokkos_var, _kokkos_var.var()); + + datum.set_local_parallel(_thread(tid, 0), _thread.size(0)); + + datum.do_derivatives(_computing_jacobian); + + bc.computeResidualInternal(bc, datum); +} + +template +KOKKOS_FUNCTION void +ADIntegratedBC::computeResidualInternal(const Derived & bc, AssemblyDatum & datum) const +{ + for (unsigned int i = datum.local_thread_id(); i < datum.n_dofs(); i += datum.num_local_threads()) + { + ADReal local_re = 0; + + for (unsigned int qp = 0; qp < datum.n_qps(); ++qp) + local_re += datum.JxW(qp) * bc.template computeQpResidual(i, qp, datum); + + if (_computing_residual) + accumulateTaggedElementalResidual(local_re.value(), datum.elem().id, i); + if (_computing_jacobian) + accumulateTaggedElementalMatrix(local_re.derivatives(), datum, i); + } +} + +} // namespace Moose::Kokkos diff --git a/framework/include/kokkos/bcs/KokkosADNeumannBC.h b/framework/include/kokkos/bcs/KokkosADNeumannBC.h new file mode 100644 index 000000000000..c5cd63ddba53 --- /dev/null +++ b/framework/include/kokkos/bcs/KokkosADNeumannBC.h @@ -0,0 +1,37 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosADIntegratedBC.h" + +class KokkosADNeumannBC : public Moose::Kokkos::ADIntegratedBC +{ +public: + static InputParameters validParams(); + + KokkosADNeumannBC(const InputParameters & parameters); + + template + KOKKOS_FUNCTION Moose::Kokkos::ADReal + computeQpResidual(const unsigned int i, const unsigned int qp, AssemblyDatum & datum) const; + +private: + /// Value of grad(u) on the boundary. + const Real _value; +}; + +template +KOKKOS_FUNCTION Moose::Kokkos::ADReal +KokkosADNeumannBC::computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const +{ + return -_test(datum, i, qp) * _value; +} diff --git a/framework/include/kokkos/bcs/KokkosADNodalBC.h b/framework/include/kokkos/bcs/KokkosADNodalBC.h new file mode 100644 index 000000000000..72763a1c1130 --- /dev/null +++ b/framework/include/kokkos/bcs/KokkosADNodalBC.h @@ -0,0 +1,99 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosNodalBCBase.h" + +namespace Moose::Kokkos +{ + +/** + * The base class for a user to derive their own Kokkos nodal boundary conditions using automatic + * differentiation (AD). + * + * The user should define computeQpResidual() as inlined public methods in their derived class (not + * virtual override). The signature of computeQpResidual() expected to be defined in the derived + * class is as follows: + * + * @tparam Derived The object type + * @param qp The dummy quadrature point index (= 0) + * @param datum The AssemblyDatum object of the current thread + * @returns The residual contribution + * + * template + * KOKKOS_FUNCTION Moose::Kokkos::ADReal computeQpResidual(const unsigned int qp, + * AssemblyDatum & datum) const; + * + * Note that computeQpJacobian() and computeQpOffDiagJacobian() are unused for AD nodal boundary + * conditions. + */ +class ADNodalBC : public NodalBCBase +{ +public: + static InputParameters validParams(); + + /** + * Constructor + */ + ADNodalBC(const InputParameters & parameters); + + virtual void computeResidual() override; + virtual void computeJacobian() override; + virtual void computeResidualAndJacobian() override; + + /** + * The parallel computation entry function called by Kokkos + */ + template + KOKKOS_FUNCTION void operator()(ResidualLoop, const ThreadID tid, const Derived & bc) const; + +protected: + /** + * Dispatch parallel calculation + */ + virtual void dispatch(); + + /** + * Current solution at nodes + */ + const ADVariableValue _u; + /** + * Whether computing residual + */ + bool _computing_residual = false; + /** + * Whether computing Jacobian + */ + bool _computing_jacobian = false; +}; + +template +KOKKOS_FUNCTION void +ADNodalBC::operator()(ResidualLoop, const ThreadID tid, const Derived & bc) const +{ + auto node = kokkosBoundaryNodeID(tid); + auto & sys = kokkosSystem(_kokkos_var.sys()); + + if (!sys.isNodalDefined(node, _kokkos_var.var())) + return; + + AssemblyDatum datum(node, kokkosAssembly(), kokkosSystems(), _kokkos_var, _kokkos_var.var()); + + datum.do_derivatives(_computing_jacobian); + + ADReal local_re = bc.template computeQpResidual(0, datum); + + if (_computing_residual) + accumulateTaggedNodalResidual(false, local_re.value(), node); + if (_computing_jacobian) + accumulateTaggedNodalMatrix(false, local_re.derivatives(), node); +} + +} // namespace Moose::Kokkos diff --git a/framework/include/kokkos/bcs/KokkosDirichletBC.h b/framework/include/kokkos/bcs/KokkosDirichletBC.h index 4d154fd4a6a9..dfae6b32270b 100644 --- a/framework/include/kokkos/bcs/KokkosDirichletBC.h +++ b/framework/include/kokkos/bcs/KokkosDirichletBC.h @@ -11,12 +11,13 @@ #include "KokkosDirichletBCBase.h" -class KokkosDirichletBC : public Moose::Kokkos::DirichletBCBase +template +class KokkosDirichletBCTempl : public Moose::Kokkos::DirichletBCBaseTempl { public: static InputParameters validParams(); - KokkosDirichletBC(const InputParameters & parameters); + KokkosDirichletBCTempl(const InputParameters & parameters); KOKKOS_FUNCTION Real computeValue(const unsigned int /* qp */, AssemblyDatum & /* datum */) const { @@ -26,3 +27,6 @@ class KokkosDirichletBC : public Moose::Kokkos::DirichletBCBase protected: const Moose::Kokkos::Scalar _value; }; + +typedef KokkosDirichletBCTempl KokkosDirichletBC; +typedef KokkosDirichletBCTempl KokkosADDirichletBC; diff --git a/framework/include/kokkos/bcs/KokkosDirichletBCBase.h b/framework/include/kokkos/bcs/KokkosDirichletBCBase.h index 850c41b33eb5..f9f0d7b7a5e1 100644 --- a/framework/include/kokkos/bcs/KokkosDirichletBCBase.h +++ b/framework/include/kokkos/bcs/KokkosDirichletBCBase.h @@ -10,6 +10,7 @@ #pragma once #include "KokkosNodalBC.h" +#include "KokkosADNodalBC.h" namespace Moose::Kokkos { @@ -17,15 +18,18 @@ namespace Moose::Kokkos /** * The base Kokkos boundary condition of a Dirichlet type */ -class DirichletBCBase : public NodalBC +template +class DirichletBCBaseTempl : public std::conditional_t { + using real_type = std::conditional_t; + public: static InputParameters validParams(); /** * Constructor */ - DirichletBCBase(const InputParameters & parameters); + DirichletBCBaseTempl(const InputParameters & parameters); /** * Get whether the value is to be preset @@ -53,9 +57,19 @@ class DirichletBCBase : public NodalBC KOKKOS_FUNCTION void operator()(PresetLoop, const ThreadID tid, const Derived & bc) const; template - KOKKOS_FUNCTION Real computeQpResidual(const unsigned int qp, AssemblyDatum & datum) const; + KOKKOS_FUNCTION auto computeQpResidual(const unsigned int qp, AssemblyDatum & datum) const; + + using Base = std::conditional_t; + using Base::operator(); - using NodalBC::operator(); +protected: + using Base::_kokkos_var; + using Base::_u; + using Base::kokkosAssembly; + using Base::kokkosBoundaryNodeID; + using Base::kokkosSystem; + using Base::kokkosSystems; + using Base::numKokkosBoundaryNodes; private: /** @@ -68,9 +82,10 @@ class DirichletBCBase : public NodalBC TagID _solution_tag; }; +template template KOKKOS_FUNCTION void -DirichletBCBase::operator()(PresetLoop, const ThreadID tid, const Derived & bc) const +DirichletBCBaseTempl::operator()(PresetLoop, const ThreadID tid, const Derived & bc) const { auto node = kokkosBoundaryNodeID(tid); auto & sys = kokkosSystem(_kokkos_var.sys()); @@ -84,17 +99,25 @@ DirichletBCBase::operator()(PresetLoop, const ThreadID tid, const Derived & bc) sys.getVectorDofValue(dof, _solution_tag) = bc.computeValue(0, datum); } +template template -KOKKOS_FUNCTION Real -DirichletBCBase::computeQpResidual(const unsigned int qp, AssemblyDatum & datum) const +KOKKOS_FUNCTION auto +DirichletBCBaseTempl::computeQpResidual(const unsigned int qp, AssemblyDatum & datum) const { auto bc = static_cast(this); - return _u(datum, qp) - bc->computeValue(qp, datum); + return _u(datum, qp) - real_type(bc->computeValue(qp, datum)); } +typedef DirichletBCBaseTempl DirichletBCBase; +typedef DirichletBCBaseTempl ADDirichletBCBase; + } // namespace Moose::Kokkos #define registerKokkosDirichletBC(app, classname) \ - registerKokkosBoundaryCondition(app, classname); \ + registerKokkosResidualObject(app, classname); \ + registerKokkosAdditionalOperation(classname, PresetLoop) + +#define registerKokkosADDirichletBC(app, classname) \ + registerKokkosADResidualObject(app, classname); \ registerKokkosAdditionalOperation(classname, PresetLoop) diff --git a/framework/include/kokkos/kernels/KokkosADCoupledTimeDerivative.h b/framework/include/kokkos/kernels/KokkosADCoupledTimeDerivative.h new file mode 100644 index 000000000000..00865a75a09d --- /dev/null +++ b/framework/include/kokkos/kernels/KokkosADCoupledTimeDerivative.h @@ -0,0 +1,39 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosADKernel.h" + +/** + * This calculates the time derivative for a coupled variable + **/ +class KokkosADCoupledTimeDerivative : public Moose::Kokkos::ADKernel +{ +public: + static InputParameters validParams(); + + KokkosADCoupledTimeDerivative(const InputParameters & parameters); + + template + KOKKOS_FUNCTION Moose::Kokkos::ADReal + computeQpResidual(const unsigned int i, const unsigned int qp, AssemblyDatum & datum) const; + +protected: + const Moose::Kokkos::ADVariableValue _v_dot; +}; + +template +KOKKOS_FUNCTION Moose::Kokkos::ADReal +KokkosADCoupledTimeDerivative::computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const +{ + return _test(datum, i, qp) * _v_dot(datum, qp); +} diff --git a/framework/include/kokkos/kernels/KokkosADDiffusion.h b/framework/include/kokkos/kernels/KokkosADDiffusion.h new file mode 100644 index 000000000000..39d1b1a27cd7 --- /dev/null +++ b/framework/include/kokkos/kernels/KokkosADDiffusion.h @@ -0,0 +1,33 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosADKernel.h" + +class KokkosADDiffusion : public Moose::Kokkos::ADKernel +{ +public: + static InputParameters validParams(); + + KokkosADDiffusion(const InputParameters & parameters); + + template + KOKKOS_FUNCTION Moose::Kokkos::ADReal + computeQpResidual(const unsigned int i, const unsigned int qp, AssemblyDatum & datum) const; +}; + +template +KOKKOS_FUNCTION Moose::Kokkos::ADReal +KokkosADDiffusion::computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const +{ + return _grad_u(datum, qp) * _grad_test(datum, i, qp); +} diff --git a/framework/include/kokkos/kernels/KokkosADKernel.h b/framework/include/kokkos/kernels/KokkosADKernel.h new file mode 100644 index 000000000000..90c2008b1fcb --- /dev/null +++ b/framework/include/kokkos/kernels/KokkosADKernel.h @@ -0,0 +1,144 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosKernelBase.h" + +namespace Moose::Kokkos +{ + +/** + * The base class for a user to derive their own Kokkos kernels using automatic differentiation + * (AD). + * + * The user should define computeQpResidual() as inlined public method in their derived class (not + * virtual override). The signature of computeQpResidual() expected to be defined in the derived + * class is as follows: + * + * @tparam Derived The object type + * @param i The element-local DOF index + * @param qp The local quadrature point index + * @param datum The AssemblyDatum object of the current thread + * @returns The residual contribution + * + * template + * KOKKOS_FUNCTION Moose::Kokkos::ADReal computeQpResidual(const unsigned int i, + * const unsigned int qp, + * AssemblyDatum & datum) const; + * + * Note that computeQpJacobian() and computeQpOffDiagJacobian() are unused for AD kernels. + */ +class ADKernel : public KernelBase +{ +public: + static InputParameters validParams(); + + /** + * Constructor + */ + ADKernel(const InputParameters & parameters); + + virtual void computeResidual() override; + virtual void computeJacobian() override; + virtual void computeResidualAndJacobian() override; + + /** + * The parallel computation entry function called by Kokkos + */ + template + KOKKOS_FUNCTION void operator()(ResidualLoop, const ThreadID tid, const Derived & kernel) const; + + /** + * Compute residual + * @param kernel The kernel object of the final derived type + * @param datum The AssemblyDatum object of the current thread + */ + template + KOKKOS_FUNCTION void computeResidualInternal(const Derived & kernel, AssemblyDatum & datum) const; + +protected: + /** + * Dispatch parallel calculation + */ + virtual void dispatch(); + + /** + * Current test function + */ + const ADVariableTestValue _test; + /** + * Gradient of the current test function + */ + const ADVariableTestGradient _grad_test; + /** + * Current shape function + */ + const ADVariablePhiValue _phi; + /** + * Gradient of the current shape function + */ + const ADVariablePhiGradient _grad_phi; + /** + * Current solution at quadrature points + */ + const ADVariableValue _u; + /** + * Gradient of the current solution at quadrature points + */ + const ADVariableGradient _grad_u; + /** + * Whether computing residual + */ + bool _computing_residual = false; + /** + * Whether computing Jacobian + */ + bool _computing_jacobian = false; +}; + +template +KOKKOS_FUNCTION void +ADKernel::operator()(ResidualLoop, const ThreadID tid, const Derived & kernel) const +{ + auto elem = kokkosBlockElementID(_thread(tid, 1)); + + AssemblyDatum datum(elem, + libMesh::invalid_uint, + kokkosAssembly(), + kokkosSystems(), + _kokkos_var, + _kokkos_var.var()); + + datum.set_local_parallel(_thread(tid, 0), _thread.size(0)); + + datum.do_derivatives(_computing_jacobian); + + kernel.computeResidualInternal(kernel, datum); +} + +template +KOKKOS_FUNCTION void +ADKernel::computeResidualInternal(const Derived & kernel, AssemblyDatum & datum) const +{ + for (unsigned int i = datum.local_thread_id(); i < datum.n_dofs(); i += datum.num_local_threads()) + { + ADReal local_re = 0; + + for (unsigned int qp = 0; qp < datum.n_qps(); ++qp) + local_re += datum.JxW(qp) * kernel.template computeQpResidual(i, qp, datum); + + if (_computing_residual) + accumulateTaggedElementalResidual(local_re.value(), datum.elem().id, i); + if (_computing_jacobian) + accumulateTaggedElementalMatrix(local_re.derivatives(), datum, i); + } +} + +} // namespace Moose::Kokkos diff --git a/framework/include/kokkos/kernels/KokkosADTimeDerivative.h b/framework/include/kokkos/kernels/KokkosADTimeDerivative.h new file mode 100644 index 000000000000..b9fee32c6a18 --- /dev/null +++ b/framework/include/kokkos/kernels/KokkosADTimeDerivative.h @@ -0,0 +1,33 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosADTimeKernel.h" + +class KokkosADTimeDerivative : public Moose::Kokkos::ADTimeKernel +{ +public: + static InputParameters validParams(); + + KokkosADTimeDerivative(const InputParameters & parameters); + + template + KOKKOS_FUNCTION Moose::Kokkos::ADReal + computeQpResidual(const unsigned int i, const unsigned int qp, AssemblyDatum & datum) const; +}; + +template +KOKKOS_FUNCTION Moose::Kokkos::ADReal +KokkosADTimeDerivative::computeQpResidual(const unsigned int i, + const unsigned int qp, + AssemblyDatum & datum) const +{ + return _test(datum, i, qp) * _u_dot(datum, qp); +} diff --git a/framework/include/kokkos/kernels/KokkosADTimeKernel.h b/framework/include/kokkos/kernels/KokkosADTimeKernel.h new file mode 100644 index 000000000000..30dda62b341d --- /dev/null +++ b/framework/include/kokkos/kernels/KokkosADTimeKernel.h @@ -0,0 +1,37 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "KokkosADKernel.h" + +namespace Moose::Kokkos +{ + +/** + * The base class for Kokkos time-derivative kernels + */ +class ADTimeKernel : public ADKernel +{ +public: + static InputParameters validParams(); + + /** + * Constructor + */ + ADTimeKernel(const InputParameters & parameters); + +protected: + /** + * Time derivative of the current solution at quadrature points + */ + const ADVariableValue _u_dot; +}; + +} // namespace Moose::Kokkos diff --git a/framework/include/kokkos/systems/KokkosMatrix.h b/framework/include/kokkos/systems/KokkosMatrix.h index bfc0e7904a56..acf724f204c6 100644 --- a/framework/include/kokkos/systems/KokkosMatrix.h +++ b/framework/include/kokkos/systems/KokkosMatrix.h @@ -36,6 +36,10 @@ class Matrix * Destructor */ ~Matrix() { destroy(); } + /** + * Get PETSc matrix handle + */ + Mat mat() { return _matrix; } /** * Free all data and reset */ diff --git a/framework/include/kokkos/systems/KokkosSystem.h b/framework/include/kokkos/systems/KokkosSystem.h index ee118a39824b..e0209174ca3b 100644 --- a/framework/include/kokkos/systems/KokkosSystem.h +++ b/framework/include/kokkos/systems/KokkosSystem.h @@ -287,10 +287,20 @@ class System : public MeshHolder, public AssemblyHolder * @param tag The vector tag * @returns The DOF value */ - KOKKOS_FUNCTION Real & getVectorDofValue(dof_id_type dof, TagID tag) const + KOKKOS_FUNCTION Real & getVectorDofValue(const dof_id_type dof, const TagID tag) const { return _vectors[tag][dof]; } + /** + * Get the DOF value of a tagged vector for automatic differentiation (AD) + * @param dof The local DOF index + * @param tag The vector tag + * @param seed The derivative seed + * @returns The DOF AD value with optional seed derivative + */ + KOKKOS_FUNCTION ADReal getVectorDofADValue(const dof_id_type dof, + const TagID tag, + const Real seed) const; /** * Get the quadrature value of a variable from a tagged vector * @param info The element information object @@ -299,11 +309,29 @@ class System : public MeshHolder, public AssemblyHolder * @param tag The vector tag * @returns The quadrature value */ - KOKKOS_FUNCTION Real & - getVectorQpValue(ElementInfo info, dof_id_type qp, unsigned int var, TagID tag) const + KOKKOS_FUNCTION Real & getVectorQpValue(const ElementInfo info, + const dof_id_type qp, + const unsigned int var, + const TagID tag) const { return _qp_solutions[tag](info.subdomain, var)[qp]; } + /** + * Get the quadrature value of a variable from a tagged vector for automatic differentiation (AD) + * @param info The element information object + * @param offset The offset into the global quadrature point index + * @param qp The local quadrature point index + * @param var The variable number + * @param tag The vector tag + * @param seed The derivative seed + * @returns The quadrature AD value + */ + KOKKOS_FUNCTION ADReal getVectorQpADValue(const ElementInfo info, + const dof_id_type offset, + const dof_id_type qp, + const unsigned int var, + const TagID tag, + const Real seed) const; /** * Get the quadrature gradient of a variable from a tagged vector * @param info The element information object @@ -312,11 +340,32 @@ class System : public MeshHolder, public AssemblyHolder * @param tag The vector tag * @returns The quadrature gradient */ - KOKKOS_FUNCTION Real3 & - getVectorQpGrad(ElementInfo info, dof_id_type qp, unsigned int var, TagID tag) const + KOKKOS_FUNCTION Real3 & getVectorQpGrad(const ElementInfo info, + const dof_id_type qp, + const unsigned int var, + const TagID tag) const { return _qp_solutions_grad[tag](info.subdomain, var)[qp]; } + /** + * Get the quadrature gradient of a variable from a tagged vector for automatic differentiation + * (AD) + * @param info The element information object + * @param jacobian The inverse Jacobian matrix + * @param offset The offset into the global quadrature point index + * @param qp The local quadrature point index + * @param var The variable number + * @param tag The vector tag + * @param seed The derivative seed + * @returns The quadrature AD gradient + */ + KOKKOS_FUNCTION ADReal3 getVectorQpADGrad(const ElementInfo info, + const Real33 jacobian, + const dof_id_type offset, + const dof_id_type qp, + const unsigned int var, + const TagID tag, + const Real seed) const; /** * Get the face quadrature value of a variable from a tagged vector * @param info The element information object @@ -331,6 +380,23 @@ class System : public MeshHolder, public AssemblyHolder const unsigned int qp, const unsigned int var, const TagID tag) const; + /** + * Get the face quadrature value of a variable from a tagged vector for automatic differentiation + * (AD) + * @param info The element information object + * @param side The side index + * @param qp The local quadrature point index + * @param var The vriable number + * @param tag The vector tag + * @param seed The derivative seed + * @returns The face quadrature AD value + */ + KOKKOS_FUNCTION ADReal getVectorQpADValueFace(const ElementInfo info, + const unsigned int side, + const unsigned int qp, + const unsigned int var, + const TagID tag, + const Real seed) const; /** * Get the face quadrature gradient of a variable from a tagged vector * @param info The element information object @@ -347,6 +413,25 @@ class System : public MeshHolder, public AssemblyHolder const unsigned int qp, const unsigned int var, const TagID tag) const; + /** + * Get the face quadrature gradient of a variable from a tagged vector for automatic + * differentiation (AD) + * @param info The element information object + * @param side The side index + * @param jacobian The inverse Jacobian matrix + * @param qp The local quadrature point index + * @param var The variable number + * @param tag The vector tag + * @param seed The derivative seed + * @returns The face quadrature AD gradient + */ + KOKKOS_FUNCTION ADReal3 getVectorQpADGradFace(const ElementInfo info, + const unsigned int side, + const Real33 jacobian, + const unsigned int qp, + const unsigned int var, + const TagID tag, + const Real seed) const; /** * Get an entry from a tagged matrix * @param row The local row index @@ -512,9 +597,75 @@ class System : public MeshHolder, public AssemblyHolder }; #ifdef MOOSE_KOKKOS_SCOPE +KOKKOS_FUNCTION inline ADReal +System::getVectorDofADValue(const dof_id_type dof, TagID tag, const Real seed) const +{ + ADReal value = _vectors[tag][dof]; + + if (seed != 0) + value.derivatives().insert(_local_to_global_dof_index[dof]) = seed; + + return value; +} + +KOKKOS_FUNCTION inline ADReal +System::getVectorQpADValue(const ElementInfo info, + const dof_id_type offset, + const dof_id_type qp, + const unsigned int var, + const TagID tag, + const Real seed) const +{ + ADReal value = 0; + + if (seed == 0) + value = getVectorQpValue(info, offset + qp, var, tag); + else + { + auto fe = _var_fe_types[var]; + auto n_dofs = kokkosAssembly().getNumDofs(info.type, fe); + auto & phi = kokkosAssembly().getPhi(info.subdomain, info.type, fe); + + for (unsigned int i = 0; i < n_dofs; ++i) + value += getVectorDofADValue(getElemLocalDofIndex(info.id, i, var), tag, seed) * phi(i, qp); + } + + return value; +} + +KOKKOS_FUNCTION inline ADReal3 +System::getVectorQpADGrad(const ElementInfo info, + const Real33 jacobian, + const dof_id_type offset, + const dof_id_type qp, + const unsigned int var, + const TagID tag, + const Real seed) const +{ + ADReal3 grad; + + if (seed == 0) + grad = getVectorQpGrad(info, offset + qp, var, tag); + else + { + auto fe = _var_fe_types[var]; + auto n_dofs = kokkosAssembly().getNumDofs(info.type, fe); + auto & grad_phi = kokkosAssembly().getGradPhi(info.subdomain, info.type, fe); + + for (unsigned int i = 0; i < n_dofs; ++i) + grad += getVectorDofADValue(getElemLocalDofIndex(info.id, i, var), tag, seed) * + (jacobian * grad_phi(i, qp)); + } + + return grad; +} + KOKKOS_FUNCTION inline Real -System::getVectorQpValueFace( - ElementInfo info, unsigned int side, unsigned int qp, unsigned int var, TagID tag) const +System::getVectorQpValueFace(const ElementInfo info, + const unsigned int side, + const unsigned int qp, + const unsigned int var, + const TagID tag) const { auto fe = _var_fe_types[var]; auto n_dofs = kokkosAssembly().getNumDofs(info.type, fe); @@ -527,13 +678,34 @@ System::getVectorQpValueFace( return value; } + +KOKKOS_FUNCTION inline ADReal +System::getVectorQpADValueFace(const ElementInfo info, + const unsigned int side, + const unsigned int qp, + const unsigned int var, + const TagID tag, + const Real seed) const +{ + auto fe = _var_fe_types[var]; + auto n_dofs = kokkosAssembly().getNumDofs(info.type, fe); + auto & phi = kokkosAssembly().getPhiFace(info.subdomain, info.type, fe)(side); + + ADReal value = 0; + + for (unsigned int i = 0; i < n_dofs; ++i) + value += getVectorDofADValue(getElemLocalDofIndex(info.id, i, var), tag, seed) * phi(i, qp); + + return value; +} + KOKKOS_FUNCTION inline Real3 -System::getVectorQpGradFace(ElementInfo info, - unsigned int side, - Real33 jacobian, - unsigned int qp, - unsigned int var, - TagID tag) const +System::getVectorQpGradFace(const ElementInfo info, + const unsigned int side, + const Real33 jacobian, + const unsigned int qp, + const unsigned int var, + const TagID tag) const { auto fe = _var_fe_types[var]; auto n_dofs = kokkosAssembly().getNumDofs(info.type, fe); @@ -547,6 +719,28 @@ System::getVectorQpGradFace(ElementInfo info, return grad; } + +KOKKOS_FUNCTION inline ADReal3 +System::getVectorQpADGradFace(const ElementInfo info, + const unsigned int side, + const Real33 jacobian, + const unsigned int qp, + const unsigned int var, + const TagID tag, + const Real seed) const +{ + auto fe = _var_fe_types[var]; + auto n_dofs = kokkosAssembly().getNumDofs(info.type, fe); + auto & grad_phi = kokkosAssembly().getGradPhiFace(info.subdomain, info.type, fe)(side); + + ADReal3 grad = ADReal(0); + + for (unsigned int i = 0; i < n_dofs; ++i) + grad += getVectorDofADValue(getElemLocalDofIndex(info.id, i, var), tag, seed) * + (jacobian * grad_phi(i, qp)); + + return grad; +} #endif /** diff --git a/framework/include/kokkos/utils/KokkosADReal.h b/framework/include/kokkos/utils/KokkosADReal.h new file mode 100644 index 000000000000..08336318e044 --- /dev/null +++ b/framework/include/kokkos/utils/KokkosADReal.h @@ -0,0 +1,30 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "ADRealForward.h" + +namespace Moose::Kokkos +{ + +using MetaPhysicL::KokkosSemiDynamicSparseNumberArray; + +typedef KokkosSemiDynamicSparseNumberArray> + DNDerivativeType; + +template +using DNDerivativeSize = + KokkosSemiDynamicSparseNumberArray>; + +typedef DualNumber ADReal; + +} // namespace Moose::Kokkos diff --git a/framework/include/systems/NonlinearSystemBase.h b/framework/include/systems/NonlinearSystemBase.h index 3b13512533dd..25c5eaa0ba2d 100644 --- a/framework/include/systems/NonlinearSystemBase.h +++ b/framework/include/systems/NonlinearSystemBase.h @@ -347,6 +347,11 @@ class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface void computeResidualAndJacobianInternal(const std::set & vector_tags, const std::set & matrix_tags); +#ifdef MOOSE_KOKKOS_ENABLED + void computeKokkosResidualAndJacobian(const std::set & vector_tags, + const std::set & matrix_tags); +#endif + /** * Form a residual vector for a given tag */ @@ -775,34 +780,43 @@ class NonlinearSystemBase : public SolverSystem, public PerfGraphInterface */ void computeResidualInternal(const std::set & tags); +#ifdef MOOSE_KOKKOS_ENABLED /** * Compute residual with Kokkos objects */ -#ifdef MOOSE_KOKKOS_ENABLED void computeKokkosResidual(const std::set & tags); - void computeKokkosNodalBCs(const std::set & tags); + /** + * Compute Kokkos nodal BCs + */ + void computeKokkosNodalBCsResidual(const std::set & tags); #endif /** * Enforces nodal boundary conditions. The boundary condition will be implemented * in the residual using all the tags in the system. */ - void computeNodalBCs(NumericVector & residual); + void computeNodalBCsResidual(NumericVector & residual); /** * Form a residual for BCs that at least has one of the given tags. */ - void computeNodalBCs(NumericVector & residual, const std::set & tags); + void computeNodalBCsResidual(NumericVector & residual, const std::set & tags); /** * Form multiple tag-associated residual vectors for the given tags. */ - void computeNodalBCs(const std::set & tags); + void computeNodalBCsResidual(const std::set & tags); /** - * compute the residual and Jacobian for nodal boundary conditions + * Compute the Jacobian for nodal boundary conditions */ - void computeNodalBCsResidualAndJacobian(); + void computeNodalBCsJacobian(const std::set & tags); + + /** + * Compute the residual and Jacobian together for nodal boundary conditions + */ + void computeNodalBCsResidualAndJacobian(const std::set & vector_tags, + const std::set & matrix_tags); /** * Form multiple matrices for all the tags. Users should not call this func directly. diff --git a/framework/kokkos.mk b/framework/kokkos.mk index 46c86b60c917..8c7c010235b5 100644 --- a/framework/kokkos.mk +++ b/framework/kokkos.mk @@ -110,7 +110,7 @@ else CXXFLAGS += -DMOOSE_ENABLE_KOKKOS_GPU=1 endif -KOKKOS_CXXFLAGS += -DMOOSE_KOKKOS_SCOPE=1 -fPIC +KOKKOS_CXXFLAGS += -DMOOSE_KOKKOS_SCOPE=1 -DMETAPHYSICL_KOKKOS_COMPILATION=1 -fPIC KOKKOS_LDFLAGS += $(libmesh_LDFLAGS) KOKKOS_INCLUDE = $(libmesh_INCLUDE) KOKKOS_LIBS = $(libmesh_LIBS) diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C index ae09fd1e342f..025294918932 100644 --- a/framework/src/base/Assembly.C +++ b/framework/src/base/Assembly.C @@ -3787,6 +3787,16 @@ Assembly::elementVolume(const Elem * elem) const return vol; } +void +Assembly::saveLocalADArray(std::vector & re, + unsigned int i, + unsigned int ntest, + const ADRealEigenVector & v) const +{ + for (unsigned int j = 0; j < v.size(); ++j, i += ntest) + re[i] += v(j); +} + void Assembly::addCachedJacobian(GlobalDataKey) { diff --git a/framework/src/interfaces/TaggingInterface.C b/framework/src/interfaces/TaggingInterface.C index 29b679be2c52..2bb4bb7518d2 100644 --- a/framework/src/interfaces/TaggingInterface.C +++ b/framework/src/interfaces/TaggingInterface.C @@ -460,6 +460,24 @@ TaggingInterface::assignTaggedLocalMatrix() *ke = _local_ke; } +void +TaggingInterface::addResiduals(Assembly & assembly, const ADResidualsPacket & packet) +{ + addResiduals(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor); +} + +void +TaggingInterface::addResidualsAndJacobian(Assembly & assembly, const ADResidualsPacket & packet) +{ + addResidualsAndJacobian(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor); +} + +void +TaggingInterface::addJacobian(Assembly & assembly, const ADResidualsPacket & packet) +{ + addJacobian(assembly, packet.residuals, packet.dof_indices, packet.scaling_factor); +} + #ifndef NDEBUG void TaggingInterface::checkForNans() const diff --git a/framework/src/kokkos/base/KokkosResidualObject.K b/framework/src/kokkos/base/KokkosResidualObject.K index 20fbc513a1f2..84cbb127ba99 100644 --- a/framework/src/kokkos/base/KokkosResidualObject.K +++ b/framework/src/kokkos/base/KokkosResidualObject.K @@ -52,6 +52,9 @@ ResidualObject::ResidualObject(const InputParameters & parameters, _dt(TransientInterface::_dt), _dt_old(TransientInterface::_dt_old) { + if (_var.isVector()) + paramError("variable", "Kokkos residual objects do not support vector variables yet."); + _kokkos_var.init(_var); } diff --git a/framework/src/kokkos/base/KokkosVariable.K b/framework/src/kokkos/base/KokkosVariable.K index b00b73837611..e7bbfca3c45a 100644 --- a/framework/src/kokkos/base/KokkosVariable.K +++ b/framework/src/kokkos/base/KokkosVariable.K @@ -16,13 +16,13 @@ namespace Moose::Kokkos { void -Variable::init(const MooseVariableBase & variable, const TagName & tag_name) +Variable::init(const MooseVariableFieldBase & variable, const TagName & tag_name) { init(variable, variable.sys().feProblem().getVectorTagID(tag_name)); } void -Variable::init(const MooseVariableBase & variable, const TagID tag) +Variable::init(const MooseVariableFieldBase & variable, const TagID tag) { _initialized = true; _coupled = true; @@ -30,23 +30,54 @@ Variable::init(const MooseVariableBase & variable, const TagID tag) _components = variable.count(); _tag = tag; + _moose_var.createHost(_components); _var.create(_components); _sys.create(_components); for (unsigned int comp = 0; comp < _components; ++comp) { + _moose_var[comp] = &variable; _var[comp] = variable.number() + comp; _sys[comp] = variable.sys().number(); } _var.copyToDevice(); _sys.copyToDevice(); + + const auto & problem = variable.sys().feProblem(); + + if (problem.vectorTagExists(Moose::SOLUTION_DOT_TAG)) + _dot = tag == problem.getVectorTagID(Moose::SOLUTION_DOT_TAG); + + if (problem.vectorTagExists(Moose::OLD_SOLUTION_TAG)) + _old = tag == problem.getVectorTagID(Moose::OLD_SOLUTION_TAG); + + if (problem.vectorTagExists(Moose::OLDER_SOLUTION_TAG)) + _old = _old || tag == problem.getVectorTagID(Moose::OLDER_SOLUTION_TAG); +} + +void +Variable::init(const std::vector & variables, const TagName & tag_name) +{ + init(std::vector(variables.begin(), variables.end()), + variables[0]->sys().feProblem().getVectorTagID(tag_name)); +} + +void +Variable::init(const std::vector & variables, + const TagName & tag_name) +{ + init(variables, variables[0]->sys().feProblem().getVectorTagID(tag_name)); +} + +void +Variable::init(const std::vector & variables, const TagID tag) +{ + init(std::vector(variables.begin(), variables.end()), tag); } void -Variable::init(const std::vector & variables, - const TagID tag, - CoupleableKey) +Variable::init(const std::vector & variables, const TagID tag) { _initialized = true; _coupled = true; @@ -54,6 +85,7 @@ Variable::init(const std::vector & variables, _components = variables.size(); _tag = tag; + _moose_var.createHost(_components); _var.create(_components); _sys.create(_components); @@ -61,12 +93,24 @@ Variable::init(const std::vector & variables, { _nodal = _nodal && variables[comp]->isNodal(); + _moose_var[comp] = variables[comp]; _var[comp] = variables[comp]->number(); _sys[comp] = variables[comp]->sys().number(); } _var.copyToDevice(); _sys.copyToDevice(); + + const auto & problem = variables[0]->sys().feProblem(); + + if (problem.vectorTagExists(Moose::SOLUTION_DOT_TAG)) + _dot = tag == problem.getVectorTagID(Moose::SOLUTION_DOT_TAG); + + if (problem.vectorTagExists(Moose::OLD_SOLUTION_TAG)) + _old = tag == problem.getVectorTagID(Moose::OLD_SOLUTION_TAG); + + if (problem.vectorTagExists(Moose::OLDER_SOLUTION_TAG)) + _old = _old || tag == problem.getVectorTagID(Moose::OLDER_SOLUTION_TAG); } void diff --git a/framework/src/kokkos/bcs/KokkosADCoupledVarNeumannBC.K b/framework/src/kokkos/bcs/KokkosADCoupledVarNeumannBC.K new file mode 100644 index 000000000000..0fd9fa4e31c7 --- /dev/null +++ b/framework/src/kokkos/bcs/KokkosADCoupledVarNeumannBC.K @@ -0,0 +1,34 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADCoupledVarNeumannBC.h" + +registerKokkosADResidualObject("MooseApp", KokkosADCoupledVarNeumannBC); + +InputParameters +KokkosADCoupledVarNeumannBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addRequiredCoupledVar("v", "Coupled variable setting the gradient on the boundary."); + params.addCoupledVar("scale_factor", 1., "Scale factor to multiply the heat flux with"); + params.addParam( + "coef", 1.0, "Coefficent ($\\sigma$) multiplier for the coupled force term."); + params.addClassDescription("Imposes the integrated boundary condition " + "$\\frac{\\partial u}{\\partial n}=v$, " + "where $v$ is a variable."); + return params; +} + +KokkosADCoupledVarNeumannBC::KokkosADCoupledVarNeumannBC(const InputParameters & parameters) + : ADIntegratedBC(parameters), + _coupled_var(kokkosADCoupledValue("v")), + _coef(getParam("coef")), + _scale_factor(kokkosADCoupledValue("scale_factor")) +{ +} diff --git a/framework/src/kokkos/bcs/KokkosADIntegratedBC.K b/framework/src/kokkos/bcs/KokkosADIntegratedBC.K new file mode 100644 index 000000000000..89a6a1f31a3d --- /dev/null +++ b/framework/src/kokkos/bcs/KokkosADIntegratedBC.K @@ -0,0 +1,74 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADIntegratedBC.h" + +namespace Moose::Kokkos +{ + +InputParameters +ADIntegratedBC::validParams() +{ + InputParameters params = IntegratedBCBase::validParams(); + return params; +} + +ADIntegratedBC::ADIntegratedBC(const InputParameters & parameters) + : IntegratedBCBase(parameters, Moose::VarFieldType::VAR_FIELD_STANDARD), + _test(), + _grad_test(), + _phi(), + _grad_phi(), + _u(_var), + _grad_u(_var) +{ + addMooseVariableDependency(&_var); +} + +void +ADIntegratedBC::computeResidual() +{ + _computing_residual = true; + _computing_jacobian = false; + + dispatch(); +} + +void +ADIntegratedBC::computeJacobian() +{ + _computing_residual = false; + _computing_jacobian = true; + + dispatch(); +} + +void +ADIntegratedBC::computeResidualAndJacobian() +{ + _computing_residual = true; + _computing_jacobian = true; + + dispatch(); +} + +void +ADIntegratedBC::dispatch() +{ + _thread.resize(_num_local_threads, numKokkosBoundarySides()); + + Policy policy(0, _thread.size()); + + if (!_residual_dispatcher) + _residual_dispatcher = DispatcherRegistry::build(this, type()); + + _residual_dispatcher->parallelFor(policy); +} + +} // namespace Moose::Kokkos diff --git a/framework/src/kokkos/bcs/KokkosADNeumannBC.K b/framework/src/kokkos/bcs/KokkosADNeumannBC.K new file mode 100644 index 000000000000..2b0f10226d79 --- /dev/null +++ b/framework/src/kokkos/bcs/KokkosADNeumannBC.K @@ -0,0 +1,32 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADNeumannBC.h" + +registerKokkosADResidualObject("MooseApp", KokkosADNeumannBC); + +InputParameters +KokkosADNeumannBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addParam("value", + 0.0, + "For a Laplacian problem, the value of the gradient dotted with the " + "normals on the boundary."); + params.declareControllable("value"); + params.addClassDescription("Imposes the integrated boundary condition " + "$\\frac{\\partial u}{\\partial n}=h$, " + "where $h$ is a constant, controllable value."); + return params; +} + +KokkosADNeumannBC::KokkosADNeumannBC(const InputParameters & parameters) + : ADIntegratedBC(parameters), _value(getParam("value")) +{ +} diff --git a/framework/src/kokkos/bcs/KokkosADNodalBC.K b/framework/src/kokkos/bcs/KokkosADNodalBC.K new file mode 100644 index 000000000000..e1b74090a411 --- /dev/null +++ b/framework/src/kokkos/bcs/KokkosADNodalBC.K @@ -0,0 +1,67 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADNodalBC.h" + +namespace Moose::Kokkos +{ + +InputParameters +ADNodalBC::validParams() +{ + InputParameters params = NodalBCBase::validParams(); + return params; +} + +ADNodalBC::ADNodalBC(const InputParameters & parameters) + : NodalBCBase(parameters, Moose::VarFieldType::VAR_FIELD_STANDARD), + _u(_var, Moose::SOLUTION_TAG, true) +{ + addMooseVariableDependency(&_var); +} + +void +ADNodalBC::computeResidual() +{ + _computing_residual = true; + _computing_jacobian = false; + + dispatch(); +} + +void +ADNodalBC::computeJacobian() +{ + _computing_residual = false; + _computing_jacobian = true; + + dispatch(); +} + +void +ADNodalBC::computeResidualAndJacobian() +{ + _computing_residual = true; + _computing_jacobian = true; + + dispatch(); +} + +void +ADNodalBC::dispatch() +{ + Policy policy(0, numKokkosBoundaryNodes()); + + if (!_residual_dispatcher) + _residual_dispatcher = DispatcherRegistry::build(this, type()); + + _residual_dispatcher->parallelFor(policy); +} + +} // namespace Moose::Kokkos diff --git a/framework/src/kokkos/bcs/KokkosDirectionalNeumannBC.K b/framework/src/kokkos/bcs/KokkosDirectionalNeumannBC.K index 39bf9413031d..15de677c2bd6 100644 --- a/framework/src/kokkos/bcs/KokkosDirectionalNeumannBC.K +++ b/framework/src/kokkos/bcs/KokkosDirectionalNeumannBC.K @@ -9,7 +9,7 @@ #include "KokkosDirectionalNeumannBC.h" -registerKokkosBoundaryCondition("MooseApp", KokkosDirectionalNeumannBC); +registerKokkosResidualObject("MooseApp", KokkosDirectionalNeumannBC); InputParameters KokkosDirectionalNeumannBC::validParams() diff --git a/framework/src/kokkos/bcs/KokkosDirichletBC.K b/framework/src/kokkos/bcs/KokkosDirichletBC.K index 01da93bbffef..c8e584adc02d 100644 --- a/framework/src/kokkos/bcs/KokkosDirichletBC.K +++ b/framework/src/kokkos/bcs/KokkosDirichletBC.K @@ -10,11 +10,13 @@ #include "KokkosDirichletBC.h" registerKokkosDirichletBC("MooseApp", KokkosDirichletBC); +registerKokkosADDirichletBC("MooseApp", KokkosADDirichletBC); +template InputParameters -KokkosDirichletBC::validParams() +KokkosDirichletBCTempl::validParams() { - InputParameters params = DirichletBCBase::validParams(); + InputParameters params = Moose::Kokkos::DirichletBCBaseTempl::validParams(); params.addRequiredParam("value", "Value of the BC"); params.declareControllable("value"); params.addClassDescription("Imposes the essential boundary condition $u=g$, where $g$ " @@ -22,7 +24,12 @@ KokkosDirichletBC::validParams() return params; } -KokkosDirichletBC::KokkosDirichletBC(const InputParameters & parameters) - : DirichletBCBase(parameters), _value(getParam("value")) +template +KokkosDirichletBCTempl::KokkosDirichletBCTempl(const InputParameters & parameters) + : Moose::Kokkos::DirichletBCBaseTempl(parameters), + _value(this->template getParam("value")) { } + +template class KokkosDirichletBCTempl; +template class KokkosDirichletBCTempl; diff --git a/framework/src/kokkos/bcs/KokkosDirichletBCBase.K b/framework/src/kokkos/bcs/KokkosDirichletBCBase.K index e0b11667d2be..121b3ffce13c 100644 --- a/framework/src/kokkos/bcs/KokkosDirichletBCBase.K +++ b/framework/src/kokkos/bcs/KokkosDirichletBCBase.K @@ -12,30 +12,36 @@ namespace Moose::Kokkos { +template InputParameters -DirichletBCBase::validParams() +DirichletBCBaseTempl::validParams() { - InputParameters params = NodalBC::validParams(); + InputParameters params = Base::validParams(); params.addParam( "preset", true, "Whether or not to preset the BC (apply the value before the solve begins)."); return params; } -DirichletBCBase::DirichletBCBase(const InputParameters & parameters) - : NodalBC(parameters), _preset(getParam("preset")) +template +DirichletBCBaseTempl::DirichletBCBaseTempl(const InputParameters & parameters) + : Base(parameters), _preset(this->template getParam("preset")) { } +template void -DirichletBCBase::presetSolution(TagID tag) +DirichletBCBaseTempl::presetSolution(TagID tag) { _solution_tag = tag; Policy policy(0, numKokkosBoundaryNodes()); - auto dispatcher = DispatcherRegistry::build(this, type()); + auto dispatcher = DispatcherRegistry::build(this, this->type()); dispatcher->parallelFor(policy); } +template class DirichletBCBaseTempl; +template class DirichletBCBaseTempl; + } // namespace Moose::Kokkos diff --git a/framework/src/kokkos/bcs/KokkosEigenDirichletBC.K b/framework/src/kokkos/bcs/KokkosEigenDirichletBC.K index 9244c7f72542..34e42924c7ff 100644 --- a/framework/src/kokkos/bcs/KokkosEigenDirichletBC.K +++ b/framework/src/kokkos/bcs/KokkosEigenDirichletBC.K @@ -9,7 +9,7 @@ #include "KokkosEigenDirichletBC.h" -registerKokkosBoundaryCondition("MooseApp", KokkosEigenDirichletBC); +registerKokkosResidualObject("MooseApp", KokkosEigenDirichletBC); InputParameters KokkosEigenDirichletBC::validParams() diff --git a/framework/src/kokkos/interfaces/KokkosCoupleable.K b/framework/src/kokkos/interfaces/KokkosCoupleable.K index 737b2d185740..040c75f8da1c 100644 --- a/framework/src/kokkos/interfaces/KokkosCoupleable.K +++ b/framework/src/kokkos/interfaces/KokkosCoupleable.K @@ -13,6 +13,8 @@ #include "SystemBase.h" #include "FEProblemBase.h" +using CoupleableKey = Moose::Kokkos::Variable::CoupleableKey; + Moose::Kokkos::Variable Coupleable::kokkosCoupledVectorTagVariable(const std::string & var_name, const std::string & tag_name, @@ -34,10 +36,10 @@ Coupleable::kokkosCoupledVectorTagVariable(const std::string & var_name, const_cast(this)->addFEVariableCoupleableVectorTag(tag); - variable.init({var}, tag, {}); + variable.init(*var, tag); } else - variable.init({_c_parameters.defaultCoupledValue(var_name, comp)}, {}); + variable.init({_c_parameters.defaultCoupledValue(var_name, comp)}, CoupleableKey{}); return variable; } @@ -52,7 +54,7 @@ Coupleable::kokkosCoupledVectorTagVariables(const std::string & var_name, if (isCoupled(var_name)) { - std::vector vars; + std::vector vars; for (unsigned int comp = 0; comp < components; ++comp) { @@ -71,7 +73,7 @@ Coupleable::kokkosCoupledVectorTagVariables(const std::string & var_name, const_cast(this)->addFEVariableCoupleableVectorTag(tag); - variable.init({vars}, tag, {}); + variable.init(vars, tag); } else { @@ -80,7 +82,7 @@ Coupleable::kokkosCoupledVectorTagVariables(const std::string & var_name, for (unsigned int comp = 0; comp < components; ++comp) default_values[comp] = _c_parameters.defaultCoupledValue(var_name, comp); - variable.init(default_values, {}); + variable.init(default_values, CoupleableKey{}); } return variable; @@ -91,7 +93,7 @@ Coupleable::kokkosZeroVariable() const { Moose::Kokkos::Variable variable; - variable.init({0}, {}); + variable.init({0}, CoupleableKey{}); return variable; } @@ -459,6 +461,371 @@ Coupleable::kokkosCoupledNodalDots(const std::string & var_name) const return kokkosCoupledVectorTagNodalValuesByName(var_name, Moose::SOLUTION_DOT_TAG); } +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagValueByName(const std::string & var_name, + const std::string & tag_name, + unsigned int comp) const +{ + auto variable = kokkosCoupledVectorTagVariable(var_name, tag_name, comp); + + return Moose::Kokkos::ADVariableValue(variable); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagValuesByName(const std::string & var_name, + const std::string & tag_name) const +{ + auto variable = kokkosCoupledVectorTagVariables(var_name, tag_name); + + return Moose::Kokkos::ADVariableValue(variable); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledVectorTagGradientByName(const std::string & var_name, + const std::string & tag_name, + unsigned int comp) const +{ + auto variable = kokkosCoupledVectorTagVariable(var_name, tag_name, comp); + + return Moose::Kokkos::ADVariableGradient(variable); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledVectorTagGradientsByName(const std::string & var_name, + const std::string & tag_name) const +{ + auto variable = kokkosCoupledVectorTagVariables(var_name, tag_name); + + return Moose::Kokkos::ADVariableGradient(variable); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagNodalValueByName(const std::string & var_name, + const std::string & tag_name, + unsigned int comp) const +{ + auto variable = kokkosCoupledVectorTagVariable(var_name, tag_name, comp); + + if (!variable.nodal()) + mooseError("Cannot get nodal values from the coupled variable '", + var_name, + "', because the associated variable is not nodal."); + + return Moose::Kokkos::ADVariableValue(variable, true); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagNodalValuesByName(const std::string & var_name, + const std::string & tag_name) const +{ + auto variable = kokkosCoupledVectorTagVariables(var_name, tag_name); + + if (!variable.nodal()) + mooseError("Cannot get nodal values from the coupled variable '", + var_name, + "', because the associated variable is not nodal."); + + return Moose::Kokkos::ADVariableValue(variable, true); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagDofValueByName(const std::string & var_name, + const std::string & tag_name, + unsigned int comp) const +{ + auto variable = kokkosCoupledVectorTagVariable(var_name, tag_name, comp); + + return Moose::Kokkos::ADVariableValue(variable, true); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagDofValuesByName(const std::string & var_name, + const std::string & tag_name) const +{ + auto variable = kokkosCoupledVectorTagVariables(var_name, tag_name); + + return Moose::Kokkos::ADVariableValue(variable, true); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagValue(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagValueByName(var_name, tag_name, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagValues(const std::string & var_name, + const std::string & tag_param_name) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagValuesByName(var_name, tag_name); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledVectorTagGradient(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagGradientByName(var_name, tag_name, comp); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledVectorTagGradients(const std::string & var_name, + const std::string & tag_param_name) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagGradientsByName(var_name, tag_name); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagNodalValue(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagNodalValueByName(var_name, tag_name, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagNodalValues(const std::string & var_name, + const std::string & tag_param_name) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagNodalValuesByName(var_name, tag_name); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagDofValue(const std::string & var_name, + const std::string & tag_param_name, + unsigned int comp) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagDofValueByName(var_name, tag_name, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledVectorTagDofValues(const std::string & var_name, + const std::string & tag_param_name) const +{ + if (!_c_parameters.isParamValid(tag_param_name)) + mooseError("Tag name parameter '", tag_param_name, "' is invalid"); + + TagName tag_name = _c_parameters.get(tag_param_name); + + return kokkosADCoupledVectorTagDofValuesByName(var_name, tag_name); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledValue(const std::string & var_name, unsigned int comp) const +{ + return _c_nodal ? kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::SOLUTION_TAG, comp) + : kokkosADCoupledVectorTagValueByName(var_name, Moose::SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledValues(const std::string & var_name) const +{ + return _c_nodal ? kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::SOLUTION_TAG) + : kokkosADCoupledVectorTagValuesByName(var_name, Moose::SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledGradient(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagGradientByName(var_name, Moose::SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledGradients(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagGradientsByName(var_name, Moose::SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalValue(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalValues(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDofValue(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagDofValueByName(var_name, Moose::SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDofValues(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagDofValuesByName(var_name, Moose::SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledValueOld(const std::string & var_name, unsigned int comp) const +{ + return _c_nodal + ? kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::OLD_SOLUTION_TAG, comp) + : kokkosADCoupledVectorTagValueByName(var_name, Moose::OLD_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledValuesOld(const std::string & var_name) const +{ + return _c_nodal ? kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::OLD_SOLUTION_TAG) + : kokkosADCoupledVectorTagValuesByName(var_name, Moose::OLD_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledGradientOld(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagGradientByName(var_name, Moose::OLD_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledGradientsOld(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagGradientsByName(var_name, Moose::OLD_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalValueOld(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::OLD_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalValuesOld(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::OLD_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDofValueOld(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagDofValueByName(var_name, Moose::OLD_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDofValuesOld(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagDofValuesByName(var_name, Moose::OLD_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledValueOlder(const std::string & var_name, unsigned int comp) const +{ + return _c_nodal + ? kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::OLDER_SOLUTION_TAG, comp) + : kokkosADCoupledVectorTagValueByName(var_name, Moose::OLDER_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledValuesOlder(const std::string & var_name) const +{ + return _c_nodal ? kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::OLDER_SOLUTION_TAG) + : kokkosADCoupledVectorTagValuesByName(var_name, Moose::OLDER_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledGradientOlder(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagGradientByName(var_name, Moose::OLDER_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableGradient +Coupleable::kokkosADCoupledGradientsOlder(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagGradientsByName(var_name, Moose::OLDER_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalValueOlder(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::OLDER_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalValuesOlder(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::OLDER_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDofValueOlder(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagDofValueByName(var_name, Moose::OLDER_SOLUTION_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDofValuesOlder(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagDofValuesByName(var_name, Moose::OLDER_SOLUTION_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDot(const std::string & var_name, unsigned int comp) const +{ + return _c_nodal + ? kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::SOLUTION_DOT_TAG, comp) + : kokkosADCoupledVectorTagValueByName(var_name, Moose::SOLUTION_DOT_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledDots(const std::string & var_name) const +{ + return _c_nodal ? kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::SOLUTION_DOT_TAG) + : kokkosADCoupledVectorTagValuesByName(var_name, Moose::SOLUTION_DOT_TAG); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalDot(const std::string & var_name, unsigned int comp) const +{ + return kokkosADCoupledVectorTagNodalValueByName(var_name, Moose::SOLUTION_DOT_TAG, comp); +} + +Moose::Kokkos::ADVariableValue +Coupleable::kokkosADCoupledNodalDots(const std::string & var_name) const +{ + return kokkosADCoupledVectorTagNodalValuesByName(var_name, Moose::SOLUTION_DOT_TAG); +} + Moose::Kokkos::Scalar Coupleable::kokkosCoupledDotDu(const std::string & var_name, unsigned int comp) const { diff --git a/framework/src/kokkos/kernels/KokkosADCoupledTimeDerivative.K b/framework/src/kokkos/kernels/KokkosADCoupledTimeDerivative.K new file mode 100644 index 000000000000..f9c864b71abe --- /dev/null +++ b/framework/src/kokkos/kernels/KokkosADCoupledTimeDerivative.K @@ -0,0 +1,27 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADCoupledTimeDerivative.h" + +registerKokkosADResidualObject("MooseApp", KokkosADCoupledTimeDerivative); + +InputParameters +KokkosADCoupledTimeDerivative::validParams() +{ + InputParameters params = ADKernel::validParams(); + params.addClassDescription("Time derivative Kernel that acts on a coupled variable. Weak form: " + "$(\\psi_i, \\frac{\\partial v_h}{\\partial t})$."); + params.addRequiredCoupledVar("v", "Coupled variable"); + return params; +} + +KokkosADCoupledTimeDerivative::KokkosADCoupledTimeDerivative(const InputParameters & parameters) + : ADKernel(parameters), _v_dot(kokkosADCoupledDot("v")) +{ +} diff --git a/framework/src/kokkos/kernels/KokkosADDiffusion.K b/framework/src/kokkos/kernels/KokkosADDiffusion.K new file mode 100644 index 000000000000..3036525def84 --- /dev/null +++ b/framework/src/kokkos/kernels/KokkosADDiffusion.K @@ -0,0 +1,23 @@ +//* This file is part of the MOOSE framework +//* https://mooseframework.inl.gov +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADDiffusion.h" + +registerKokkosADResidualObject("MooseApp", KokkosADDiffusion); + +InputParameters +KokkosADDiffusion::validParams() +{ + auto params = ADKernel::validParams(); + params.addClassDescription("Same as `KokkosDiffusion` in terms of physics/residual, but the " + "Jacobian is computed using forward automatic differentiation"); + return params; +} + +KokkosADDiffusion::KokkosADDiffusion(const InputParameters & parameters) : ADKernel(parameters) {} diff --git a/framework/src/kokkos/kernels/KokkosADKernel.K b/framework/src/kokkos/kernels/KokkosADKernel.K new file mode 100644 index 000000000000..0bb46458610d --- /dev/null +++ b/framework/src/kokkos/kernels/KokkosADKernel.K @@ -0,0 +1,74 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADKernel.h" + +namespace Moose::Kokkos +{ + +InputParameters +ADKernel::validParams() +{ + InputParameters params = KernelBase::validParams(); + return params; +} + +ADKernel::ADKernel(const InputParameters & parameters) + : KernelBase(parameters, Moose::VarFieldType::VAR_FIELD_STANDARD), + _test(), + _grad_test(), + _phi(), + _grad_phi(), + _u(_var), + _grad_u(_var) +{ + addMooseVariableDependency(&_var); +} + +void +ADKernel::computeResidual() +{ + _computing_residual = true; + _computing_jacobian = false; + + dispatch(); +} + +void +ADKernel::computeJacobian() +{ + _computing_residual = false; + _computing_jacobian = true; + + dispatch(); +} + +void +ADKernel::computeResidualAndJacobian() +{ + _computing_residual = true; + _computing_jacobian = true; + + dispatch(); +} + +void +ADKernel::dispatch() +{ + _thread.resize(_num_local_threads, numKokkosBlockElements()); + + Policy policy(0, _thread.size()); + + if (!_residual_dispatcher) + _residual_dispatcher = DispatcherRegistry::build(this, type()); + + _residual_dispatcher->parallelFor(policy); +} + +} // namespace Moose::Kokkos diff --git a/framework/src/kokkos/kernels/KokkosADTimeDerivative.K b/framework/src/kokkos/kernels/KokkosADTimeDerivative.K new file mode 100644 index 000000000000..d9a001756ae0 --- /dev/null +++ b/framework/src/kokkos/kernels/KokkosADTimeDerivative.K @@ -0,0 +1,26 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADTimeDerivative.h" + +registerKokkosADResidualObject("MooseApp", KokkosADTimeDerivative); + +InputParameters +KokkosADTimeDerivative::validParams() +{ + InputParameters params = ADTimeKernel::validParams(); + params.addClassDescription("The time derivative operator with the weak form of $(\\psi_i, " + "\\frac{\\partial u_h}{\\partial t})$."); + return params; +} + +KokkosADTimeDerivative::KokkosADTimeDerivative(const InputParameters & parameters) + : ADTimeKernel(parameters) +{ +} diff --git a/framework/src/kokkos/kernels/KokkosADTimeKernel.K b/framework/src/kokkos/kernels/KokkosADTimeKernel.K new file mode 100644 index 000000000000..5ab594801f5a --- /dev/null +++ b/framework/src/kokkos/kernels/KokkosADTimeKernel.K @@ -0,0 +1,31 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "KokkosADTimeKernel.h" + +namespace Moose::Kokkos +{ + +InputParameters +ADTimeKernel::validParams() +{ + InputParameters params = ADKernel::validParams(); + + params.set("vector_tags") = "time"; + params.set("matrix_tags") = "system time"; + + return params; +} + +ADTimeKernel::ADTimeKernel(const InputParameters & parameters) + : ADKernel(parameters), _u_dot(_var, Moose::SOLUTION_DOT_TAG) +{ +} + +} // namespace Moose::Kokkos diff --git a/framework/src/kokkos/systems/KokkosAuxiliarySystem.K b/framework/src/kokkos/systems/KokkosAuxiliarySystem.K index 9326c5037eb8..91350a0774e9 100644 --- a/framework/src/kokkos/systems/KokkosAuxiliarySystem.K +++ b/framework/src/kokkos/systems/KokkosAuxiliarySystem.K @@ -31,7 +31,8 @@ AuxiliarySystem::kokkosCompute(ExecFlagType type) { TIME_SECTION("computeKokkosAuxKernel", 3); - if (!_kokkos_elemental_aux_storage[type].size() && !_kokkos_nodal_aux_storage[type].size()) + if (!_kokkos_elemental_aux_storage[type].hasActiveObjects() && + !_kokkos_nodal_aux_storage[type].hasActiveObjects()) return; auto & systems = _fe_problem.getKokkosSystems(); @@ -104,10 +105,10 @@ AuxiliarySystem::kokkosCompute(ExecFlagType type) // Compute auxiliary kernels - for (auto nodal : _kokkos_nodal_aux_storage[type].getActiveObjects()) + for (auto & nodal : _kokkos_nodal_aux_storage[type].getActiveObjects()) nodal->compute(); - for (auto elemental : _kokkos_elemental_aux_storage[type].getActiveObjects()) + for (auto & elemental : _kokkos_elemental_aux_storage[type].getActiveObjects()) elemental->compute(); } diff --git a/framework/src/kokkos/systems/KokkosNonlinearSystemBase.K b/framework/src/kokkos/systems/KokkosNonlinearSystemBase.K index ff1afa436e27..4fa9d276883b 100644 --- a/framework/src/kokkos/systems/KokkosNonlinearSystemBase.K +++ b/framework/src/kokkos/systems/KokkosNonlinearSystemBase.K @@ -112,7 +112,7 @@ NonlinearSystemBase::setKokkosInitialSolution() systems.copyToDevice(); - for (auto nbc : _kokkos_preset_nodal_bcs.getActiveObjects()) + for (auto & nbc : _kokkos_preset_nodal_bcs.getActiveObjects()) std::static_pointer_cast(nbc)->presetSolution(tag); Kokkos::fence(); @@ -168,7 +168,8 @@ NonlinearSystemBase::computeKokkosResidual(const std::set & tags) needed_fe_var_vector_tags); } - // Copy vectors and cache variable values at element quadature points + // Copy solution vectors and residuals and cache variable values at element quadature + // points for (auto & system : systems) { @@ -204,13 +205,13 @@ NonlinearSystemBase::computeKokkosResidual(const std::set & tags) // Compute kernels - for (auto kernel : kernels.getActiveObjects()) + for (auto & kernel : kernels.getActiveObjects()) kernel->computeResidual(); - for (auto nodal_kernel : nodal_kernels.getActiveObjects()) + for (auto & nodal_kernel : nodal_kernels.getActiveObjects()) nodal_kernel->computeResidual(); - for (auto ibc : integrated_bcs.getActiveObjects()) + for (auto & ibc : integrated_bcs.getActiveObjects()) ibc->computeResidual(); } @@ -235,15 +236,15 @@ NonlinearSystemBase::computeKokkosResidual(const std::set & tags) } void -NonlinearSystemBase::computeKokkosNodalBCs(const std::set & tags) +NonlinearSystemBase::computeKokkosNodalBCsResidual(const std::set & tags) { - TIME_SECTION("computeKokkosNodalBC", 1); + TIME_SECTION("computeKokkosNodalBCsResidual", 1); // Get warehouses const auto & nodal_bcs = _kokkos_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0); - if (!nodal_bcs.size()) + if (!nodal_bcs.hasActiveObjects()) return; // Resolve dependencies @@ -259,7 +260,7 @@ NonlinearSystemBase::computeKokkosNodalBCs(const std::set & tags) nodal_bcs.updateFEVariableCoupledVectorTagDependency(needed_fe_var_vector_tags); - // Copy vectors + // Copy solution vectors and residuals for (auto & system : systems) { @@ -278,7 +279,7 @@ NonlinearSystemBase::computeKokkosNodalBCs(const std::set & tags) // Compute kernels - for (auto nbc : nodal_bcs.getActiveObjects()) + for (auto & nbc : nodal_bcs.getActiveObjects()) nbc->computeResidual(); } @@ -311,7 +312,8 @@ NonlinearSystemBase::computeKokkosJacobian(const std::set & tags) const auto & integrated_bcs = _kokkos_integrated_bcs.getMatrixTagsObjectWarehouse(tags, 0); const auto & nodal_bcs = _kokkos_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0); - if (!kernels.size() && !nodal_kernels.size() && !integrated_bcs.size() && !nodal_bcs.size()) + if (!kernels.hasActiveObjects() && !nodal_kernels.hasActiveObjects() && + !integrated_bcs.hasActiveObjects() && !nodal_bcs.hasActiveObjects()) return; // Resolve dependencies @@ -348,7 +350,8 @@ NonlinearSystemBase::computeKokkosJacobian(const std::set & tags) needed_fe_var_vector_tags); } - // Copy vectors and cache variable values at element quadature points + // Copy solution vectors, initialize matrices, and cache variable values at element quadature + // points for (auto & system : systems) { @@ -384,16 +387,16 @@ NonlinearSystemBase::computeKokkosJacobian(const std::set & tags) // Compute kernels - for (auto kernel : kernels.getActiveObjects()) + for (auto & kernel : kernels.getActiveObjects()) kernel->computeJacobian(); - for (auto nodal_kernel : nodal_kernels.getActiveObjects()) + for (auto & nodal_kernel : nodal_kernels.getActiveObjects()) nodal_kernel->computeJacobian(); - for (auto ibc : integrated_bcs.getActiveObjects()) + for (auto & ibc : integrated_bcs.getActiveObjects()) ibc->computeJacobian(); - for (auto nbc : nodal_bcs.getActiveObjects()) + for (auto & nbc : nodal_bcs.getActiveObjects()) nbc->computeJacobian(); } @@ -416,3 +419,120 @@ NonlinearSystemBase::computeKokkosJacobian(const std::set & tags) system.clearActiveSolutionTags(); } } + +void +NonlinearSystemBase::computeKokkosResidualAndJacobian(const std::set & vector_tags, + const std::set & matrix_tags) +{ + TIME_SECTION("computeKokkosResidualAndJacobian", 1); + + if (!_kokkos_kernels.hasActiveObjects() && !_kokkos_nodal_kernels.hasActiveObjects() && + !_kokkos_integrated_bcs.hasActiveObjects() && !_kokkos_nodal_bcs.hasActiveObjects()) + return; + + // Resolve dependencies + + auto & systems = _fe_problem.getKokkosSystems(); + + systems[number()].setActiveResidualTags(vector_tags); + systems[number()].setActiveMatrixTags(matrix_tags); + + std::set needed_moose_vars; + std::set needed_fe_var_vector_tags; + std::unordered_set needed_mat_props; + + for (auto tag : _fe_problem.getVectorTags(Moose::VECTOR_TAG_SOLUTION)) + needed_fe_var_vector_tags.insert(tag._id); + + _kokkos_kernels.updateVariableDependency(needed_moose_vars); + _kokkos_kernels.updateFEVariableCoupledVectorTagDependency(needed_fe_var_vector_tags); + _kokkos_kernels.updateMatPropDependency(needed_mat_props); + + _kokkos_nodal_kernels.updateVariableDependency(needed_moose_vars); + _kokkos_nodal_kernels.updateFEVariableCoupledVectorTagDependency(needed_fe_var_vector_tags); + + _kokkos_integrated_bcs.updateVariableDependency(needed_moose_vars); + _kokkos_integrated_bcs.updateFEVariableCoupledVectorTagDependency(needed_fe_var_vector_tags); + _kokkos_integrated_bcs.updateMatPropDependency(needed_mat_props); + + _kokkos_nodal_bcs.updateVariableDependency(needed_moose_vars); + _kokkos_nodal_bcs.updateFEVariableCoupledVectorTagDependency(needed_fe_var_vector_tags); + + if (needed_mat_props.size()) + { + _fe_problem.getKokkosMaterialsWarehouse().updateVariableDependency(needed_moose_vars); + _fe_problem.getKokkosMaterialsWarehouse().updateFEVariableCoupledVectorTagDependency( + needed_fe_var_vector_tags); + } + + // Copy solution vectors and residuals, initialize matrices, and cache variable values at element + // quadature points + + for (auto & system : systems) + { + system.setActiveVariables(needed_moose_vars); + system.setActiveSolutionTags(needed_fe_var_vector_tags); + + { + TIME_SECTION("KokkosPreallocateAndCopy", 1); + system.sync(Moose::Kokkos::MemcpyType::HOST_TO_DEVICE); + } + { + TIME_SECTION("KokkosReinit", 1); + system.reinit(); + } + } + + systems.copyToDevice(); + + { + TIME_SECTION("KokkosMaterial", 1); + + // Compute material properties + + if (needed_mat_props.size()) + { + _fe_problem.prepareKokkosMaterials(needed_mat_props); + _fe_problem.reinitKokkosMaterials(); + } + } + + { + TIME_SECTION("KokkosKernel", 1); + + // Compute kernels + + for (auto & kernel : _kokkos_kernels.getActiveObjects()) + kernel->computeResidualAndJacobian(); + + for (auto & nodal_kernel : _kokkos_nodal_kernels.getActiveObjects()) + nodal_kernel->computeResidualAndJacobian(); + + for (auto & ibc : _kokkos_integrated_bcs.getActiveObjects()) + ibc->computeResidualAndJacobian(); + + /// Nodal BC residuals are computed separately + for (auto & nbc : _kokkos_nodal_bcs.getActiveObjects()) + nbc->computeJacobian(); + } + + // Close and restore vectors and matrices + + { + TIME_SECTION("KokkosCopy", 1); + + for (auto & system : systems) + system.sync(Moose::Kokkos::MemcpyType::DEVICE_TO_HOST); + } + + // Clear + + systems[number()].clearActiveResidualTags(); + systems[number()].clearActiveMatrixTags(); + + for (auto & system : systems) + { + system.clearActiveVariables(); + system.clearActiveSolutionTags(); + } +} diff --git a/framework/src/systems/NonlinearSystemBase.C b/framework/src/systems/NonlinearSystemBase.C index 1e4ad58587ff..ba1b52e3d11e 100644 --- a/framework/src/systems/NonlinearSystemBase.C +++ b/framework/src/systems/NonlinearSystemBase.C @@ -875,7 +875,7 @@ NonlinearSystemBase::computeResidualTags(const std::set & tags) // We don't want to do nodal bcs or anything else return; - computeNodalBCs(tags); + computeNodalBCsResidual(tags); closeTaggedVectors(tags); // If we are debugging residuals we need one more assignment to have the ghosted copy up to @@ -933,7 +933,7 @@ NonlinearSystemBase::computeResidualAndJacobianTags(const std::set & vect residual.close(); } - computeNodalBCsResidualAndJacobian(); + computeNodalBCsResidualAndJacobian(vector_tags, matrix_tags); closeTaggedVectors(vector_tags); closeTaggedMatrices(matrix_tags); } @@ -1775,13 +1775,6 @@ NonlinearSystemBase::computeResidualInternal(const std::set & tags) residualSetup(); -#ifdef MOOSE_KOKKOS_ENABLED - if (_fe_problem.hasKokkosResidualObjects()) - computeKokkosResidual(tags); -#endif - - const auto vector_tag_data = _fe_problem.getVectorTags(tags); - // Residual contributions from UOs - for now this is used for ray tracing // and ray kernels that contribute to the residual (think line sources) std::vector uos; @@ -1803,6 +1796,11 @@ NonlinearSystemBase::computeResidualInternal(const std::set & tags) for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++) _fe_problem.reinitScalars(tid); +#ifdef MOOSE_KOKKOS_ENABLED + if (_fe_problem.hasKokkosResidualObjects()) + computeKokkosResidual(tags); +#endif + // residual contributions from the domain PARALLEL_TRY { @@ -2045,6 +2043,11 @@ NonlinearSystemBase::computeResidualAndJacobianInternal(const std::set & for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++) _fe_problem.reinitScalars(tid); +#ifdef MOOSE_KOKKOS_ENABLED + if (_fe_problem.hasKokkosResidualObjects()) + computeKokkosResidualAndJacobian(vector_tags, matrix_tags); +#endif + // residual contributions from the domain PARALLEL_TRY { @@ -2087,7 +2090,7 @@ NonlinearSystemBase::computeResidualAndJacobianInternal(const std::set & } void -NonlinearSystemBase::computeNodalBCs(NumericVector & residual) +NonlinearSystemBase::computeNodalBCsResidual(NumericVector & residual) { _nl_vector_tags.clear(); @@ -2096,26 +2099,27 @@ NonlinearSystemBase::computeNodalBCs(NumericVector & residual) _nl_vector_tags.insert(residual_vector_tag._id); associateVectorToTag(residual, residualVectorTag()); - computeNodalBCs(residual, _nl_vector_tags); + computeNodalBCsResidual(residual, _nl_vector_tags); disassociateVectorFromTag(residual, residualVectorTag()); } void -NonlinearSystemBase::computeNodalBCs(NumericVector & residual, const std::set & tags) +NonlinearSystemBase::computeNodalBCsResidual(NumericVector & residual, + const std::set & tags) { associateVectorToTag(residual, residualVectorTag()); - computeNodalBCs(tags); + computeNodalBCsResidual(tags); disassociateVectorFromTag(residual, residualVectorTag()); } void -NonlinearSystemBase::computeNodalBCs(const std::set & tags) +NonlinearSystemBase::computeNodalBCsResidual(const std::set & tags) { #ifdef MOOSE_KOKKOS_ENABLED if (_fe_problem.hasKokkosResidualObjects()) - computeKokkosNodalBCs(tags); + computeKokkosNodalBCsResidual(tags); #endif // We need to close the diag_save_in variables on the aux system before NodalBCBases clear the @@ -2134,7 +2138,7 @@ NonlinearSystemBase::computeNodalBCs(const std::set & tags) nbc_warehouse = &(_nodal_bcs.getVectorTagsObjectWarehouse(tags, 0)); // Return early if there is no nodal kernel - if (!nbc_warehouse->size()) + if (!nbc_warehouse->hasActiveObjects()) return; PARALLEL_TRY @@ -2172,8 +2176,139 @@ NonlinearSystemBase::computeNodalBCs(const std::set & tags) } void -NonlinearSystemBase::computeNodalBCsResidualAndJacobian() +NonlinearSystemBase::computeNodalBCsJacobian(const std::set & tags) +{ + // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs + // on boundary nodes + if (_has_diag_save_in) + _fe_problem.getAuxiliarySystem().solution().close(); + + MooseObjectWarehouse * nbc_warehouse; + + // Select nodal kernels + if (tags.size() == _fe_problem.numMatrixTags() || !tags.size()) + nbc_warehouse = &_nodal_bcs; + else if (tags.size() == 1) + nbc_warehouse = &(_nodal_bcs.getMatrixTagObjectWarehouse(*(tags.begin()), 0)); + else + nbc_warehouse = &(_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0)); + + // Return early if there is no nodal kernel + if (!nbc_warehouse->hasActiveObjects()) + return; + + PARALLEL_TRY + { + // We may be switching from add to set. Moreover, we rely on a call to MatZeroRows to enforce + // the nodal boundary condition constraints which requires that the matrix be truly assembled + // as opposed to just flushed. Consequently we can't do the following despite any desire to + // keep our initial sparsity pattern honored (see https://gitlab.com/petsc/petsc/-/issues/852) + // + // flushTaggedMatrices(tags); + closeTaggedMatrices(tags); + + // Cache the information about which BCs are coupled to which + // variables, so we don't have to figure it out for each node. + std::map> bc_involved_vars; + const std::set & all_boundary_ids = _mesh.getBoundaryIDs(); + for (const auto & bid : all_boundary_ids) + { + // Get reference to all the NodalBCs for this ID. This is only + // safe if there are NodalBCBases there to be gotten... + if (nbc_warehouse->hasActiveBoundaryObjects(bid)) + { + const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid); + for (const auto & bc : bcs) + { + const std::vector & coupled_moose_vars = bc->getCoupledMooseVars(); + + // Create the set of "involved" MOOSE nonlinear vars, which includes all coupled vars + // and the BC's own variable + std::set & var_set = bc_involved_vars[bc->name()]; + for (const auto & coupled_var : coupled_moose_vars) + if (coupled_var->kind() == Moose::VAR_SOLVER) + var_set.insert(coupled_var->number()); + + var_set.insert(bc->variable().number()); + } + } + } + + // reinit scalar variables again. This reinit does not re-fill any of the scalar variable + // solution arrays because that was done above. It only will reorder the derivative + // information for AD calculations to be suitable for NodalBC calculations + for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++) + _fe_problem.reinitScalars(tid, true); + + // Get variable coupling list. We do all the NodalBCBase stuff on + // thread 0... The couplingEntries() data structure determines + // which variables are "coupled" as far as the preconditioner is + // concerned, not what variables a boundary condition specifically + // depends on. + auto & coupling_entries = _fe_problem.couplingEntries(/*_tid=*/0, this->number()); + + // Compute Jacobians for NodalBCBases + const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange(); + for (const auto & bnode : bnd_nodes) + { + BoundaryID boundary_id = bnode->_bnd_id; + Node * node = bnode->_node; + + if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) && + node->processor_id() == processor_id()) + { + _fe_problem.reinitNodeFace(node, boundary_id, 0); + + const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id); + for (const auto & bc : bcs) + { + // Get the set of involved MOOSE vars for this BC + std::set & var_set = bc_involved_vars[bc->name()]; + + // Loop over all the variables whose Jacobian blocks are + // actually being computed, call computeOffDiagJacobian() + // for each one which is actually coupled (otherwise the + // value is zero.) + for (const auto & it : coupling_entries) + { + unsigned int ivar = it.first->number(), jvar = it.second->number(); + + // We are only going to call computeOffDiagJacobian() if: + // 1.) the BC's variable is ivar + // 2.) jvar is "involved" with the BC (including jvar==ivar), and + // 3.) the BC should apply. + if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply()) + bc->computeOffDiagJacobian(jvar); + } + + const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars(); + for (const auto & jvariable : coupled_scalar_vars) + if (hasScalarVariable(jvariable->name())) + bc->computeOffDiagJacobianScalar(jvariable->number()); + } + } + } // end loop over boundary nodes + + // Set the cached NodalBCBase values in the Jacobian matrix + _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{}); + } + PARALLEL_CATCH; +} + +void +NonlinearSystemBase::computeNodalBCsResidualAndJacobian( + [[maybe_unused]] const std::set & vector_tags, + [[maybe_unused]] const std::set & matrix_tags) { +#ifdef MOOSE_KOKKOS_ENABLED + if (_fe_problem.hasKokkosResidualObjects()) + computeKokkosNodalBCsResidual(vector_tags); +#endif + + // Return early if there is no nodal kernel + if (!_nodal_bcs.hasActiveObjects()) + return; + PARALLEL_TRY { const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange(); @@ -2884,11 +3019,6 @@ NonlinearSystemBase::computeJacobianInternal(const std::set & tags) jacobianSetup(); -#ifdef MOOSE_KOKKOS_ENABLED - if (_fe_problem.hasKokkosResidualObjects()) - computeKokkosJacobian(tags); -#endif - // Jacobian contributions from UOs - for now this is used for ray tracing // and ray kernels that contribute to the Jacobian (think line sources) std::vector uos; @@ -2910,6 +3040,11 @@ NonlinearSystemBase::computeJacobianInternal(const std::set & tags) for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++) _fe_problem.reinitScalars(tid); +#ifdef MOOSE_KOKKOS_ENABLED + if (_fe_problem.hasKokkosResidualObjects()) + computeKokkosJacobian(tags); +#endif + PARALLEL_TRY { // We would like to compute ScalarKernels, block NodalKernels, FVFluxKernels, and mortar objects @@ -3097,121 +3232,7 @@ NonlinearSystemBase::computeJacobianInternal(const std::set & tags) } PARALLEL_CATCH; - // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs - // on boundary nodes - if (_has_diag_save_in) - _fe_problem.getAuxiliarySystem().solution().close(); - - PARALLEL_TRY - { - MooseObjectWarehouse * nbc_warehouse; - // Select nodal kernels - if (tags.size() == _fe_problem.numMatrixTags() || !tags.size()) - nbc_warehouse = &_nodal_bcs; - else if (tags.size() == 1) - nbc_warehouse = &(_nodal_bcs.getMatrixTagObjectWarehouse(*(tags.begin()), 0)); - else - nbc_warehouse = &(_nodal_bcs.getMatrixTagsObjectWarehouse(tags, 0)); - - if (nbc_warehouse->hasActiveObjects()) - { - // We may be switching from add to set. Moreover, we rely on a call to MatZeroRows to enforce - // the nodal boundary condition constraints which requires that the matrix be truly assembled - // as opposed to just flushed. Consequently we can't do the following despite any desire to - // keep our initial sparsity pattern honored (see https://gitlab.com/petsc/petsc/-/issues/852) - // - // flushTaggedMatrices(tags); - closeTaggedMatrices(tags); - - // Cache the information about which BCs are coupled to which - // variables, so we don't have to figure it out for each node. - std::map> bc_involved_vars; - const std::set & all_boundary_ids = _mesh.getBoundaryIDs(); - for (const auto & bid : all_boundary_ids) - { - // Get reference to all the NodalBCs for this ID. This is only - // safe if there are NodalBCBases there to be gotten... - if (nbc_warehouse->hasActiveBoundaryObjects(bid)) - { - const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(bid); - for (const auto & bc : bcs) - { - const std::vector & coupled_moose_vars = - bc->getCoupledMooseVars(); - - // Create the set of "involved" MOOSE nonlinear vars, which includes all coupled vars - // and the BC's own variable - std::set & var_set = bc_involved_vars[bc->name()]; - for (const auto & coupled_var : coupled_moose_vars) - if (coupled_var->kind() == Moose::VAR_SOLVER) - var_set.insert(coupled_var->number()); - - var_set.insert(bc->variable().number()); - } - } - } - - // reinit scalar variables again. This reinit does not re-fill any of the scalar variable - // solution arrays because that was done above. It only will reorder the derivative - // information for AD calculations to be suitable for NodalBC calculations - for (unsigned int tid = 0; tid < libMesh::n_threads(); tid++) - _fe_problem.reinitScalars(tid, true); - - // Get variable coupling list. We do all the NodalBCBase stuff on - // thread 0... The couplingEntries() data structure determines - // which variables are "coupled" as far as the preconditioner is - // concerned, not what variables a boundary condition specifically - // depends on. - auto & coupling_entries = _fe_problem.couplingEntries(/*_tid=*/0, this->number()); - - // Compute Jacobians for NodalBCBases - const ConstBndNodeRange & bnd_nodes = _fe_problem.getCurrentAlgebraicBndNodeRange(); - for (const auto & bnode : bnd_nodes) - { - BoundaryID boundary_id = bnode->_bnd_id; - Node * node = bnode->_node; - - if (nbc_warehouse->hasActiveBoundaryObjects(boundary_id) && - node->processor_id() == processor_id()) - { - _fe_problem.reinitNodeFace(node, boundary_id, 0); - - const auto & bcs = nbc_warehouse->getActiveBoundaryObjects(boundary_id); - for (const auto & bc : bcs) - { - // Get the set of involved MOOSE vars for this BC - std::set & var_set = bc_involved_vars[bc->name()]; - - // Loop over all the variables whose Jacobian blocks are - // actually being computed, call computeOffDiagJacobian() - // for each one which is actually coupled (otherwise the - // value is zero.) - for (const auto & it : coupling_entries) - { - unsigned int ivar = it.first->number(), jvar = it.second->number(); - - // We are only going to call computeOffDiagJacobian() if: - // 1.) the BC's variable is ivar - // 2.) jvar is "involved" with the BC (including jvar==ivar), and - // 3.) the BC should apply. - if ((bc->variable().number() == ivar) && var_set.count(jvar) && bc->shouldApply()) - bc->computeOffDiagJacobian(jvar); - } - - const auto & coupled_scalar_vars = bc->getCoupledMooseScalarVars(); - for (const auto & jvariable : coupled_scalar_vars) - if (hasScalarVariable(jvariable->name())) - bc->computeOffDiagJacobianScalar(jvariable->number()); - } - } - } // end loop over boundary nodes - - // Set the cached NodalBCBase values in the Jacobian matrix - _fe_problem.assembly(0, number()).setCachedJacobian(Assembly::GlobalDataKey{}); - } - } - PARALLEL_CATCH; - + computeNodalBCsJacobian(tags); closeTaggedMatrices(tags); // We need to close the save_in variables on the aux system before NodalBCBases clear the dofs diff --git a/modules/doc/content/newsletter/2019/2019_05.md b/modules/doc/content/newsletter/2019/2019_05.md index d8d48dfee014..a55ec0d3f0b7 100644 --- a/modules/doc/content/newsletter/2019/2019_05.md +++ b/modules/doc/content/newsletter/2019/2019_05.md @@ -50,7 +50,7 @@ The need to include the `` template argument for AD objects is no long the MOOSEDocs syntax calls. For example, the following command can be used to list the available parameters for the `ADDiffusion` object. -!listing ADDiffusion.md line=syntax parameters +!listing source/kernels/ADDiffusion.md line=syntax parameters ## What is a Requirement? diff --git a/modules/heat_transfer/src/kokkos/kernels/KokkosHeatConductionTimeDerivative.K b/modules/heat_transfer/src/kokkos/kernels/KokkosHeatConductionTimeDerivative.K index 1217bc09f52b..c48cb9bfa7c3 100644 --- a/modules/heat_transfer/src/kokkos/kernels/KokkosHeatConductionTimeDerivative.K +++ b/modules/heat_transfer/src/kokkos/kernels/KokkosHeatConductionTimeDerivative.K @@ -9,7 +9,7 @@ #include "KokkosHeatConductionTimeDerivative.h" -registerKokkosKernel("HeatTransferApp", KokkosHeatConductionTimeDerivative); +registerKokkosResidualObject("HeatTransferApp", KokkosHeatConductionTimeDerivative); InputParameters KokkosHeatConductionTimeDerivative::validParams() diff --git a/test/include/kokkos/bcs/KokkosCoupledDirichletBC.h b/test/include/kokkos/bcs/KokkosCoupledDirichletBC.h index 1ac35c45a2af..23c1faab1886 100644 --- a/test/include/kokkos/bcs/KokkosCoupledDirichletBC.h +++ b/test/include/kokkos/bcs/KokkosCoupledDirichletBC.h @@ -49,8 +49,10 @@ template KOKKOS_FUNCTION Real KokkosCoupledDirichletBC::computeQpResidual(const unsigned int qp, AssemblyDatum & datum) const { - return _c * _u(datum, qp) + _u(datum, qp) * _u(datum, qp) + _v(datum, qp) * _v(datum, qp) - - _value; + Real u = _u(datum, qp); + Real v = _v(datum, qp); + + return _c * u + u * u + v * v - _value; } template diff --git a/test/tests/kokkos/bcs/ad_coupled_var_neumann/gold/kokkos_ad_coupled_var_neumann_nl_out.e b/test/tests/kokkos/bcs/ad_coupled_var_neumann/gold/kokkos_ad_coupled_var_neumann_nl_out.e new file mode 100644 index 000000000000..fb97c6c7d9db Binary files /dev/null and b/test/tests/kokkos/bcs/ad_coupled_var_neumann/gold/kokkos_ad_coupled_var_neumann_nl_out.e differ diff --git a/test/tests/kokkos/bcs/ad_coupled_var_neumann/gold/kokkos_ad_coupled_var_neumann_out.e b/test/tests/kokkos/bcs/ad_coupled_var_neumann/gold/kokkos_ad_coupled_var_neumann_out.e new file mode 100644 index 000000000000..8bd19633d2d8 Binary files /dev/null and b/test/tests/kokkos/bcs/ad_coupled_var_neumann/gold/kokkos_ad_coupled_var_neumann_out.e differ diff --git a/test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann.i b/test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann.i new file mode 100644 index 000000000000..f68b674e7ff7 --- /dev/null +++ b/test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann.i @@ -0,0 +1,69 @@ +[Mesh] + type = GeneratedMesh + dim = 2 + xmin = 0 + xmax = 1 + ymin = 0 + ymax = 1 + nx = 10 + ny = 10 +[] + +[Variables] + [u] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [diff] + type = KokkosADDiffusion + variable = u + [] +[] + +[AuxVariables] + [coupled_bc_var] + [] +[] + +[ICs] + [coupled_bc_var] + type = FunctionIC + variable = coupled_bc_var + function = set_coupled_bc_var + [] +[] + +[Functions] + [set_coupled_bc_var] + type = ParsedFunction + expression = 'y - 0.5' + [] +[] + +[BCs] + [left] + type = KokkosADDirichletBC + variable = u + boundary = 3 + value = 0 + [] + + [right] + type = KokkosADCoupledVarNeumannBC + variable = u + boundary = 1 + v = coupled_bc_var + [] +[] + +[Executioner] + type = Steady + residual_and_jacobian_together = true +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann_nl.i b/test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann_nl.i new file mode 100644 index 000000000000..51df5d747d1b --- /dev/null +++ b/test/tests/kokkos/bcs/ad_coupled_var_neumann/kokkos_ad_coupled_var_neumann_nl.i @@ -0,0 +1,64 @@ +[Mesh] + type = GeneratedMesh + dim = 1 + nx = 20 +[] + +[Variables] + [u][] + [v][] +[] + +[Kernels] + [diff] + type = KokkosADDiffusion + variable = u + [] + [diff_v] + type = KokkosADDiffusion + variable = v + [] +[] + +[BCs] + [left] + type = KokkosADDirichletBC + variable = u + boundary = 'left' + value = 0 + [] + [right] + type = KokkosADCoupledVarNeumannBC + variable = u + boundary = 'right' + v = v + [] + [v_left] + type = KokkosADDirichletBC + variable = v + boundary = 'left' + value = 0 + [] + [v_right] + type = KokkosADDirichletBC + variable = v + boundary = 'right' + value = 1 + [] +[] + +[Preconditioning] + [smp] + type = SMP + full = true + [] +[] + +[Executioner] + type = Steady + residual_and_jacobian_together = true +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/kokkos/bcs/ad_coupled_var_neumann/tests b/test/tests/kokkos/bcs/ad_coupled_var_neumann/tests new file mode 100644 index 000000000000..746d4c7642fb --- /dev/null +++ b/test/tests/kokkos/bcs/ad_coupled_var_neumann/tests @@ -0,0 +1,34 @@ +[Tests] + issues = '#30655' + design = 'KokkosADCoupledVarNeumannBC.md' + + [test] + type = 'Exodiff' + input = 'kokkos_ad_coupled_var_neumann.i' + exodiff = 'kokkos_ad_coupled_var_neumann_out.e' + requirement = 'The Kokkos system shall support coupled Neumann type boundary conditions using automatic differentiation.' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + + [nonlinear] + requirement = 'When coupling nonlinear variables into a Neumann type boundary condition using automatic differentiation, the Kokkos system shall' + [exo] + type = Exodiff + input = 'kokkos_ad_coupled_var_neumann_nl.i' + exodiff = 'kokkos_ad_coupled_var_neumann_nl_out.e' + detail = 'generate accurate results' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + [jac] + type = PetscJacobianTester + run_sim = True + input = 'kokkos_ad_coupled_var_neumann_nl.i' + difference_tol = 1e-6 + detail = 'generate a perfect Jacobian' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + [] +[] diff --git a/test/tests/kokkos/kernels/2d_diffusion/kokkos_2d_diffusion_test.i b/test/tests/kokkos/kernels/2d_diffusion/kokkos_2d_diffusion_test.i index 4c4fd888434d..88dc89dde162 100644 --- a/test/tests/kokkos/kernels/2d_diffusion/kokkos_2d_diffusion_test.i +++ b/test/tests/kokkos/kernels/2d_diffusion/kokkos_2d_diffusion_test.i @@ -17,8 +17,6 @@ [] [Variables] - active = 'u' - [u] order = FIRST family = LAGRANGE @@ -26,8 +24,6 @@ [] [Kernels] - active = 'diff' - [diff] type = KokkosDiffusion variable = u @@ -35,9 +31,6 @@ [] [BCs] - # BCs cannot be preset due to Jacobian test - active = 'left right' - [left] type = KokkosDirichletBC variable = u diff --git a/test/tests/kokkos/kernels/ad_diffusion/gold/kokkos_2d_ad_diffusion_neumannbc_test_out.e b/test/tests/kokkos/kernels/ad_diffusion/gold/kokkos_2d_ad_diffusion_neumannbc_test_out.e new file mode 100644 index 000000000000..f20f4a3cfad7 Binary files /dev/null and b/test/tests/kokkos/kernels/ad_diffusion/gold/kokkos_2d_ad_diffusion_neumannbc_test_out.e differ diff --git a/test/tests/kokkos/kernels/ad_diffusion/gold/kokkos_2d_ad_diffusion_test_out.e b/test/tests/kokkos/kernels/ad_diffusion/gold/kokkos_2d_ad_diffusion_test_out.e new file mode 100644 index 000000000000..20c28eedbc73 Binary files /dev/null and b/test/tests/kokkos/kernels/ad_diffusion/gold/kokkos_2d_ad_diffusion_test_out.e differ diff --git a/test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_neumannbc_test.i b/test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_neumannbc_test.i new file mode 100644 index 000000000000..e75377b3d6e7 --- /dev/null +++ b/test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_neumannbc_test.i @@ -0,0 +1,48 @@ +[Mesh] + [square] + type = GeneratedMeshGenerator + nx = 2 + ny = 2 + dim = 2 + [] +[] + +[Variables] + [u] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [diff] + type = KokkosADDiffusion + variable = u + [] +[] + +[BCs] + [left] + type = KokkosADDirichletBC + variable = u + boundary = 3 + value = 0 + [] + [right] + type = KokkosADNeumannBC + variable = u + boundary = 1 + value = 1 + [] +[] + +[Executioner] + type = Steady + + solve_type = NEWTON + residual_and_jacobian_together = true +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_test.i b/test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_test.i new file mode 100644 index 000000000000..2b2314f9fae5 --- /dev/null +++ b/test/tests/kokkos/kernels/ad_diffusion/kokkos_2d_ad_diffusion_test.i @@ -0,0 +1,48 @@ +[Mesh] + [square] + type = GeneratedMeshGenerator + nx = 2 + ny = 2 + dim = 2 + [] +[] + +[Variables] + [u] + order = FIRST + family = LAGRANGE + [] +[] + +[Kernels] + [diff] + type = KokkosADDiffusion + variable = u + [] +[] + +[BCs] + [left] + type = KokkosADDirichletBC + variable = u + boundary = 3 + value = 0 + [] + [right] + type = KokkosADDirichletBC + variable = u + boundary = 1 + value = 1 + [] +[] + +[Executioner] + type = Steady + + solve_type = NEWTON + residual_and_jacobian_together = true +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/kokkos/kernels/ad_diffusion/tests b/test/tests/kokkos/kernels/ad_diffusion/tests new file mode 100644 index 000000000000..0492fb5d358f --- /dev/null +++ b/test/tests/kokkos/kernels/ad_diffusion/tests @@ -0,0 +1,32 @@ +[Tests] + issues = '#32699' + + [dirichlet] + type = Exodiff + input = 'kokkos_2d_ad_diffusion_test.i' + exodiff = 'kokkos_2d_ad_diffusion_test_out.e' + requirement = 'The Kokkos system shall provide an ability to solve a 2D diffusion problem with Dirichlet boundary conditions using automatic differentation.' + design = 'KokkosADDiffusion.md KokkosADDirichletBC.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + [neumann] + type = Exodiff + input = 'kokkos_2d_ad_diffusion_neumannbc_test.i' + exodiff = 'kokkos_2d_ad_diffusion_neumannbc_test_out.e' + requirement = 'The Kokkos system shall provide an ability to solve a 2D diffusion problem with Neumann boundary conditions using automatic differentation.' + design = 'KokkosADDiffusion.md KokkosADNeumannBC.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + [jacobian] + type = PetscJacobianTester + input = 'kokkos_2d_ad_diffusion_test.i' + cli_args = 'Outputs/exodus=false' + recover = false + requirement = 'The Kokkos system shall provide an exact Jacobian using automatic differentation.' + design = 'framework_stp.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] +[] diff --git a/test/tests/kokkos/kernels/ad_time_derivative/gold/kokkos_ad_coupled_time_derivative_test_out.e b/test/tests/kokkos/kernels/ad_time_derivative/gold/kokkos_ad_coupled_time_derivative_test_out.e new file mode 100644 index 000000000000..4a7538b92f82 Binary files /dev/null and b/test/tests/kokkos/kernels/ad_time_derivative/gold/kokkos_ad_coupled_time_derivative_test_out.e differ diff --git a/test/tests/kokkos/kernels/ad_time_derivative/gold/kokkos_ad_time_derivative_test_out.e b/test/tests/kokkos/kernels/ad_time_derivative/gold/kokkos_ad_time_derivative_test_out.e new file mode 100644 index 000000000000..20f8c257834e Binary files /dev/null and b/test/tests/kokkos/kernels/ad_time_derivative/gold/kokkos_ad_time_derivative_test_out.e differ diff --git a/test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_coupled_time_derivative_test.i b/test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_coupled_time_derivative_test.i new file mode 100644 index 000000000000..8aa019c56215 --- /dev/null +++ b/test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_coupled_time_derivative_test.i @@ -0,0 +1,60 @@ +[Mesh] + type = GeneratedMesh + nx = 5 + ny = 5 + dim = 2 +[] + +[Variables] + [u] + [] + [v] + [] +[] + +[Kernels] + [time_u] + type = KokkosADTimeDerivative + variable = u + [] + [fn_u] + type = KokkosBodyForce + variable = u + [] + [time_v] + type = KokkosADCoupledTimeDerivative + variable = v + v = u + [] + [diff_v] + type = KokkosADDiffusion + variable = v + [] +[] + +[BCs] + [left] + type = KokkosADDirichletBC + variable = v + boundary = 'left' + value = 0 + [] + [right] + type = KokkosADDirichletBC + variable = v + boundary = 'right' + value = 1 + [] +[] + +[Executioner] + type = Transient + num_steps = 1 + end_time = 10 + solve_type = 'NEWTON' + residual_and_jacobian_together = true +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_time_derivative_test.i b/test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_time_derivative_test.i new file mode 100644 index 000000000000..31ba10071675 --- /dev/null +++ b/test/tests/kokkos/kernels/ad_time_derivative/kokkos_ad_time_derivative_test.i @@ -0,0 +1,51 @@ +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 10 + ny = 10 +[] + +[Variables] + [u] + [] +[] + +[Kernels] + [td] + type = KokkosADTimeDerivative + variable = u + [] + [diff] + type = KokkosADDiffusion + variable = u + [] +[] + +[BCs] + [left] + type = KokkosADDirichletBC + variable = u + boundary = left + value = 0 + [] + [right] + type = KokkosADDirichletBC + variable = u + boundary = right + value = 1 + [] +[] + +[Executioner] + type = Transient + num_steps = 5 + dt = 1e-2 + solve_type = NEWTON + petsc_options_iname = '-pc_type -pc_hypre_type' + petsc_options_value = 'hypre boomeramg' + residual_and_jacobian_together = true +[] + +[Outputs] + exodus = true +[] diff --git a/test/tests/kokkos/kernels/ad_time_derivative/tests b/test/tests/kokkos/kernels/ad_time_derivative/tests new file mode 100644 index 000000000000..5b512c5f8137 --- /dev/null +++ b/test/tests/kokkos/kernels/ad_time_derivative/tests @@ -0,0 +1,58 @@ +[Tests] + issues = '#32699' + + [test] + type = 'Exodiff' + input = 'kokkos_ad_time_derivative_test.i' + exodiff = 'kokkos_ad_time_derivative_test_out.e' + requirement = 'The Kokkos system shall provide a time derivative kernel using automatic differentation.' + design = 'KokkosADTimeDerivative.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + [test_jacobian] + type = 'PetscJacobianTester' + input = 'kokkos_ad_time_derivative_test.i' + cli_args = 'Executioner/num_steps=1 Outputs/exodus=false' + run_sim = True + ratio_tol = 1e-7 + difference_tol = 1e-6 + requirement = 'The Jacobian from a Kokkos-based time derivative kernel shall be perfect.' + design = 'KokkosADTimeDerivative.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + + [coupled] + type = 'Exodiff' + input = 'kokkos_ad_coupled_time_derivative_test.i' + exodiff = 'kokkos_ad_coupled_time_derivative_test_out.e' + requirement = 'The Kokkos system shall provide a coupled time derivative kernel using automatic differentation.' + design = 'KokkosADCoupledTimeDerivative.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + [coupled_jacobian] + type = 'PetscJacobianTester' + input = 'kokkos_ad_coupled_time_derivative_test.i' + cli_args = 'Outputs/exodus=false' + run_sim = True + ratio_tol = 1e-7 + difference_tol = 1e-6 + requirement = 'The Jacobian from a Kokkos-based coupled variable time derivative kernel shall be perfect.' + design = 'KokkosADCoupledTimeDerivative.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] + + [separate] + type = 'Exodiff' + input = 'kokkos_ad_coupled_time_derivative_test.i' + exodiff = 'kokkos_ad_coupled_time_derivative_test_out.e' + cli_args = 'Executioner/residual_and_jacobian_together=false' + requirement = 'The Kokkos system shall be able to compute residual and Jacobian separately with AD kernels.' + design = 'Executioner.md' + capabilities = 'kokkos' + compute_devices = 'cpu cuda' + [] +[] diff --git a/tutorials/darcy_thermo_mech/doc/content/workshop/systems/kernels.md b/tutorials/darcy_thermo_mech/doc/content/workshop/systems/kernels.md index 286c7fc246d0..3aecabdb6f8b 100644 --- a/tutorials/darcy_thermo_mech/doc/content/workshop/systems/kernels.md +++ b/tutorials/darcy_thermo_mech/doc/content/workshop/systems/kernels.md @@ -65,7 +65,7 @@ where $\psi_i$ are the test functions and $u_h$ is the finite element solution. ## ADDiffusion.h -!listing ADDiffusion.h +!listing include/kernels/ADDiffusion.h !---