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             : #pragma once
      11             : 
      12             : #include "MooseTypes.h"
      13             : #include "KDTree.h"
      14             : 
      15             : #include "libmesh/id_types.h"
      16             : #include <unordered_map>
      17             : #include "libmesh/parallel.h"
      18             : #include "libmesh/fe_type.h"
      19             : #include "libmesh/point.h"
      20             : #include "libmesh/replicated_mesh.h"
      21             : #include "libmesh/equation_systems.h"
      22             : #include "libmesh/point_locator_base.h"
      23             : #include "libmesh/exodusII_io.h"
      24             : #include <vector>
      25             : #include <functional>
      26             : 
      27             : namespace libMesh
      28             : {
      29             : class System;
      30             : class DofMap;
      31             : }
      32             : 
      33             : using libMesh::RealGradient;
      34             : 
      35             : /**
      36             :  * Utility class to use an Exodus mesh to define controllable parameters for optimization problems
      37             :  * This class will:
      38             :  *  - Ensure that controllable parameters defined by the mesh are correctly ordered on optimiation
      39             :  * main app and forward or adjoint sub-apps
      40             :  *  - Define the parameter space (nodal/elemental and shape function)
      41             :  *  - Interpolate parameter values to the forward problem mesh using nodal/elemental shape function
      42             :  */
      43             : class ParameterMesh
      44             : {
      45             : public:
      46             :   ParameterMesh(const libMesh::FEType & param_type,
      47             :                 const std::string & exodus_mesh,
      48             :                 const bool find_closest = false,
      49             :                 const unsigned int kdtree_candidates = 5);
      50             : 
      51             :   /// Enumerations for regularization computations
      52             :   enum class RegularizationType
      53             :   {
      54             :     L2_GRADIENT,
      55             :     // Future regularization types can be added here:
      56             :     // L1,
      57             :     // H1,
      58             :     // TV (Total Variation)
      59             :   };
      60             : 
      61             :   /**
      62             :    * @return the number of parameters read from the mesh for a single timestep
      63             :    */
      64    24412015 :   dof_id_type size() const { return _param_dofs; }
      65             :   /**
      66             :    * Interpolate parameters onto the computational mesh
      67             :    * getIndexAndWeight is only used by ParameterMeshFunction
      68             :    * @param pt location to compute elemnent dof_indices weights
      69             :    * @param dof_indices return dof indices for element containing pt
      70             :    * @param weights returns element shape function weights at pt
      71             :    */
      72             :   void getIndexAndWeight(const Point & pt,
      73             :                          std::vector<dof_id_type> & dof_indices,
      74             :                          std::vector<Real> & weights) const;
      75             :   /**
      76             :    * Performs inner products of parameters with functions on the computational mesh
      77             :    * getIndexAndWeight is only used by ParameterMeshFunction
      78             :    * @param pt location to compute elemnent dof_indices weights
      79             :    * @param dof_indices return dof indices for element containing pt
      80             :    * @param weights returns element shape function gradient weights at pt
      81             :    */
      82             :   void getIndexAndWeight(const Point & pt,
      83             :                          std::vector<dof_id_type> & dof_indices,
      84             :                          std::vector<RealGradient> & weights) const;
      85             :   /**
      86             :    * Computes regularization objective value for a given regularization type
      87             :    * @param parameter_values  vector of parameter values to compute regularization for
      88             :    * @param reg_type  type of regularization (L2_GRADIENT, etc.)
      89             :    * @return scalar objective value
      90             :    */
      91             :   Real computeRegularizationObjective(const std::vector<Real> & parameter_values,
      92             :                                       RegularizationType reg_type) const;
      93             : 
      94             :   /**
      95             :    * Computes regularization gradient for a given regularization type
      96             :    * @param parameter_values  vector of parameter values to compute gradient for
      97             :    * @param reg_type  type of regularization (L2_GRADIENT, etc.)
      98             :    * @return vector of gradient values (same size as parameter_values)
      99             :    */
     100             :   std::vector<Real> computeRegularizationGradient(const std::vector<Real> & parameter_values,
     101             :                                                   RegularizationType reg_type) const;
     102             : 
     103             : private:
     104             :   /**
     105             :    * Template method containing the element loop for regularization computations
     106             :    * @param parameter_values  vector of parameter values
     107             :    * @param reg_type  type of regularization
     108             :    * @return result of type T (Real for objective, std::vector<Real> for gradient)
     109             :    */
     110             :   template <typename T>
     111             :   T computeRegularizationLoop(const std::vector<Real> & parameter_values,
     112             :                               RegularizationType reg_type) const;
     113             : 
     114             : protected:
     115             :   libMesh::Parallel::Communicator _communicator;
     116             :   libMesh::ReplicatedMesh _mesh;
     117             :   /// Find closest projection points
     118             :   const bool _find_closest;
     119             :   std::unique_ptr<libMesh::EquationSystems> _eq;
     120             :   libMesh::System * _sys;
     121             :   std::unique_ptr<libMesh::PointLocatorBase> _point_locator;
     122             :   std::unique_ptr<libMesh::ExodusII_IO> _exodusII_io;
     123             : 
     124             :   dof_id_type _param_dofs;
     125             : 
     126             :   /// Node-based KDTree optimization
     127             :   std::vector<Point> _mesh_nodes;
     128             :   std::unique_ptr<KDTree> _node_kdtree;
     129             :   std::unordered_map<dof_id_type, std::set<const libMesh::Elem *>> _node_to_elements;
     130             :   unsigned int _kdtree_candidates;
     131             : 
     132             : private:
     133             :   /**
     134             :    * Returns the point on the parameter mesh that is projected from the test point
     135             :    * @param p test point
     136             :    * @return Point
     137             :    */
     138             :   Point projectToMesh(const Point & p) const;
     139             : 
     140             :   /**
     141             :    * Find closest point on the element to the given point
     142             :    * @param elem
     143             :    * @param p
     144             :    * @return Point
     145             :    */
     146             :   Point closestPoint(const Elem & elem, const Point & p) const;
     147             : 
     148             :   /**
     149             :    * Compute regularization objective for a single quadrature point
     150             :    * This is the main function users should modify to add new regularization types for objectives
     151             :    * @param parameter_values  all parameter values
     152             :    * @param phi  shape function values (full array)
     153             :    * @param dphi  shape function gradients (full array)
     154             :    * @param qp  quadrature point index
     155             :    * @param dof_indices  element DOF indices
     156             :    * @param JxW  quadrature weights array
     157             :    * @param reg_type  type of regularization to compute
     158             :    * @return contribution to objective function
     159             :    */
     160             :   Real computeRegularizationQp(const std::vector<Real> & parameter_values,
     161             :                                const std::vector<std::vector<Real>> & phi,
     162             :                                const std::vector<std::vector<RealGradient>> & dphi,
     163             :                                const unsigned int qp,
     164             :                                const std::vector<dof_id_type> & dof_indices,
     165             :                                const std::vector<Real> & JxW,
     166             :                                RegularizationType reg_type) const;
     167             : 
     168             :   /**
     169             :    * Compute regularization gradient for a single quadrature point
     170             :    * This is the main function users should modify to add new regularization types for gradients
     171             :    * @param parameter_values  all parameter values
     172             :    * @param phi  shape function values (full array)
     173             :    * @param dphi  shape function gradients (full array)
     174             :    * @param qp  quadrature point index
     175             :    * @param dof_indices  element DOF indices
     176             :    * @param JxW  quadrature weights array
     177             :    * @param reg_type  type of regularization to compute
     178             :    * @param gradient  gradient vector to update
     179             :    */
     180             :   void computeRegularizationGradientQp(const std::vector<Real> & parameter_values,
     181             :                                        const std::vector<std::vector<Real>> & phi,
     182             :                                        const std::vector<std::vector<RealGradient>> & dphi,
     183             :                                        const unsigned int qp,
     184             :                                        const std::vector<dof_id_type> & dof_indices,
     185             :                                        const std::vector<Real> & JxW,
     186             :                                        RegularizationType reg_type,
     187             :                                        std::vector<Real> & gradient) const;
     188             :   // Cached values for gradient computations
     189             :   const unsigned short int _param_var_id;
     190             :   const libMesh::DofMap * _dof_map;
     191             :   const libMesh::FEType _fe_type;
     192             : };