compliance sensitivity computed with a derivative parsed material closes

idaholab#3
recuero · Oct 13, 2023 · ce561e5 · ce561e5
1 parent ae3f67e
commit ce561e5
Show file tree

Hide file tree

Showing 6 changed files with 319 additions and 30 deletions.
diff --git a/include/materials/ComplianceSensitivity2.h b/include/materials/ComplianceSensitivity2.h
@@ -0,0 +1,29 @@
+//* 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 "StrainEnergyDensity.h"
+#include "MathUtils.h"
+
+class ComplianceSensitivity2 : public StrainEnergyDensity
+{
+public:
+  static InputParameters validParams();
+
+  ComplianceSensitivity2(const InputParameters & parameters);
+
+  virtual void computeQpProperties() override;
+
+protected:
+  MaterialProperty<Real> & _sensitivity;
+  const VariableValue & _design_density;
+  const MaterialPropertyName _design_density_name;
+  const MaterialProperty<Real> & _dEdp;
+};
diff --git a/src/materials/ComplianceSensitivity.C b/src/materials/ComplianceSensitivity.C
@@ -45,7 +45,6 @@ ComplianceSensitivity::computeQpProperties()
   Real compliance = _strain_energy_density[_qp];
 
   // Compute the derivative of the compliance with respect to the design density
-  // _power-2 because StrainEnergyDensity needed to be divided by the _design_density
   Real derivative =
       -_power * (_E - _Emin) * MathUtils::pow(_design_density[_qp], _power - 1) * compliance;
 

diff --git a/src/materials/ComplianceSensitivity2.C b/src/materials/ComplianceSensitivity2.C
@@ -0,0 +1,42 @@
+//* 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 "ComplianceSensitivity2.h"
+
+registerMooseObject("troutApp", ComplianceSensitivity2);
+
+InputParameters
+ComplianceSensitivity2::validParams()
+{
+  InputParameters params = StrainEnergyDensity::validParams();
+  params.addClassDescription("Computes compliance sensitivity needed for SIMP method.");
+  params.addRequiredCoupledVar("design_density", "Design density variable name.");
+  params.addRequiredParam<MaterialPropertyName>("youngs_modulus",
+                                                "DerivativeParsedMaterial for Youngs modulus.");
+
+  return params;
+}
+
+ComplianceSensitivity2::ComplianceSensitivity2(const InputParameters & parameters)
+  : StrainEnergyDensity(parameters),
+    _sensitivity(declareProperty<Real>(_base_name + "sensitivity")),
+    _design_density(coupledValue("design_density")),
+    _design_density_name(coupledName("design_density", 0)),
+    _dEdp(getMaterialPropertyDerivativeByName<Real>(
+        getParam<MaterialPropertyName>("youngs_modulus"), _design_density_name))
+{
+}
+
+void
+ComplianceSensitivity2::computeQpProperties()
+{
+  // Call the parent class's method to compute the strain energy density
+  StrainEnergyDensity::computeQpProperties();
+  _sensitivity[_qp] = -_dEdp[_qp] * _strain_energy_density[_qp];
+}
diff --git a/test/tests/materials/ComplianceSensitivity/2d_mbb_ComplianceSensitivity2.i b/test/tests/materials/ComplianceSensitivity/2d_mbb_ComplianceSensitivity2.i
@@ -0,0 +1,178 @@
+vol_frac = 0.5
+E0 = 1
+Emin = 1e-8
+power = 2
+[GlobalParams]
+  displacements = 'disp_x disp_y'
+[]
+
+[Mesh]
+  [MeshGenerator]
+    type = GeneratedMeshGenerator
+    dim = 2
+    nx = 60
+    ny = 20
+    xmin = 0
+    xmax = 30
+    ymin = 0
+    ymax = 10
+  []
+  [node]
+    type = ExtraNodesetGenerator
+    input = MeshGenerator
+    new_boundary = hold
+    nodes = 0
+  []
+  [push]
+    type = ExtraNodesetGenerator
+    input = node
+    new_boundary = push
+    coord = '30 10 0'
+  []
+[]
+
+[Variables]
+  [disp_x]
+  []
+  [disp_y]
+  []
+[]
+
+[AuxVariables]
+  [Dc]
+    family = MONOMIAL
+    order = CONSTANT
+    initial_condition = -1.0
+  []
+  [mat_den]
+    family = MONOMIAL
+    order = CONSTANT
+    initial_condition = ${vol_frac}
+  []
+[]
+
+[Modules/TensorMechanics/Master]
+  [all]
+    strain = SMALL
+    add_variables = true
+    incremental = false
+  []
+[]
+
+[BCs]
+  [no_x]
+    type = DirichletBC
+    variable = disp_y
+    boundary = hold
+    value = 0.0
+  []
+  [no_y]
+    type = DirichletBC
+    variable = disp_x
+    boundary = right
+    value = 0.0
+  []
+
+[]
+[NodalKernels]
+  [push]
+    type = NodalGravity
+    variable = disp_y
+    boundary = push
+    gravity_value = -1
+    mass = 1
+  []
+[]
+[Materials]
+  [elasticity_tensor]
+    type = ComputeVariableIsotropicElasticityTensor
+    youngs_modulus = E_phys
+    poissons_ratio = poissons_ratio
+    args = 'mat_den'
+  []
+  [E_phys]
+    type = DerivativeParsedMaterial
+    # Emin + (density^penal) * (E0 - Emin)
+    expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
+    coupled_variables = 'mat_den'
+    property_name = E_phys
+  []
+
+  [poissons_ratio]
+    type = GenericConstantMaterial
+    prop_names = poissons_ratio
+    prop_values = 0.3
+  []
+  [stress]
+    type = ComputeLinearElasticStress
+  []
+  [dc]
+    type = ComplianceSensitivity2
+    design_density = mat_den
+    youngs_modulus = E_phys
+  []
+[]
+
+[Preconditioning]
+  [smp]
+    type = SMP
+    full = true
+  []
+[]
+[UserObjects]
+  [rad_avg]
+    type = RadialAverage
+    radius = 1.2
+    weights = linear
+    prop_name = sensitivity
+    execute_on = TIMESTEP_END
+    force_preaux = true
+  []
+  [update]
+    type = DensityUpdate
+    density_sensitivity = Dc
+    design_density = mat_den
+    power = ${power}
+    volume_fraction = ${vol_frac}
+    execute_on = TIMESTEP_BEGIN
+  []
+  [calc_sense]
+    type = SensitivityFilter
+    density_sensitivity = Dc
+    design_density = mat_den
+    filter_UO = rad_avg
+    execute_on = TIMESTEP_END
+    force_postaux = true
+  []
+[]
+
+[Executioner]
+  type = Transient
+  solve_type = NEWTON
+  petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
+  petsc_options_value = 'lu superlu_dist'
+  nl_abs_tol = 1e-8
+  dt = 1.0
+  num_steps = 30
+[]
+
+[Outputs]
+  exodus = true
+  [out]
+    type = CSV
+    execute_on = 'TIMESTEP_END'
+  []
+  print_linear_residuals = false
+[]
+
+[Postprocessors]
+  [total_vol]
+    type = ElementIntegralVariablePostprocessor
+    variable = mat_den
+    execute_on = 'INITIAL TIMESTEP_END'
+  []
+  [sensitivity]
+    type = ElementIntegralMaterialProperty
+    mat_prop = sensitivity
+  []
+[]
diff --git a/test/tests/materials/ComplianceSensitivity/gold/2d_mbb_ComplianceSensitivity2_out.csv b/test/tests/materials/ComplianceSensitivity/gold/2d_mbb_ComplianceSensitivity2_out.csv
@@ -0,0 +1,31 @@
+time,sensitivity,total_vol
+1,-229.1081493809,149.94613363944
+2,-161.4398863729,149.94819980152
+3,-162.27902819402,149.98391415911
+4,-167.39933759651,149.98572736093
+5,-176.06320742248,150.02154996884
+6,-189.26683075254,150.00751584886
+7,-195.7421073052,150.00094635493
+8,-193.52592659895,149.98817937627
+9,-189.7735786412,149.97932559696
+10,-185.55847044291,150.01812407771
+11,-182.66390735576,150.00831525335
+12,-181.02328473312,150.0153331652
+13,-179.98445432083,149.99243925594
+14,-179.02578014307,150.00060655735
+15,-178.29074239331,150.01089055168
+16,-177.78149904281,149.99588065028
+17,-177.36043213031,150.00457727799
+18,-177.04485802053,149.99754495777
+19,-176.79977652138,150.00350973286
+20,-176.60187423747,150.00089504899
+21,-176.43335893786,149.99990472782
+22,-176.28633568933,150.00638022164
+23,-176.15480897719,150.00450883322
+24,-176.02410972416,149.99781773574
+25,-175.91027708796,149.99474107163
+26,-175.81519603888,149.99840879314
+27,-175.72048156499,150.00028728316
+28,-175.6064063585,150.00539638885
+29,-175.52813806785,149.99831938762
+30,-175.46202104732,150.00443732769
diff --git a/test/tests/materials/ComplianceSensitivity/tests b/test/tests/materials/ComplianceSensitivity/tests
@@ -1,31 +1,41 @@
 [Tests]
-    issues = '#1'
-    [2d_mbb]
-      type = 'CSVDiff'
-      input = '2d_mbb.i'
-      csvdiff = '2d_mbb_out.csv'
-      requirement = 'The system shall be able to compute the compliance sensitivty correctly for a 2D topology optimization problem.'
-      heavy = true
-    []
-    [3d_mbb]
-      type = 'CSVDiff'
-      input = '3d_mbb.i'
-      csvdiff = '3d_mbb_out.csv'
-      requirement = 'The system shall be able to compute the compliance sensitivty correctly for a 3D topology optimization problem.'
-      heavy = true
-    []
-    [2d_mbb_pde]
-      type = 'CSVDiff'
-      input = '2d_mbb_pde.i'
-      csvdiff = '2d_mbb_pde_out.csv'
-      requirement = 'The system shall be able to compute the compliance sensitivty correctly for a 2D topology optimization problem with a PDE filter.'
-      heavy = true
-    []
-    [2d_mbb_pde_amr]
-      type = 'CSVDiff'
-      input = '2d_mbb_pde_amr.i'
-      csvdiff = '2d_mbb_pde_amr_out.csv'
-      requirement = 'The system shall be able to compute the compliance sensitivty correctly for a 2D topology optimization problem with AMR.'
-      heavy = true
-    []
+  issues = '#1'
+  requirement = 'The system shall be able to compute the compliance sensitivty correctly for a '
+  [2d_mbb]
+    type = 'CSVDiff'
+    input = '2d_mbb.i'
+    csvdiff = '2d_mbb_out.csv'
+    detail = '2D topology optimization problem.'
+    heavy = true
   []
+  [3d_mbb]
+    type = 'CSVDiff'
+    min_parallel = 4
+    input = '3d_mbb.i'
+    csvdiff = '3d_mbb_out.csv'
+    detail = '3D topology optimization problem.'
+    heavy = true
+  []
+  [2d_mbb_pde]
+    type = 'CSVDiff'
+    input = '2d_mbb_pde.i'
+    csvdiff = '2d_mbb_pde_out.csv'
+    detail = '2D topology optimization problem with a PDE filter.'
+    heavy = true
+  []
+  [2d_mbb_pde_amr]
+    type = 'CSVDiff'
+    input = '2d_mbb_pde_amr.i'
+    csvdiff = '2d_mbb_pde_amr_out.csv'
+    detail = '2D topology optimization problem with AMR.'
+    heavy = true
+  []
+  [2d_mbb_cs2]
+    issues = '#3'
+    type = 'CSVDiff'
+    input = '2d_mbb_ComplianceSensitivity2.i'
+    csvdiff = '2d_mbb_ComplianceSensitivity2_out.csv'
+    rel_err = 1e-1
+    detail = '2D topology optimization problem using a derivativeParsedMaterial.'
+  []
+[]