Merge pull request #78 from dewenyushu/elastic_inversion

Add elastic inversion capability using batch materials

doc/content/source/userobjects/BatchStressGrad.md

@@ -0,0 +1,22 @@
+# BatchStressGrad
+
+## Description
+
+This `UserObject` computes double contraction of the elasticity tensor derivative and the forward mechanical strain, i.e.,
+\begin{equation}
+\underbrace{\frac{\partial \boldsymbol{C}}{\partial p}}_{\text{elasticity tensor derivative}} \underbrace{(\boldsymbol{L} u)}_{\text{forward strain}}
+\end{equation}
+as a batch material.
+Here, $\boldsymbol{C}$ is the elasticity tensor, $p$ is the interested parameter, $\boldsymbol{L}$ is the symmetric strain operator for the elasticity problem, and $\boldsymbol{L} u$ is the strain from the forward problem. This object requires the elasticity tensor derivative material property (i.e., $\frac{\partial \boldsymbol{C}}{\partial p}$) as its input.
+
+This object is used together with [AdjointStrainStressGradInnerProduct](/AdjointStrainStressGradInnerProduct.md) in a linear elastic inversion problem.
+
+## Example Input File Syntax
+
+!listing test/tests/bimaterial_elastic_inversion_batch_mat/grad.i block=UserObjects/stress_grad_lambda
+
+!syntax parameters /UserObjects/BatchStressGrad
+
+!syntax inputs /UserObjects/BatchStressGrad
+
+!syntax children /UserObjects/BatchStressGrad

doc/content/source/vectorpostprocessors/AdjointStrainStressGradInnerProduct.md

@@ -0,0 +1,31 @@
+# AdjointStrainStressGradInnerProduct
+
+!syntax description /VectorPostprocessors/AdjointStrainStressGradInnerProduct
+
+## Description
+
+This `VectorPostprocessor` is utilized in a mechanical inverse optimization problem (see [Inverse Optimization Theory](/InvOptTheory.md)) for computing the gradient of objective function with respect to a parameter of interest.
+
+Specifically, the objective function is defined as the sum of the regularization and the displacement mismatch (between the simulation and the measurement data at specific locations) as follows:
+\begin{equation}
+  F({u}, \boldsymbol{p}) = ||\boldsymbol{S}{u} - \boldsymbol{u}_m ||^2 + R({u}, \boldsymbol{p}),
+\end{equation}
+where ${u}$ denotes the displacement from the forward problem, $\boldsymbol{S}$ is the sampling operator, $\boldsymbol{u}_m$ are the measurements, $\boldsymbol{p} = [p_1, p_2, ..., p_n]$ is the vector of parameters,  and $R$ is the regularization term.
+
+The gradient of objective function with respect to one scalar parameter ($p_i$) is
+\begin{equation}
+  \frac{\text{d}{F}}{\text{d}p_i} = - \int \underbrace{(\boldsymbol{L}\lambda)^T}_{\text{adjoint strain}}  \underbrace{\frac{\partial \boldsymbol{C}}{\partial p_i} (\boldsymbol{L} u)}_{\text{gradient of stress}} \text{d}\Omega,
+\end{equation}
+where ${F}$ is the objective function, $\boldsymbol{L}$ is the symmetric strain operator for the elasticity problem, $\lambda$ is the adjoint displacement, $\boldsymbol{C}$ is the elasticity tensor, and $p_i$ is the interested scalar parameter.
+
+Specifically, this object takes the `adjoint strain` and the `gradient of stress w.r.t. parameter` as inputs, computes the double contraction of the two inputs, and integrate over the entire domain. Specifically, the `adjoint strain` is the strain from the adjoint problem, the `gradient of stress w.r.t. parameter` should be a batch material (e.g., [BatchStressGrad](/BatchStressGrad.md)).
+
+## Example Input File Syntax
+
+!listing test/tests/bimaterial_elastic_inversion_batch_mat/grad.i block=VectorPostprocessors/grad_lambda
+
+!syntax parameters /VectorPostprocessors/AdjointStrainStressGradInnerProduct
+
+!syntax inputs /VectorPostprocessors/AdjointStrainStressGradInnerProduct
+
+!syntax children /VectorPostprocessors/AdjointStrainStressGradInnerProduct

include/userobjects/BatchStressGrad.h

@@ -0,0 +1,33 @@
+//* 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 "BatchMaterial.h"
+
+typedef BatchMaterial<
+    // tuple representation
+    BatchMaterialUtils::TupleStd,
+    // output data type
+    RankTwoTensor,
+    // gathered input data types:
+    BatchMaterialUtils::GatherMatProp<RankFourTensor>,
+    BatchMaterialUtils::GatherMatProp<RankTwoTensor>>
+
+    BatchStressGradParent;
+
+class BatchStressGrad : public BatchStressGradParent
+{
+public:
+  static InputParameters validParams();
+
+  BatchStressGrad(const InputParameters & params);
+
+  void batchCompute() override;
+};

include/vectorpostprocessors/AdjointStrainStressGradInnerProduct.h

@@ -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
+
+#pragma once
+
+#include "ElementOptimizationFunctionInnerProduct.h"
+#include "BatchStressGrad.h"
+
+class AdjointStrainStressGradInnerProduct : public ElementOptimizationFunctionInnerProduct
+{
+public:
+  static InputParameters validParams();
+
+  AdjointStrainStressGradInnerProduct(const InputParameters & parameters);
+
+protected:
+  virtual Real computeQpInnerProduct() override;
+  /// Base name of the material system
+  const std::string _base_name;
+  /// Holds adjoint strain at current quadrature points
+  const MaterialProperty<RankTwoTensor> & _adjoint_strain;
+  /// UO that holds gradient of stress wrt material parameter at current quadrature points
+  const BatchStressGrad & _stress_grad_uo;
+  const BatchStressGrad::OutputVector & _stress_grad_uo_output;
+};

src/userobjects/BatchStressGrad.C

@@ -0,0 +1,58 @@
+//* 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 "BatchStressGrad.h"
+#include "libmesh/int_range.h"
+
+registerMooseObject("isopodApp", BatchStressGrad);
+
+InputParameters
+BatchStressGrad::validParams()
+{
+  auto params = BatchStressGradParent::validParams();
+  params.addRequiredParam<MaterialPropertyName>(
+      "elasticity_tensor_derivative",
+      "Name of the elasticity tensor derivative material property.");
+  params.setDocString(
+      "execution_order_group",
+      "BatchStressGrad userObject needs to be completely executed before vectorPostprocessors.");
+  params.set<int>("execution_order_group") = -1;
+  params.addClassDescription("Compute the double contraction of the elasticity tensor derivative "
+                             "and the forward mechanical strain");
+  return params;
+}
+
+BatchStressGrad::BatchStressGrad(const InputParameters & params)
+  : BatchStressGradParent(params,
+                          // here we pass the derivative of elasticity tensor wrt to the parameter
+                          "elasticity_tensor_derivative",
+                          // here we pass in the forward strain
+                          "forward_mechanical_strain")
+{
+}
+
+/*
+  Note: The following calculation is currently specialized for a linear elastic inverse optimization
+  problem. This will be swapped out later on when we can get the gradient of stress w.r.t. the
+  interested parameter from NEML for general material models.
+*/
+void
+BatchStressGrad::batchCompute()
+{
+  for (const auto i : index_range(_input_data))
+  {
+    const auto & input = _input_data[i];
+    auto & output = _output_data[i];
+
+    const auto & elasticity_dev = std::get<0>(input);
+    const auto & strain = std::get<1>(input);
+
+    output = elasticity_dev * strain;
+  }
+}

src/vectorpostprocessors/AdjointStrainStressGradInnerProduct.C

@@ -0,0 +1,49 @@
+//* 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 "AdjointStrainStressGradInnerProduct.h"
+
+registerMooseObject("isopodApp", AdjointStrainStressGradInnerProduct);
+
+InputParameters
+AdjointStrainStressGradInnerProduct::validParams()
+{
+  InputParameters params = ElementOptimizationFunctionInnerProduct::validParams();
+  params.addRequiredParam<UserObjectName>(
+      "stress_grad_name",
+      "Name of the stress gradient user object with respect to the material parameter");
+
+  params.addRequiredParam<MaterialPropertyName>(
+      "adjoint_strain_name", "Name of the strain property in the adjoint problem");
+
+  params.addClassDescription(
+      "Compute the parameter gradient for linear elastic material inversion");
+  return params;
+}
+AdjointStrainStressGradInnerProduct::AdjointStrainStressGradInnerProduct(
+    const InputParameters & parameters)
+  : ElementOptimizationFunctionInnerProduct(parameters),
+    _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
+    _adjoint_strain(getMaterialPropertyByName<RankTwoTensor>(
+        getParam<MaterialPropertyName>("adjoint_strain_name"))),
+    _stress_grad_uo(getUserObject<BatchStressGrad>("stress_grad_name")),
+    _stress_grad_uo_output(_stress_grad_uo.getOutputData())
+{
+}
+
+Real
+AdjointStrainStressGradInnerProduct::computeQpInnerProduct()
+{
+  if (!_stress_grad_uo.outputReady())
+    mooseError("Stress gradient batch material property is not ready to output.");
+
+  const auto index = _stress_grad_uo.getIndex(_current_elem->id());
+
+  return -_adjoint_strain[_qp].doubleContraction(_stress_grad_uo_output[index + _qp]);
+}

test/tests/bimaterial_elastic_inversion_batch_mat/gold/main_out_OptimizationReporter_0001.csv

@@ -0,0 +1,10 @@
+grad_lambda,grad_mu,lambda,measurement_time,measurement_values,measurement_xcoord,measurement_ycoord,measurement_zcoord,misfit_values,mu,simulation_values
+5.1698615078396e-11,5.2394323056988e-11,4.9999752463318,0,4.46462,-1,-1,0,9.2479240265675e-08,1.0000001055133,4.4646200924792
+5.097187284528e-12,3.5417284419815e-11,4.0000666870441,0,6.447155,-1,0,0,5.2641227377137e-08,2.0000064790597,6.4471550526412
+-2.7212788369455e-12,8.6258766757591e-12,2.9997143457064,0,8.434803,-1,1,0,-1.0044113984975e-07,2.9999940961679,8.4348028995589
+0,0,0,0,4.176264,0,-1,0,-1.14358395642e-07,0,4.1762638856416
+0,0,0,0,6.172984,0,0,0,-1.7726937251439e-07,0,6.1729838227306
+0,0,0,0,8.200859,0,1,0,1.9391855410333e-07,0,8.2008591939186
+0,0,0,0,4.049477,1,-1,0,2.2853767767117e-08,0,4.0494770228538
+0,0,0,0,6.052499,1,0,0,1.1289934853664e-07,0,6.0524991128993
+0,0,0,0,8.079385,1,1,0,-8.6613697547477e-08,0,8.0793849133863