Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "DensityUpdate.h"
      11             : #include <algorithm>
      12             : 
      13             : registerMooseObject("OptimizationApp", DensityUpdate);
      14             : 
      15             : InputParameters
      16          82 : DensityUpdate::validParams()
      17             : {
      18          82 :   InputParameters params = ElementUserObject::validParams();
      19          82 :   params.addClassDescription(
      20             :       "Compute updated densities based on sensitivities using an optimality criteria method to "
      21             :       "keep the volume constraint satisified.");
      22         164 :   params.addRequiredCoupledVar("design_density", "Design density variable name.");
      23         164 :   params.addRequiredParam<VariableName>("density_sensitivity",
      24             :                                         "Name of the density_sensitivity variable.");
      25         164 :   params.addRequiredParam<Real>("volume_fraction", "Volume Fraction");
      26         164 :   params.addParam<Real>("bisection_lower_bound", 0, "Lower bound for the bisection algorithm.");
      27         164 :   params.addParam<Real>("bisection_upper_bound", 1e16, "Upper bound for the bisection algorithm.");
      28          82 :   return params;
      29           0 : }
      30             : 
      31          44 : DensityUpdate::DensityUpdate(const InputParameters & parameters)
      32             :   : ElementUserObject(parameters),
      33          44 :     _mesh(_subproblem.mesh()),
      34          44 :     _density_sensitivity_name(getParam<VariableName>("density_sensitivity")),
      35          44 :     _design_density(&writableVariable("design_density")),
      36          44 :     _density_sensitivity(&_subproblem.getStandardVariable(_tid, _density_sensitivity_name)),
      37          88 :     _volume_fraction(getParam<Real>("volume_fraction")),
      38          88 :     _lower_bound(getParam<Real>("bisection_lower_bound")),
      39         132 :     _upper_bound(getParam<Real>("bisection_upper_bound"))
      40             : {
      41          44 :   if (!dynamic_cast<MooseVariableFE<Real> *>(_design_density))
      42           0 :     paramError("design_density", "Design density must be a finite element variable");
      43          44 : }
      44             : 
      45             : void
      46         512 : DensityUpdate::timestepSetup()
      47             : {
      48         512 :   gatherElementData();
      49         512 :   performOptimCritLoop();
      50         512 : }
      51             : 
      52             : void
      53       83776 : DensityUpdate::execute()
      54             : {
      55             :   // Grab the element data for each id
      56       83776 :   auto elem_data_iter = _elem_data_map.find(_current_elem->id());
      57             : 
      58             :   // Check if the element data is not null
      59       83776 :   if (elem_data_iter != _elem_data_map.end())
      60             :   {
      61             :     ElementData & elem_data = elem_data_iter->second;
      62       83776 :     dynamic_cast<MooseVariableFE<Real> *>(_design_density)->setNodalValue(elem_data.new_density);
      63             :   }
      64             :   else
      65             :   {
      66           0 :     mooseError("Element data not found for the current element id.");
      67             :   }
      68       83776 : }
      69             : 
      70             : void
      71         512 : DensityUpdate::gatherElementData()
      72             : {
      73             :   _elem_data_map.clear();
      74         512 :   _total_allowable_volume = 0;
      75             : 
      76        1024 :   for (const auto & sub_id : blockIDs())
      77      103280 :     for (const auto & elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(sub_id))
      78             :     {
      79      102256 :       dof_id_type elem_id = elem->id();
      80             :       ElementData data = ElementData(
      81      102256 :           dynamic_cast<MooseVariableFE<Real> *>(_design_density)->getElementalValue(elem),
      82      102256 :           dynamic_cast<const MooseVariableFE<Real> *>(_density_sensitivity)
      83             :               ->getElementalValue(elem),
      84      102256 :           elem->volume(),
      85      306768 :           0);
      86      102256 :       _elem_data_map[elem_id] = data;
      87      102256 :       _total_allowable_volume += elem->volume();
      88         512 :     }
      89             : 
      90         512 :   _communicator.sum(_total_allowable_volume);
      91         512 :   _total_allowable_volume *= _volume_fraction;
      92         512 : }
      93             : 
      94             : void
      95         512 : DensityUpdate::performOptimCritLoop()
      96             : {
      97             :   // Initialize the lower and upper bounds for the bisection method
      98         512 :   Real l1 = _lower_bound;
      99         512 :   Real l2 = _upper_bound;
     100             :   bool perform_loop = true;
     101             :   // Loop until the relative difference between l1 and l2 is less than a small tolerance
     102       31914 :   while (perform_loop)
     103             :   {
     104             :     // Compute the midpoint between l1 and l2
     105       31402 :     Real lmid = 0.5 * (l2 + l1);
     106             : 
     107             :     // Initialize the current total volume
     108       31402 :     Real curr_total_volume = 0;
     109             :     // Loop over all elements
     110     6302898 :     for (auto && [id, elem_data] : _elem_data_map)
     111             :     {
     112             :       // Compute the updated density for the current element
     113     6271496 :       Real new_density = computeUpdatedDensity(elem_data.old_density, elem_data.sensitivity, lmid);
     114             :       // Update the current filtered density for the current element
     115     6271496 :       elem_data.new_density = new_density;
     116             :       // Update the current total volume
     117     6271496 :       curr_total_volume += new_density * elem_data.volume;
     118             :     }
     119             : 
     120             :     // Sum the current total volume across all processors
     121       31402 :     _communicator.sum(curr_total_volume);
     122             : 
     123             :     // Update l1 or l2 based on whether the current total volume is greater than the total
     124             :     // allowable volume
     125       31402 :     if (curr_total_volume > _total_allowable_volume)
     126             :       l1 = lmid;
     127             :     else
     128             :       l2 = lmid;
     129             : 
     130             :     // Determine whether to continue the loop based on the relative difference between l1 and l2
     131       31402 :     perform_loop = (l2 - l1) / (l1 + l2) > 1e-3;
     132             :   }
     133         512 : }
     134             : 
     135             : // Method to compute the updated density for an element
     136             : Real
     137     6271496 : DensityUpdate::computeUpdatedDensity(Real current_density, Real dc, Real lmid)
     138             : {
     139             :   // Define the maximum allowable change in density
     140             :   Real move = 0.5;
     141             :   // Compute the updated density based on the current density, the sensitivity, and the midpoint
     142             :   // value
     143             :   Real updated_density = std::max(
     144     6271496 :       0.0,
     145     6271496 :       std::min(1.0,
     146     6271496 :                std::min(current_density + move,
     147    12542992 :                         std::max(current_density - move,
     148    12510784 :                                  std::max(1e-5, current_density) * std::sqrt(-dc / lmid)))));
     149             :   // Return the updated density
     150     6271496 :   return updated_density;
     151             : }