From b3cbe435928e96c23f31824dba4418b0b8f70ce0 Mon Sep 17 00:00:00 2001 From: ttruster Date: Thu, 10 Nov 2022 23:39:07 -0500 Subject: [PATCH] Add separate assembly case for two scalar variables Add a bunch more comments to document which terms are done where Verified with test case --- framework/include/kernels/KernelScalarBase.h | 2 +- framework/src/kernels/KernelScalarBase.C | 13 +- ...ogenizedTotalLagrangianStressDivergenceA.h | 156 ++++++ ...ogenizedTotalLagrangianStressDivergenceS.h | 20 +- .../ComputeHomogenizedLagrangianStrainA.h | 61 +++ ...ogenizedTotalLagrangianStressDivergenceA.C | 477 ++++++++++++++++++ ...ogenizedTotalLagrangianStressDivergenceS.C | 56 +- .../ComputeHomogenizedLagrangianStrainA.C | 87 ++++ .../2dscalar.i => scalar_kernel/2dsole.i} | 49 +- .../scalar_kernel/gold/2dsole_out.csv | 7 + .../total/homogenization/scalar_kernel/tests | 15 +- 11 files changed, 884 insertions(+), 59 deletions(-) create mode 100644 modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceA.h create mode 100644 modules/tensor_mechanics/test/include/materials/ComputeHomogenizedLagrangianStrainA.h create mode 100644 modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceA.C create mode 100644 modules/tensor_mechanics/test/src/materials/ComputeHomogenizedLagrangianStrainA.C rename modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/{small-tests/2dscalar.i => scalar_kernel/2dsole.i} (82%) create mode 100644 modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/gold/2dsole_out.csv diff --git a/framework/include/kernels/KernelScalarBase.h b/framework/include/kernels/KernelScalarBase.h index 2d0459124d88..75706aee9b3e 100644 --- a/framework/include/kernels/KernelScalarBase.h +++ b/framework/include/kernels/KernelScalarBase.h @@ -43,7 +43,7 @@ class KernelScalarBase : public Kernel virtual void computeOffDiagJacobian(unsigned int jvar_num) override; /** * Computes jacobian block with respect to a scalar variable - * @param jvar, the number of the (other) scalar variable + * @param svar, the number of the (other) scalar variable */ void computeOffDiagJacobianScalar(unsigned int svar_num) override; diff --git a/framework/src/kernels/KernelScalarBase.C b/framework/src/kernels/KernelScalarBase.C index 38348a9704c0..e846017293b2 100644 --- a/framework/src/kernels/KernelScalarBase.C +++ b/framework/src/kernels/KernelScalarBase.C @@ -202,20 +202,23 @@ KernelScalarBase::computeOffDiagJacobianScalar(const unsigned int svar_num) } else if (svar_num == _kappa_var) // column for this kernel's scalar variable { - // Perform assembly using method in Kernel + // Perform assembly using method in Kernel; works for simple cases but not general // Kernel::computeOffDiagJacobianScalar(svar_num); // d-_var-residual / d-_kappa - // Perform assembly using DenseMatrix like d-_kappa_var-residual / d-_var + // Perform assembly using local_ke like d-_kappa_var-residual / d-_var computeOffDiagJacobianScalarLocal(svar_num); // d-_var-residual / d-_kappa computeScalarJacobian(); // d-_kappa-residual / d-_kappa } else // some other column for scalar variable { - Kernel::computeOffDiagJacobianScalar(svar_num); // d-_var-residual / d-jvar - computeScalarOffDiagJacobianScalar(svar_num); // d-_kappa-residual / d-jvar + // Perform assembly using method in Kernel; works for simple cases but not general + // Kernel::computeOffDiagJacobianScalar(svar_num); // d-_var-residual / d-jvar + // Perform assembly using local_ke like d-_kappa_var-residual / d-_var + computeOffDiagJacobianScalarLocal(svar_num); // d-_var-residual / d-svar + computeScalarOffDiagJacobianScalar(svar_num); // d-_kappa-residual / d-svar } } else - Kernel::computeOffDiagJacobianScalar(svar_num); // d-_var-residual / d-jvar + Kernel::computeOffDiagJacobianScalar(svar_num); // d-_var-residual / d-svar } void diff --git a/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceA.h b/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceA.h new file mode 100644 index 000000000000..db197f018e2d --- /dev/null +++ b/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceA.h @@ -0,0 +1,156 @@ +//* 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 "TotalLagrangianStressDivergenceS.h" + +// Helpers common to the whole homogenization system +namespace HomogenizationA +{ +/// Moose constraint type, for input +const MultiMooseEnum constraintType("strain stress none"); +/// Constraint type: stress/PK stress or strain/deformation gradient +enum class ConstraintType +{ + Strain, + Stress, + None +}; +typedef std::map, std::pair> + ConstraintMap; +} + +/// Total Lagrangian formulation with most homogenization terms (one disp_xyz field and one scalar) +/// The macro_gradient variable is split into two scalars: the first component called '_hvar' +/// herein and all other components called '_avar' herein. For parameter _beta = 0, the primary +/// scalar (_kappa) is _hvar and the coupled scalar is _avar. For parameter _beta = 1, the primary +/// scalar (_kappa) is _avar and the coupled scalar is _hvar. Just like the primary field variable +/// (_var) is either disp_x or disp_y or disp_z depending on _alpha. +/// +/// Thus, each instance of HomogenizedTotalLagrangianStressDivergenceA acts on one field variable +/// (_disp_alpha) and one scalar variable (_hvar_beta). The job of the kernel is to assemble the +/// residual of all dofs of _disp_alpha and of all dofs of _hvar_beta (namely, selected rows). +/// Similarly, it assembles the ENTIRE row for _disp_alpha and _hvar_beta (namely the columns +/// from all dofs of all _disp field variables and all dofs of all scalar variables _hvar and +/// _avar). The rows for the other field/scalar variables are handled by other instances of the +/// kernel, which have other values of _alpha AND _beta. The logical checks ensure the proper +/// decomposition of the jobs. +/// +/// In summary, for x=disp_x etc. and h=_hvar and a=_avar, then the contributions of the instances are +/// _alpha=0, _beta=0 +/// R = [Rx, 00, 00, Rh, 00 ]^T +/// J = [Jxx, Jxy, Jxz, Jxh, 000 +/// Jhx, 000, 000, Jhh, Jha] +/// _alpha=1, _beta=0 +/// R = [00, Ry, 00, 00, 00 ]^T +/// J = [Jyx, Jyy, Jyz, Jyh, 000 +/// 000, Jhy, 000, 000, 000] +/// _alpha=2, _beta=0 +/// R = [00, 00, Rz, 00, 00 ]^T +/// J = [Jzx, Jzy, Jzz, Jzh, 000 +/// 000, 000, Jhz, 000, 000] +/// _alpha=0, _beta=1 +/// R = [00, 00, 00, 00, Ra ]^T +/// J = [000, 000, 000, 000, Jxa +/// Jax, 000, 000, Jah, Jaa] +/// _alpha=1, _beta=1 +/// R = [00, 00, 00, 00, 00 ]^T +/// J = [000, 000, 000, 000, Jya +/// 000, Jay, 000, 000, 000] +/// _alpha=2, _beta=1 +/// R = [00, 00, 00, 00, 00 ]^T +/// J = [000, 000, 000, 000, Jza +/// 000, 000, Jaz, 000, 000] +/// +/// In this manner, the full R and J are obtained with NO duplication of jobs: +/// R = [Rx, Ry, Rz, Rh, Ra ]^T +/// J = [Jxx, Jxy, Jxz, Jxh, Jxa +/// Jxy, Jyy, Jyz, Jyh, Jya +/// Jzx, Jzy, Jzz, Jzh, Jza +/// Jhx, Jhy, Jhz, Jhh, Jha +/// Jax, Jay, Jaz, Jah, Jaa] +/// +class HomogenizedTotalLagrangianStressDivergenceA : public TotalLagrangianStressDivergenceS +{ +public: + static InputParameters validParams(); + HomogenizedTotalLagrangianStressDivergenceA(const InputParameters & parameters); + +protected: + // Add overrides to base class contributions to only happen for _beta==0, to happen only once + virtual Real computeQpResidual() override; + virtual Real computeQpJacobianDisplacement(unsigned int alpha, unsigned int beta) override; + + /** + * Method for computing the scalar part of residual for _kappa + */ + virtual void computeScalarResidual() override; + + /** + * Method for computing the scalar variable part of Jacobian for d-_kappa-residual / d-_kappa + */ + virtual void computeScalarJacobian() override; + + /** + * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-jvar + * jvar is looped over all field variables, which herein is just disp_x and disp_y + */ + virtual void computeScalarOffDiagJacobian(const unsigned int jvar_num) override; + + /** + * Method for computing an off-diagonal jacobian component at quadrature points. + */ + virtual Real computeScalarQpOffDiagJacobian(const unsigned int jvar_num) override; + + /** + * Method for computing an off-diagonal jacobian component d-_var-residual / d-svar. + * svar is looped over all scalar variables, which herein is just _kappa and _kappa_other + */ + virtual void computeOffDiagJacobianScalarLocal(const unsigned int svar_num) override; + + /** + * Method for computing d-_var-residual / d-_svar at quadrature points. + */ + virtual Real computeQpOffDiagJacobianScalar(const unsigned int jvar) override; + + /** + * Method for computing an off-diagonal jacobian component d-_kappa-residual / d-svar + * svar is looped over other scalar variables, which herein is just _kappa_other + */ + virtual void computeScalarOffDiagJacobianScalar(const unsigned int svar_num) override; + +protected: + /// Which component of the scalar vector residual this constraint is responsible for + const unsigned int _beta; + + /// (Pointer to) Scalar variable this kernel operates on + const MooseVariableScalar * const _kappao_var_ptr; + + /// The unknown scalar variable ID + const unsigned int _kappao_var; + + /// Order of the scalar variable, used in several places + const unsigned int _ko_order; + + /// Reference to the current solution at the current quadrature point + const VariableValue & _kappa_other; + + /// Type of each constraint (stress or strain) for each component + HomogenizationA::ConstraintMap _cmap; + + /// The constraint type; initialize with 'none' + HomogenizationA::ConstraintType _ctype = HomogenizationA::ConstraintType::None; + + /// Used internally to iterate over each scalar component + unsigned int _m; + unsigned int _n; + unsigned int _a; + unsigned int _b; +}; diff --git a/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceS.h b/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceS.h index b71ae5cf7323..582bb0cf99ac 100644 --- a/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceS.h +++ b/modules/tensor_mechanics/test/include/kernels/HomogenizedTotalLagrangianStressDivergenceS.h @@ -27,12 +27,7 @@ typedef std::map, std::pair, std::pair> + ConstraintMap; +} + +/// Calculate the tensor corresponding to homogenization gradient +/// +/// This class takes TWO scalar fields of the correct total size, representing a +/// constant deformation gradient over the domain, and casts them into +/// a RankTwo material property. _macro_gradientA is the 1st component and +/// _macro_gradient is the rest of the components. +/// +class ComputeHomogenizedLagrangianStrainA : public Material +{ +public: + static InputParameters validParams(); + ComputeHomogenizedLagrangianStrainA(const InputParameters & parameters); + +protected: + virtual void computeQpProperties() override; + +protected: + /// The base name for material properties + const std::string _base_name; + + /// Constraint map + HomogenizationB::ConstraintMap _cmap; + + /// ScalarVariable with the field + const VariableValue & _macro_gradient; + + /// ScalarVariable with 1st component of the field + const VariableValue & _macro_gradientA; + + /// Unwrapped into a tensor + MaterialProperty & _homogenization_contribution; +}; diff --git a/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceA.C b/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceA.C new file mode 100644 index 000000000000..a5e64d6f4e3e --- /dev/null +++ b/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceA.C @@ -0,0 +1,477 @@ +//* 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 "HomogenizedTotalLagrangianStressDivergenceA.h" + +// MOOSE includes +#include "Function.h" +#include "MooseVariableScalar.h" +// #include "Assembly.h" +// #include "MooseVariableFE.h" +// #include "MooseVariableScalar.h" +// #include "SystemBase.h" + +// #include "libmesh/quadrature.h" + +registerMooseObject("TensorMechanicsTestApp", HomogenizedTotalLagrangianStressDivergenceA); + +namespace +{ +const InputParameters & +setHTLSDAParam(const InputParameters & params_in) +{ + // Reset the scalar_variable parameter to a relevant name for this physics + InputParameters & ret = const_cast(params_in); + ret.set("scalar_variable") = {params_in.get("macro_var")}; + return ret; +} +} + +InputParameters +HomogenizedTotalLagrangianStressDivergenceA::validParams() +{ + InputParameters params = TotalLagrangianStressDivergenceS::validParams(); + params.addClassDescription("Total Lagrangian stress equilibrium kernel with " + "homogenization constraint Jacobian terms"); + params.addRequiredParam("macro_var", + "Optional scalar field with the macro gradient"); + params.addRequiredCoupledVar("macro_other", "Other components of coupled scalar variable"); + params.addRequiredParam("prime_scalar", "Either 0=_var or 1=_other scalar"); + params.addRequiredParam( + "constraint_types", + HomogenizationA::constraintType, + "Type of each constraint: strain, stress, or none. The types are specified in the " + "column-major order, and there must be 9 entries in total."); + params.addRequiredParam>( + "targets", "Functions giving the targets to hit for constraint types that are not none."); + + return params; +} + +HomogenizedTotalLagrangianStressDivergenceA::HomogenizedTotalLagrangianStressDivergenceA( + const InputParameters & parameters) + : TotalLagrangianStressDivergenceS(setHTLSDAParam(parameters)), + _beta(getParam("prime_scalar")), + _kappao_var_ptr(getScalarVar("macro_other", 0)), + _kappao_var(coupledScalar("macro_other")), + _ko_order(getScalarVar("macro_other", 0)->order()), + _kappa_other(coupledScalarValue("macro_other")) +{ + // Constraint types + auto types = getParam("constraint_types"); + if (types.size() != Moose::dim * Moose::dim) + mooseError("Number of constraint types must equal dim * dim. ", types.size(), " are provided."); + + // Targets to hit + const std::vector & fnames = getParam>("targets"); + + // Prepare the constraint map + unsigned int fcount = 0; + for (const auto j : make_range(Moose::dim)) + for (const auto i : make_range(Moose::dim)) + { + const auto idx = i + Moose::dim * j; + const auto ctype = static_cast(types.get(idx)); + if (ctype != HomogenizationA::ConstraintType::None) + { + const Function * const f = &getFunctionByName(fnames[fcount++]); + _cmap[{i, j}] = {ctype, f}; + } + } +} + +Real +HomogenizedTotalLagrangianStressDivergenceA::computeQpResidual() +{ + // Assemble R_alpha if _beta==0 + if (_beta == 0) + return gradTest(_alpha).doubleContraction(_pk1[_qp]); + else + return 0.0; +} + +Real +HomogenizedTotalLagrangianStressDivergenceA::computeQpJacobianDisplacement(unsigned int alpha, + unsigned int beta) +{ + // Assemble J-alpha-beta if _beta==0 + if (_beta == 0) + return gradTest(alpha).doubleContraction(_dpk1[_qp] * gradTrial(beta)); + else + return 0.0; +} + +void +HomogenizedTotalLagrangianStressDivergenceA::computeScalarResidual() +{ + std::vector scalar_residuals(_k_order); + + // Assemble R_beta if _alpha==0 + if (_alpha == 0) + { + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + initScalarQpResidual(); + Real dV = _JxW[_qp] * _coord[_qp]; + _h = 0; // single index for residual vector; double indices for constraint tensor component + for (auto && [indices, constraint] : _cmap) + { + auto && [i, j] = indices; + auto && [ctype, ctarget] = constraint; + + // ONLY the component(s) that this constraint will contribute here; + // other one is handled in the other constraint + if (_beta == 0) + { + if (_h == 1) // only assemble first=0 component of _hvar, then break the loop + break; + } + else + { + // skip the first component (_hvar) and continue to "first" component of _avar + if (_h == 0) + { + _h++; + continue; + } + } + + // I am not great with C++ precedence; so, store the index + unsigned int r_ind = -_beta + _h; // move 1 row up if _beta=1 for the other scalar + _h++; // increment the index before we forget + if (_large_kinematics) + { + if (ctype == HomogenizationA::ConstraintType::Stress) + scalar_residuals[r_ind] += dV * (_pk1[_qp](i, j) - ctarget->value(_t, _q_point[_qp])); + else if (ctype == HomogenizationA::ConstraintType::Strain) + scalar_residuals[r_ind] += + dV * (_F[_qp](i, j) - (Real(i == j) + ctarget->value(_t, _q_point[_qp]))); + else + mooseError("Unknown constraint type in the integral!"); + } + else + { + if (ctype == HomogenizationA::ConstraintType::Stress) + scalar_residuals[r_ind] += dV * (_pk1[_qp](i, j) - ctarget->value(_t, _q_point[_qp])); + else if (ctype == HomogenizationA::ConstraintType::Strain) + scalar_residuals[r_ind] += dV * (0.5 * (_F[_qp](i, j) + _F[_qp](j, i)) - + (Real(i == j) + ctarget->value(_t, _q_point[_qp]))); + else + mooseError("Unknown constraint type in the integral!"); + } + } + } + } + + _assembly.processResiduals(scalar_residuals, + _kappa_var_ptr->dofIndices(), + _vector_tags, + _kappa_var_ptr->scalingFactor()); +} + +void +HomogenizedTotalLagrangianStressDivergenceA::computeScalarJacobian() +{ + _local_ke.resize(_k_order, _k_order); + + // Assemble J_beta_beta if _alpha==0 + if (_alpha == 0) + { + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + initScalarQpJacobian(_kappa_var); + Real dV = _JxW[_qp] * _coord[_qp]; + + _h = 0; + for (auto && [indices1, constraint1] : _cmap) + { + auto && [i, j] = indices1; + auto && ctype = constraint1.first; + + // identical logic to computeScalarResidual + if (_beta == 0) + { + if (_h == 1) + break; + } + else + { + if (_h == 0) + { + _h++; + continue; + } + } + + _l = 0; + for (auto && indices2 : _cmap) + { + auto && [a, b] = indices2.first; + + // identical logic to computeScalarResidual, but for column index + if (_beta == 0) + { + if (_l == 1) + break; + } + else + { + if (_l == 0) + { + _l++; + continue; + } + } + + unsigned int c_ind = -_beta + _l; // move 1 column left if _beta=1 for the other scalar + _l++; // increment the index before we forget + if (ctype == HomogenizationA::ConstraintType::Stress) + _local_ke(-_beta + _h, c_ind) += dV * (_dpk1[_qp](i, j, a, b)); + else if (ctype == HomogenizationA::ConstraintType::Strain) + { + if (_large_kinematics) + _local_ke(-_beta + _h, c_ind) += dV * (Real(i == a && j == b)); + else + _local_ke(-_beta + _h, c_ind) += + dV * (0.5 * Real(i == a && j == b) + 0.5 * Real(i == b && j == a)); + } + else + mooseError("Unknown constraint type in Jacobian calculator!"); + } + _h++; + } + } + } + + for (const auto & matrix_tag : _matrix_tags) + _assembly.cacheJacobianBlock(_local_ke, + _kappa_var_ptr->dofIndices(), + _kappa_var_ptr->dofIndices(), + _kappa_var_ptr->scalingFactor(), + matrix_tag); +} + +void +HomogenizedTotalLagrangianStressDivergenceA::computeScalarOffDiagJacobian( + const unsigned int jvar_num) +{ + // Assembly J_alpha_beta ONLY + if (jvar_num == _var.number()) + { + _local_ke.resize(_k_order, _test.size()); + + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + // single index for Jacobian column; double indices for constraint tensor component + unsigned int h = 0; + Real dV = _JxW[_qp] * _coord[_qp]; + for (auto && [indices, constraint] : _cmap) + { + // identical logic to computeScalarResidual + if (_beta == 0) + { + if (h == 1) + break; + } + else + { + if (h == 0) + { + h++; + continue; + } + } + // copy constraint indices to protected variables to pass to Qp routine + _m = indices.first; + _n = indices.second; + _ctype = constraint.first; + initScalarQpOffDiagJacobian(_var); + for (_j = 0; _j < _test.size(); _j++) + _local_ke(-_beta + h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num); + h++; + } + } + + for (const auto & matrix_tag : _matrix_tags) + _assembly.cacheJacobianBlock(_local_ke, + _kappa_var_ptr->dofIndices(), + _var.dofIndices(), + _kappa_var_ptr->scalingFactor(), + matrix_tag); + } +} + +void +HomogenizedTotalLagrangianStressDivergenceA::computeOffDiagJacobianScalarLocal( + const unsigned int svar_num) +{ + // Assembly J_beta_alpha ONLY + if (svar_num == _kappa_var) + { + _local_ke.resize(_test.size(), _k_order); + + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + // single index for Jacobian row; double indices for constraint tensor component + unsigned int l = 0; + Real dV = _JxW[_qp] * _coord[_qp]; + for (auto && [indices, constraint] : _cmap) + { + // identical logic to computeScalarResidual, but for column index + if (_beta == 0) + { + if (l == 1) + break; + } + else + { + if (l == 0) + { + l++; + continue; + } + } + // copy constraint indices to protected variables to pass to Qp routine + _m = indices.first; + _n = indices.second; + _ctype = constraint.first; + initScalarQpJacobian(svar_num); + for (_i = 0; _i < _test.size(); _i++) + { + _local_ke(_i, -_beta + l) += dV * computeQpOffDiagJacobianScalar(svar_num); + } + l++; + } + } + + for (const auto & matrix_tag : _matrix_tags) + _assembly.cacheJacobianBlock(_local_ke, + _var.dofIndices(), + _kappa_var_ptr->dofIndices(), + _var.scalingFactor(), + matrix_tag); + } +} + +Real +HomogenizedTotalLagrangianStressDivergenceA::computeQpOffDiagJacobianScalar(unsigned int svar_num) +{ + if (svar_num == _kappa_var) + return _dpk1[_qp].contractionKl(_m, _n, gradTest(_alpha)); + else + return 0.; +} + +Real +HomogenizedTotalLagrangianStressDivergenceA::computeScalarQpOffDiagJacobian(unsigned int jvar_num) +{ + if (jvar_num == _var.number()) + { + if (_ctype == HomogenizationA::ConstraintType::Stress) + return _dpk1[_qp].contractionIj(_m, _n, gradTrial(_alpha)); + else if (_ctype == HomogenizationA::ConstraintType::Strain) + if (_large_kinematics) + return Real(_m == _alpha) * gradTrial(_alpha)(_m, _n); + else + return 0.5 * (Real(_m == _alpha) * gradTrial(_alpha)(_m, _n) + + Real(_n == _alpha) * gradTrial(_alpha)(_n, _m)); + else + mooseError("Unknown constraint type in kernel calculation!"); + } + else + return 0.; +} + +void +HomogenizedTotalLagrangianStressDivergenceA::computeScalarOffDiagJacobianScalar( + const unsigned int svar_num) +{ + // Only do this for the other macro variable + if (svar_num == _kappao_var) + { + _local_ke.resize(_k_order, _ko_order); + + // Assemble J-kappa-kappa_other if _alpha==0 + if (_alpha == 0) + { + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + initScalarQpJacobian(_kappa_var); + Real dV = _JxW[_qp] * _coord[_qp]; + + _h = 0; + for (auto && [indices1, constraint1] : _cmap) + { + auto && [i, j] = indices1; + auto && ctype = constraint1.first; + + // identical logic to computeScalarResidual + if (_beta == 0) + { + if (_h == 1) + break; + } + else + { + if (_h == 0) + { + _h++; + continue; + } + } + + _l = 0; + for (auto && indices2 : _cmap) + { + auto && [a, b] = indices2.first; + + // OPPOSITE logic/scalar from computeScalarResidual, AND for column index + if (_beta == 1) + { + if (_l == 1) // only assemble first=0 component of _hvar, then break the loop + break; + } + else + { + if (_l == 0) // skip first component (_hvar) & continue to "first" component of _avar + { + _l++; + continue; + } + } + + unsigned int c_ind = + -(1 - _beta) + _l; // DON'T move 1 column left if _beta=1 for the other scalar + _l++; + if (ctype == HomogenizationA::ConstraintType::Stress) + _local_ke(-_beta + _h, c_ind) += dV * (_dpk1[_qp](i, j, a, b)); + else if (ctype == HomogenizationA::ConstraintType::Strain) + { + if (_large_kinematics) + _local_ke(-_beta + _h, c_ind) += dV * (Real(i == a && j == b)); + else + _local_ke(-_beta + _h, c_ind) += + dV * (0.5 * Real(i == a && j == b) + 0.5 * Real(i == b && j == a)); + } + else + mooseError("Unknown constraint type in Jacobian calculator!"); + } + _h++; + } + } + } + + for (const auto & matrix_tag : _matrix_tags) + _assembly.cacheJacobianBlock(_local_ke, + _kappa_var_ptr->dofIndices(), + _kappao_var_ptr->dofIndices(), + _kappa_var_ptr->scalingFactor(), + matrix_tag); + } +} diff --git a/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceS.C b/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceS.C index 153331bc4f46..de07d221946a 100644 --- a/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceS.C +++ b/modules/tensor_mechanics/test/src/kernels/HomogenizedTotalLagrangianStressDivergenceS.C @@ -12,12 +12,6 @@ // MOOSE includes #include "Function.h" #include "MooseVariableScalar.h" -// #include "Assembly.h" -// #include "MooseVariableFE.h" -// #include "MooseVariableScalar.h" -// #include "SystemBase.h" - -// #include "libmesh/quadrature.h" registerMooseObject("TensorMechanicsTestApp", HomogenizedTotalLagrangianStressDivergenceS); @@ -84,13 +78,14 @@ HomogenizedTotalLagrangianStressDivergenceS::computeScalarResidual() { std::vector scalar_residuals(_k_order); - if (_alpha == 0) // only do for the first component + // only assemble scalar residual once; i.e. when handling the first displacement component + if (_alpha == 0) { for (_qp = 0; _qp < _qrule->n_points(); _qp++) { initScalarQpResidual(); Real dV = _JxW[_qp] * _coord[_qp]; - _h = 0; + _h = 0; // single index for residual vector; double indices for constraint tensor component for (auto && [indices, constraint] : _cmap) { auto && [i, j] = indices; @@ -131,7 +126,8 @@ HomogenizedTotalLagrangianStressDivergenceS::computeScalarJacobian() { _local_ke.resize(_k_order, _k_order); - if (_alpha == 0) // only do for the first component + // only assemble scalar residual once; i.e. when handling the first displacement component + if (_alpha == 0) { for (_qp = 0; _qp < _qrule->n_points(); _qp++) { @@ -178,34 +174,27 @@ HomogenizedTotalLagrangianStressDivergenceS::computeScalarOffDiagJacobian( const unsigned int jvar_num) { const auto & jvar = getVariable(jvar_num); - if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_STANDARD) - { - // Get dofs and order of this variable; at least one will be _var - // const auto & jv0 = static_cast(jvar); - // const auto & loc_phi = jv0.phi(); - const auto jvar_size = jvar.phiSize(); - _local_ke.resize(_k_order, jvar_size); + // Get dofs and order of this variable; at least one will be _var + const auto jvar_size = jvar.phiSize(); + _local_ke.resize(_k_order, jvar_size); - for (_qp = 0; _qp < _qrule->n_points(); _qp++) + for (_qp = 0; _qp < _qrule->n_points(); _qp++) + { + // single index for Jacobian column; double indices for constraint tensor component + unsigned int h = 0; + Real dV = _JxW[_qp] * _coord[_qp]; + for (auto && [indices, constraint] : _cmap) { - unsigned int h = 0; - Real dV = _JxW[_qp] * _coord[_qp]; - for (auto && [indices, constraint] : _cmap) - { - _m = indices.first; - _n = indices.second; - _ctype = constraint.first; - initScalarQpOffDiagJacobian(jvar); - for (_j = 0; _j < jvar_size; _j++) - _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num); - h++; - } + // copy constraint indices to protected variables to pass to Qp routine + _m = indices.first; + _n = indices.second; + _ctype = constraint.first; + initScalarQpOffDiagJacobian(jvar); + for (_j = 0; _j < jvar_size; _j++) + _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num); + h++; } } - else if (jvar.fieldType() == Moose::VarFieldType::VAR_FIELD_ARRAY) - mooseError("Array variable cannot be coupled into Kernel Scalar currently"); - else - mooseError("Vector variable cannot be coupled into Kernel Scalar currently"); for (const auto & matrix_tag : _matrix_tags) _assembly.cacheJacobianBlock(_local_ke, @@ -231,6 +220,7 @@ HomogenizedTotalLagrangianStressDivergenceS::computeOffDiagJacobianScalarLocal( Real dV = _JxW[_qp] * _coord[_qp]; for (auto && [indices, constraint] : _cmap) { + // copy constraint indices to protected variables to pass to Qp routine _m = indices.first; _n = indices.second; _ctype = constraint.first; diff --git a/modules/tensor_mechanics/test/src/materials/ComputeHomogenizedLagrangianStrainA.C b/modules/tensor_mechanics/test/src/materials/ComputeHomogenizedLagrangianStrainA.C new file mode 100644 index 000000000000..d2e05da995b1 --- /dev/null +++ b/modules/tensor_mechanics/test/src/materials/ComputeHomogenizedLagrangianStrainA.C @@ -0,0 +1,87 @@ +//* 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 "ComputeHomogenizedLagrangianStrainA.h" + +registerMooseObject("TensorMechanicsTestApp", ComputeHomogenizedLagrangianStrainA); + +InputParameters +ComputeHomogenizedLagrangianStrainA::validParams() +{ + InputParameters params = Material::validParams(); + params.addParam("base_name", "Material property base name"); + params.addParam("homogenization_gradient_name", + "homogenization_gradient", + "Name of the constant gradient field"); + params.addRequiredParam( + "constraint_types", + HomogenizationB::constraintType, + "Type of each constraint: strain, stress, or none. The types are specified in the " + "column-major order, and there must be 9 entries in total."); + params.addRequiredParam>( + "targets", "Functions giving the targets to hit for constraint types that are not none."); + params.addRequiredCoupledVar("macro_gradientA", + "Scalar field defining the 1st component of" + "macro gradient"); + params.addRequiredCoupledVar("macro_gradient", + "Scalar field defining the other components of" + "macro gradient"); + return params; +} + +ComputeHomogenizedLagrangianStrainA::ComputeHomogenizedLagrangianStrainA( + const InputParameters & parameters) + : Material(parameters), + _base_name(isParamValid("base_name") ? getParam("base_name") + "_" : ""), + _macro_gradient(coupledScalarValue("macro_gradient")), + _macro_gradientA(coupledScalarValue("macro_gradientA")), + _homogenization_contribution( + declareProperty(_base_name + "homogenization_gradient_name")) +{ + // Constraint types + auto types = getParam("constraint_types"); + if (types.size() != Moose::dim * Moose::dim) + mooseError("Number of constraint types must equal dim * dim. ", types.size(), " are provided."); + + // Targets to hit + const std::vector & fnames = getParam>("targets"); + + // Prepare the constraint map + unsigned int fcount = 0; + for (const auto j : make_range(Moose::dim)) + for (const auto i : make_range(Moose::dim)) + { + const auto idx = i + Moose::dim * j; + const auto ctype = static_cast(types.get(idx)); + if (ctype != HomogenizationB::ConstraintType::None) + { + const Function * const f = &getFunctionByName(fnames[fcount++]); + _cmap[{i, j}] = {ctype, f}; + } + } +} + +void +ComputeHomogenizedLagrangianStrainA::computeQpProperties() +{ + _homogenization_contribution[_qp].zero(); + unsigned int count = 0; + for (auto && indices : _cmap) + { + auto && [i, j] = indices.first; + if (count == 0) + _homogenization_contribution[_qp](i, j) = _macro_gradientA[count++]; + else + { + unsigned int r_ind = count - 1; + count++; + _homogenization_contribution[_qp](i, j) = _macro_gradient[r_ind]; + } + } +} diff --git a/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/small-tests/2dscalar.i b/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/2dsole.i similarity index 82% rename from modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/small-tests/2dscalar.i rename to modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/2dsole.i index a42a749f9ee1..54d1619429ac 100644 --- a/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/small-tests/2dscalar.i +++ b/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/2dsole.i @@ -3,7 +3,6 @@ [GlobalParams] displacements = 'disp_x disp_y' large_kinematics = false - macro_gradient = hvar [] [Mesh] @@ -31,7 +30,11 @@ [] [hvar] family = SCALAR - order = THIRD + order = FIRST + [] + [hvarA] + family = SCALAR + order = SECOND [] [] @@ -108,21 +111,47 @@ [] [Kernels] - [sdx] - type = HomogenizedTotalLagrangianStressDivergenceS + [sdx0] + type = HomogenizedTotalLagrangianStressDivergenceA variable = disp_x component = 0 coupled_scalar = hvar macro_var = hvar + macro_other = hvarA + prime_scalar = 0 constraint_types = ${constraint_types} targets = ${targets} [] - [sdy] - type = HomogenizedTotalLagrangianStressDivergenceS + [sdy0] + type = HomogenizedTotalLagrangianStressDivergenceA variable = disp_y component = 1 coupled_scalar = hvar macro_var = hvar + macro_other = hvarA + prime_scalar = 0 + constraint_types = ${constraint_types} + targets = ${targets} + [] + [sdx1] + type = HomogenizedTotalLagrangianStressDivergenceA + variable = disp_x + component = 0 + coupled_scalar = hvarA + macro_var = hvarA + macro_other = hvar + prime_scalar = 1 + constraint_types = ${constraint_types} + targets = ${targets} + [] + [sdy1] + type = HomogenizedTotalLagrangianStressDivergenceA + variable = disp_y + component = 1 + coupled_scalar = hvarA + macro_var = hvarA + macro_other = hvar + prime_scalar = 1 constraint_types = ${constraint_types} targets = ${targets} [] @@ -137,6 +166,10 @@ type = NullScalarKernel variable = hvar [] + [nullA] + type = NullScalarKernel + variable = hvarA + [] [] [Functions] @@ -226,7 +259,9 @@ homogenization_gradient_names = 'homogenization_gradient' [] [compute_homogenization_gradient] - type = ComputeHomogenizedLagrangianStrainS + type = ComputeHomogenizedLagrangianStrainA + macro_gradientA = hvar + macro_gradient = hvarA constraint_types = ${constraint_types} targets = ${targets} [] diff --git a/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/gold/2dsole_out.csv b/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/gold/2dsole_out.csv new file mode 100644 index 000000000000..78352e5d7f98 --- /dev/null +++ b/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/gold/2dsole_out.csv @@ -0,0 +1,7 @@ +time,exx,exy,eyy,sxx,sxy,syy,hvar,hvarA_0,hvarA_1 +0,0,0,0,0,0,0,0,0,0 +0.2,0.00096026311741124,0.00030401429687746,-0.00080769950780423,80.000000000007,20,-39.999999999991,0.00096026311741122,0.00060802859375492,-0.00080769950780428 +0.4,0.0019205262348225,0.00060802859375492,-0.0016153990156086,160.00000000001,40,-80,0.0019205262348224,0.0012160571875098,-0.0016153990156086 +0.6,0.0028807893522336,0.00091204289063238,-0.0024230985234128,239.99999999999,60,-120.00000000001,0.0028807893522336,0.0018240857812648,-0.0024230985234128 +0.8,0.0038410524696449,0.0012160571875098,-0.0032307980312172,319.99999999999,80,-160.00000000001,0.0038410524696449,0.0024321143750197,-0.0032307980312171 +1,0.0048013155870561,0.0015200714843873,-0.0040384975390213,400,100,-199.99999999999,0.0048013155870561,0.0030401429687746,-0.0040384975390213 diff --git a/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/tests b/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/tests index f75fecdfff70..2e3e2d78f359 100644 --- a/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/tests +++ b/modules/tensor_mechanics/test/tests/lagrangian/cartesian/total/homogenization/scalar_kernel/tests @@ -2,7 +2,7 @@ issues = '#22174' design = 'source/kernels/ScalarLMKernel.md' - [2d-strain] + [2d-strainS] type = CSVDiff input = '2dscalar.i' csvdiff = '2d-strain_out.csv' @@ -11,7 +11,7 @@ requirement = "Homogenization from scalar wrapper class with strain constraints hits the targets in 2D" allow_test_objects = true [] - [2d-stress] + [2d-stressS] type = CSVDiff input = '2dscalar.i' csvdiff = '2d-stress_out.csv' @@ -20,7 +20,7 @@ requirement = "Homogenization from scalar wrapper class with stress constraints hits the targets in 2D" allow_test_objects = true [] - [2d-mixed] + [2d-mixedS] type = CSVDiff input = '2dscalar.i' csvdiff = '2d-mixed_out.csv' @@ -29,4 +29,13 @@ requirement = "Homogenization from scalar wrapper class with mixed stress and strain constraints hits the targets in 2D" allow_test_objects = true [] + [2d-stressA] + type = CSVDiff + input = '2dsole.i' + csvdiff = '2dsole_out.csv' + cli_args = "constraint_types='stress none none stress stress " + "none none none none' targets='stress11 stress12 stress22'" + requirement = "Framework scalar kernel wrapper correctly assembles scalar-to-scalar coupling Jacobian" + allow_test_objects = true + [] []