nonlinear_functional.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/global/gate.hpp>
#include <kernel/global/vector.hpp>
 
namespace FEAT
{
  namespace Global
  {
    template<typename LocalNonlinearFunctional_, typename RowMirrorType_, typename ColMirrorType_>
    class NonlinearFunctional
    {
      public:
        typedef typename LocalNonlinearFunctional_::DataType DataType;
        typedef typename LocalNonlinearFunctional_::IndexType IndexType;
 
        typedef typename LocalNonlinearFunctional_::VectorTypeL LocalVectorTypeL;
        typedef typename LocalNonlinearFunctional_::VectorTypeR LocalVectorTypeR;
        typedef typename LocalNonlinearFunctional_::FilterType LocalFilterType;
 
        typedef Vector<LocalVectorTypeL, RowMirrorType_> VectorTypeL;
        typedef Vector<LocalVectorTypeR, ColMirrorType_> VectorTypeR;
        typedef Filter<LocalFilterType, RowMirrorType_> FilterType;
 
        typedef Gate<LocalVectorTypeL, RowMirrorType_> GateRowType;
        typedef Gate<LocalVectorTypeR, ColMirrorType_> GateColType;
 
        typedef VectorTypeL GradientType;
 
        static constexpr int BlockHeight = LocalNonlinearFunctional_::BlockHeight;
        static constexpr int BlockWidth = LocalNonlinearFunctional_::BlockWidth;
 
      protected:
        GateRowType* _row_gate;
        GateColType* _col_gate;
        LocalNonlinearFunctional_& _nonlinear_functional;
 
      public:
        explicit NonlinearFunctional(GateRowType* row_gate, GateColType* col_gate,
        LocalNonlinearFunctional_& nonlinear_functional) :
          _row_gate(row_gate),
          _col_gate(col_gate),
          _nonlinear_functional(nonlinear_functional)
          {
          }
 
        NonlinearFunctional() = delete;
        NonlinearFunctional(const NonlinearFunctional&) = delete;
 
        virtual ~NonlinearFunctional()
        {
        }
 
        LocalNonlinearFunctional_& local()
        {
          return _nonlinear_functional;
        }
 
        const LocalNonlinearFunctional_& local() const
        {
          return _nonlinear_functional;
        }
 
        const Dist::Comm* get_comm() const
        {
          return _row_gate != nullptr ? _row_gate->get_comm() : nullptr;
        }
 
        Index get_num_func_evals() const
        {
          return _nonlinear_functional.get_num_func_evals();
        }
 
        Index get_num_grad_evals() const
        {
          return _nonlinear_functional.get_num_grad_evals();
        }
 
        Index get_num_hess_evals() const
        {
          return _nonlinear_functional.get_num_hess_evals();
        }
 
        void reset_num_evals()
        {
          _nonlinear_functional.reset_num_evals();
        }
 
        VectorTypeL create_vector_l() const
        {
          return VectorTypeL(_row_gate, _nonlinear_functional.create_vector_l());
        }
 
        VectorTypeR create_vector_r() const
        {
          return VectorTypeR(_col_gate, _nonlinear_functional.create_vector_r());
        }
 
        void init()
        {
          _nonlinear_functional.init_pre_sync();
          // Synchronize the necessary parts
          while(!_nonlinear_functional.sync_scalars.empty())
          {
            auto it = _nonlinear_functional.sync_scalars.begin();
            if(_col_gate != nullptr)
            {
              *((*it).second) = _col_gate->sum(*((*it).second));
            }
            _nonlinear_functional.sync_scalars.erase(it);
          }
 
          if(_col_gate != nullptr)
          {
            for(auto* it:_nonlinear_functional.sync_vecs)
            {
              _col_gate->sync_0(*it);
            }
          }
          _nonlinear_functional.init_post_sync();
        }
 
        void prepare(const VectorTypeR& vec_state, FilterType& filter)
        {
 
          // Prepare the patch local nonlinear functional, computing everything that does need synchronising
          _nonlinear_functional.prepare_pre_sync(vec_state.local(), filter.local());
 
          // Synchronize the necessary parts
          while(!_nonlinear_functional.sync_scalars.empty())
          {
            auto it = _nonlinear_functional.sync_scalars.begin();
            if(_col_gate != nullptr)
            {
              *((*it).second) = _col_gate->sum(*((*it).second));
            }
            _nonlinear_functional.sync_scalars.erase(it);
          }
 
          if(_col_gate != nullptr)
          {
            for(auto* it:_nonlinear_functional.sync_vecs)
            {
              _col_gate->sync_0(*it);
            }
          }
 
          // Prepare the rest of the functional
          _nonlinear_functional.prepare_post_sync(vec_state.local(), filter.local());
 
          // Sync the filter vector in the SlipFilter
          if(_col_gate != nullptr )
          {
            // For all slip filters...
            for(auto& it : filter.local().template at<0>())
            {
              // get the filter vector
              auto& slip_filter_vector = it.second.get_filter_vector();
 
              if(slip_filter_vector.used_elements() > 0)
              {
                // Temporary DenseVector for syncing
                LocalVectorTypeR tmp(slip_filter_vector.size(), DataType(0));
 
                auto* tmp_elements = tmp.template elements<LAFEM::Perspective::native>();
                auto* sfv_elements = slip_filter_vector.template elements<LAFEM::Perspective::native>();
 
                // Copy sparse filter vector contents to DenseVector
                for(Index isparse(0); isparse < slip_filter_vector.used_elements(); ++isparse)
                {
                  Index idense(slip_filter_vector.indices()[isparse]);
                  tmp_elements[idense] = sfv_elements[isparse];
                }
 
                _col_gate->sync_0(tmp);
                // Copy sparse filter vector contents to DenseVector
                for(Index isparse(0); isparse < slip_filter_vector.used_elements(); ++isparse)
                {
                  Index idense(slip_filter_vector.indices()[isparse]);
                  tmp_elements[idense].normalize();
                  sfv_elements[isparse] = tmp_elements[idense];
 
                }
              }
 
              else
              {
                // Temporary DenseVector for syncing
                LocalVectorTypeR tmp(slip_filter_vector.size(), DataType(0));
                _col_gate->sync_0(tmp);
              }
            }
          } // col_gate
        } // prepare()
 
        Index columns()
        {
          // Compute total number of rows and columns
          auto vec_r = create_vector_r();
          vec_r.format(DataType(1));
 
          return Index(vec_r.norm2sqr());
        }
 
        Index rows()
        {
          // Compute total number of rows and rows
          auto vec_l = create_vector_l();
          vec_l.format(DataType(1));
 
          return Index(vec_l.norm2sqr());
        }
 
        void eval_fval_grad(DataType& fval, VectorTypeL& grad)
        {
          // As the penalty term is quadratic, we have to compute the local functional values separately for summing
          // them up over all ranks
          const bool add_penalty_fval(false);
          _nonlinear_functional.eval_fval_grad(fval, grad.local(), add_penalty_fval);
 
          // Sum up over all patches
          if(_row_gate != nullptr)
            fval = _row_gate->sum(fval);
 
          // Add the penalty term
          if(get_penalty_param() > DataType(0))
          {
            DataType constraint(get_constraint());
            fval += get_penalty_param()*DataType(0.5)*Math::sqr(constraint);
          }
 
          // Synchronize the gradient vector
          grad.sync_0();
        }
 
        DataType get_penalty_param() const
        {
          return _nonlinear_functional.get_penalty_param();
        }
 
        void set_penalty_param(const DataType fac)
        {
          _nonlinear_functional.set_penalty_param(fac);
        }
 
        DataType get_constraint()
        {
          return _nonlinear_functional.get_constraint();
        }
 
        DataType compute_constraint()
        {
          return _nonlinear_functional.compute_constraint();
        }
 
    };
  } // namespace Global
} // namespace FEAT