statpdeproblem.h

Go to the documentation of this file.

#ifndef STAT_PDE_PROBLEM_H_
#define STAT_PDE_PROBLEM_H_

#include "pdeprobleminterface.h"
#include "integrator.h"
#include "parameterreader.h"
#include "statevector.h"
#include "stateproblem.h"

#include <deal.II/lac/vector.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_sparse_matrix.h>

#include "pdeproblemcontainer.h"
#include "pdeinterface.h"
#include "functionalinterface.h"
#include "dirichletdatainterface.h"
#include "dopeexception.h"
#include "newtonsolver.h"
#include "cglinearsolver.h"
#include "gmreslinearsolver.h"
#include "directlinearsolver.h"
#include "solutionextractor.h"

#include <deal.II/base/data_out_base.h>
#include <deal.II/base/utilities.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/function.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/sparse_direct.h>

#include <fstream>
namespace DOpE
{
  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    class StatPDEProblem : public PDEProblemInterface<PROBLEM, VECTOR, dealdim>
    {
      public:
        template<typename INTEGRATORDATACONT>
          StatPDEProblem(PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
              ParameterReader &param_reader, INTEGRATORDATACONT& idc,
              int base_priority = 0);

        template<typename INTEGRATORDATACONT>
          StatPDEProblem(PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
              ParameterReader &param_reader, INTEGRATORDATACONT& idc1,
              INTEGRATORDATACONT& idc2, int base_priority = 0);

        virtual
        ~StatPDEProblem();

        /******************************************************/

        static void
        declare_params(ParameterReader &param_reader);

        /******************************************************/

        void
        ReInit();

        /******************************************************/

        void
        ComputeReducedFunctionals();

        /******************************************************/

        void
        ComputeReducedFunctionals(const StateVector<VECTOR>& st_u);

        /******************************************************/

        template<class DWRC, class PDE>
          void
          ComputeRefinementIndicators(DWRC& dwrc, PDE& pde);

        /******************************************************/

        void
        StateSizeInfo(std::stringstream& out)
        {
          GetU().PrintInfos(out);
        }

        /******************************************************/

        void
        WriteToFile(const VECTOR &v, std::string name, std::string outfile,
            std::string dof_type, std::string filetype);

        /******************************************************/

        void
        WriteToFileElementwise(const Vector<double> &v, std::string name,
            std::string outfile, std::string dof_type, std::string filetype);

        /******************************************************/

        void
        WriteToFile(const ControlVector<VECTOR> &v, std::string name, std::string dof_type);

        void
          WriteToFile(const std::vector<double> &/*v*/,
            std::string /*outfile*/)
        {
          abort();
        }

      protected:
        void
        ComputeReducedState();

        /******************************************************/

        void
        ComputeDualForErrorEstimation(DOpEtypes::WeightComputation);

        /******************************************************/
        const StateVector<VECTOR> &
        GetU() const
        {
          return u_;
        }
        StateVector<VECTOR> &
        GetU()
        {
          return u_;
        }
        const StateVector<VECTOR> &
        GetZForEE() const
        {
          return z_for_ee_;
        }
        StateVector<VECTOR> &
        GetZForEE()
        {
          return z_for_ee_;
        }

        NONLINEARSOLVER&
        GetNonlinearSolver(std::string type);
        INTEGRATOR&
        GetIntegrator()
        {
          return integrator_;
        }

        const bool&
        GetBuildStateMatrix() const
        {
          return build_state_matrix_;
        }

        const bool&
        GetBuildAdjointMatrix() const
        {
          return build_adjoint_matrix_;
        }

        const bool&
        GetStateReinit() const
        {
          return state_reinit_;
        }

        const bool&
        GetAdjointReinit() const
        {
          return adjoint_reinit_;
        }


        bool&
        GetBuildStateMatrix()
        {
          return build_state_matrix_;
        }

        bool&
        GetBuildAdjointMatrix()
        {
          return build_adjoint_matrix_;
        }

        bool&
        GetStateReinit()
        {
          return state_reinit_;
        }

        bool&
        GetAdjointReinit()
        {
          return adjoint_reinit_;
        }


      private:
        void
        AddUDD()
        {
          for (auto it = this->GetUserDomainData().begin();
              it != this->GetUserDomainData().end(); it++)
              {
            this->GetIntegrator().AddDomainData(it->first, it->second);
          }
        }

        void
        DeleteUDD()
        {
          for (auto it = this->GetUserDomainData().begin();
              it != this->GetUserDomainData().end(); it++)
              {
            this->GetIntegrator().DeleteDomainData(it->first);
          }
        }

        StateVector<VECTOR> u_;
        StateVector<VECTOR> z_for_ee_;

        INTEGRATOR integrator_;
        NONLINEARSOLVER nonlinear_state_solver_;
        NONLINEARSOLVER nonlinear_adjoint_solver_;

        bool build_state_matrix_;
        bool build_adjoint_matrix_;
        bool state_reinit_;
        bool adjoint_reinit_;

        int n_patches_;

        friend class SolutionExtractor<
            StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>,
            VECTOR> ;
    };

  /*************************************************************************/
  /*****************************IMPLEMENTATION******************************/
  /*************************************************************************/
  using namespace dealii;

  /******************************************************/
  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::declare_params(
        ParameterReader &param_reader)
    {
      NONLINEARSOLVER::declare_params(param_reader);
      param_reader.SetSubsection("output parameters");
      param_reader.declare_entry("number of patches", "0",
          Patterns::Integer(0));

    }
  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    template<typename INTEGRATORDATACONT>
      StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::StatPDEProblem(
          PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
          ParameterReader &param_reader, INTEGRATORDATACONT& idc,
          int base_priority)
          : PDEProblemInterface<PROBLEM, VECTOR, dealdim>(OP, base_priority), u_(
              OP->GetSpaceTimeHandler(), state_behavior, param_reader), z_for_ee_(
              OP->GetSpaceTimeHandler(), state_behavior, param_reader), integrator_(
              idc), nonlinear_state_solver_(integrator_, param_reader), nonlinear_adjoint_solver_(
              integrator_, param_reader)
      {
        param_reader.SetSubsection("output parameters");
        n_patches_ = param_reader.get_integer("number of patches");
        //PDEProblems should be ReInited
        {
          state_reinit_ = true;
          adjoint_reinit_ = true;
        }
      }

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    template<typename INTEGRATORDATACONT>
      StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::StatPDEProblem(
          PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
          ParameterReader &param_reader, INTEGRATORDATACONT& idc,
          INTEGRATORDATACONT& idc2, int base_priority)
          : PDEProblemInterface<PROBLEM, VECTOR, dealdim>(OP, base_priority), u_(
              OP->GetSpaceTimeHandler(), state_behavior, param_reader), z_for_ee_(
              OP->GetSpaceTimeHandler(), state_behavior, param_reader), integrator_(
              idc, idc2), nonlinear_state_solver_(integrator_, param_reader), nonlinear_adjoint_solver_(
              integrator_, param_reader)
      {
        //PDEProblems should be ReInited
        {
          state_reinit_ = true;
          adjoint_reinit_ = true;
        }
      }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::~StatPDEProblem()
    {
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    NONLINEARSOLVER&
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::GetNonlinearSolver(
        std::string type)
    {
      if (type == "state")
      {
        return nonlinear_state_solver_;
      }
      else if (type == "adjoint_for_ee")
      {
        return nonlinear_adjoint_solver_;
      }
      else
      {
        throw DOpEException("No Solver for Problem type:`" + type + "' found",
            "StatReducedProblem::GetNonlinearSolver");

      }
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::ReInit()
    {
      PDEProblemInterface < PROBLEM, VECTOR, dealdim > ::ReInit();

      //Some Solvers must be reinited when called
      // Better we have subproblems, so that solver can be reinited here
      {
        state_reinit_ = true;
        adjoint_reinit_ = true;
      }

      build_state_matrix_ = true;
      build_adjoint_matrix_ = true;

      GetU().ReInit();
      GetZForEE().ReInit();
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::ComputeReducedState()
    {
      this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals());
      this->GetOutputHandler()->Write("Computing State Solution:",
          4 + this->GetBasePriority());

      this->SetProblemType("state");
      auto& problem = this->GetProblem()->GetStateProblem();
      if (state_reinit_ == true)
      {
        GetNonlinearSolver("state").ReInit(problem);
        state_reinit_ = false;
      }

      this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
      AddUDD();
      build_state_matrix_ = this->GetNonlinearSolver("state").NonlinearSolve(
          problem, (GetU().GetSpacialVector()), true, build_state_matrix_);
      DeleteUDD();
      
      this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());

      this->GetOutputHandler()->Write((GetU().GetSpacialVector()),
      "State" + this->GetPostIndex(), problem.GetDoFType());

    }

      /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::ComputeDualForErrorEstimation(
        DOpEtypes::WeightComputation weight_comp)
    {
      this->GetOutputHandler()->Write("Computing Dual for Error Estimation:",
          4 + this->GetBasePriority());
      if (weight_comp == DOpEtypes::higher_order_interpolation)
      {
        this->SetProblemType("adjoint_for_ee");
      }
      else
      {
        throw DOpEException("Unknown WeightComputation",
            "StatPDEProblem::ComputeDualForErrorEstimation");
      }

      //      auto& problem = this->GetProblem()->GetStateProblem();//Hier ist adjoint problem einzufuegen
      auto& problem = *(this->GetProblem());
      if (adjoint_reinit_ == true)
      {
        GetNonlinearSolver("adjoint_for_ee").ReInit(problem);
        adjoint_reinit_ = false;
      }

      this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());

      this->GetIntegrator().AddDomainData("state",
          &(GetU().GetSpacialVector()));
      AddUDD();

      build_adjoint_matrix_ =
          this->GetNonlinearSolver("adjoint_for_ee").NonlinearSolve(problem,
              (GetZForEE().GetSpacialVector()), true, build_adjoint_matrix_);

      this->GetIntegrator().DeleteDomainData("state");
      DeleteUDD();

      this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());

      this->GetOutputHandler()->Write((GetZForEE().GetSpacialVector()),
      "Adjoint_for_ee" + this->GetPostIndex(), problem.GetDoFType());

    }

      /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::ComputeReducedFunctionals()
    {
      this->ComputeReducedState();

      this->GetOutputHandler()->Write("Computing Functionals:",
          4 + this->GetBasePriority());

      this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());

      this->GetIntegrator().AddDomainData("state",
          &(GetU().GetSpacialVector()));
      AddUDD();

      for (unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
      {
        double ret = 0;
        bool found = false;

        this->SetProblemType("aux_functional", i);
        if (this->GetProblem()->GetFunctionalType().find("domain")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputeDomainScalar(
              *(this->GetProblem()));
        }
        if (this->GetProblem()->GetFunctionalType().find("point")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputePointScalar(
              *(this->GetProblem()));
        }
        if (this->GetProblem()->GetFunctionalType().find("boundary")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputeBoundaryScalar(
              *(this->GetProblem()));
        }
        if (this->GetProblem()->GetFunctionalType().find("face")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
        }

        if (!found)
        {
          throw DOpEException(
              "Unknown Functional Type: "
                  + this->GetProblem()->GetFunctionalType(),
              "StatPDEProblem::ComputeReducedFunctionals");
        }
        this->GetFunctionalValues()[i].push_back(ret);
        std::stringstream out;
        this->GetOutputHandler()->InitOut(out);
        out << this->GetProblem()->GetFunctionalName() << ": " << ret;
        this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
      }

      this->GetIntegrator().DeleteDomainData("state");
      DeleteUDD();

      this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::ComputeReducedFunctionals(
        const StateVector<VECTOR>& st_u)
    {
      this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals());
      GetU() = st_u;

      this->GetOutputHandler()->Write("Computing Functionals:",
          4 + this->GetBasePriority());

      this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());

      this->GetIntegrator().AddDomainData("state",
          &(GetU().GetSpacialVector()));
      AddUDD();

      for (unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
      {
        double ret = 0;
        bool found = false;

        this->SetProblemType("aux_functional", i);
        if (this->GetProblem()->GetFunctionalType().find("domain")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputeDomainScalar(
              *(this->GetProblem()));
        }
        if (this->GetProblem()->GetFunctionalType().find("point")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputePointScalar(
              *(this->GetProblem()));
        }
        if (this->GetProblem()->GetFunctionalType().find("boundary")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputeBoundaryScalar(
              *(this->GetProblem()));
        }
        if (this->GetProblem()->GetFunctionalType().find("face")
            != std::string::npos)
        {
          found = true;
          ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
        }

        if (!found)
        {
          throw DOpEException(
              "Unknown Functional Type: "
                  + this->GetProblem()->GetFunctionalType(),
              "StatPDEProblem::ComputeReducedFunctionals");
        }
        this->GetFunctionalValues()[i].push_back(ret);
        std::stringstream out;
        this->GetOutputHandler()->InitOut(out);
        out << this->GetProblem()->GetFunctionalName() << ": " << ret;
        this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
      }

      this->GetIntegrator().DeleteDomainData("state");
      DeleteUDD();

      this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    template<class DWRC, class PDE>
      void
      StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::ComputeRefinementIndicators(
          DWRC& dwrc, PDE& pde)
      {
        //Attach the ResidualModifier to the PDE.
        pde.ResidualModifier = boost::bind<void>(
            boost::mem_fn(&DWRC::ResidualModifier), boost::ref(dwrc), _1);
        pde.VectorResidualModifier = boost::bind<void>(
            boost::mem_fn(&DWRC::VectorResidualModifier), boost::ref(dwrc), _1);
        //first we reinit the dwrdatacontainer (this
        //sets the weight-vectors to their correct length)
        const unsigned int n_elements =
            this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler().get_tria().n_active_cells();
        dwrc.ReInit(n_elements);
        //If we need the dual solution, compute it
        if (dwrc.NeedDual())
          this->ComputeDualForErrorEstimation(dwrc.GetWeightComputation());

        //some output
        this->GetOutputHandler()->Write("Computing Error Indicators:",
            4 + this->GetBasePriority());

        this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());

        //add the primal and (if needed) dual solution to the integrator
        this->GetIntegrator().AddDomainData("state",
            &(GetU().GetSpacialVector()));
        if (dwrc.NeedDual())
          this->GetIntegrator().AddDomainData("adjoint_for_ee",
              &(GetZForEE().GetSpacialVector()));
        AddUDD();

        this->SetProblemType("error_evaluation");

        //prepare the weights...
        dwrc.PrepareWeights(GetU(), GetZForEE());

        //now we finally compute the refinement indicators
        this->GetIntegrator().ComputeRefinementIndicators(*this->GetProblem(),
            dwrc);

        // release the lock on the refinement indicators (see dwrcontainer.h)
        dwrc.ReleaseLock();

        const float error = dwrc.GetError();

        // clear the data
        dwrc.ClearWeightData();
        this->GetIntegrator().DeleteDomainData("state");
        if (dwrc.NeedDual())
          this->GetIntegrator().DeleteDomainData("adjoint_for_ee");
        DeleteUDD();
        this->GetProblem()->DeleteAuxiliaryFromIntegrator(
            this->GetIntegrator());

        std::stringstream out;
        this->GetOutputHandler()->InitOut(out);
        out << "Error estimate using " << dwrc.GetName();
        if (dwrc.NeedDual())
          out << " for the " << this->GetProblem()->GetFunctionalName();
        out << ": " << error;
        this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
        this->GetOutputHandler()->WriteElementwise(dwrc.GetErrorIndicators(),
            "Error_Indicators" + this->GetPostIndex(),
            this->GetProblem()->GetDoFType());

      }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::WriteToFileElementwise(
        const Vector<double> &v, std::string name, std::string outfile,
        std::string dof_type, std::string filetype)
    {
      if (dof_type == "state")
      {
        auto& data_out =
            this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
        data_out.attach_dof_handler(
            this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());

        data_out.add_data_vector(v, name);
        data_out.build_patches(n_patches_);

        std::ofstream output(outfile.c_str());

        if (filetype == ".vtk")
        {
          data_out.write_vtk(output);
        }
        else
        {
          throw DOpEException(
              "Don't know how to write filetype `" + filetype + "'!",
              "StatPDEProblem::WriteToFileElementwise");
        }
        data_out.clear();
      }
      else
      {
        throw DOpEException("No such DoFHandler `" + dof_type + "'!",
            "StatPDEProblem::WriteToFileElementwise");
      }
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::WriteToFile(
        const VECTOR &v, std::string name, std::string outfile,
        std::string dof_type, std::string filetype)
    {
      if (dof_type == "state")
      {
        auto& data_out =
            this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
        data_out.attach_dof_handler(
            this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());

        data_out.add_data_vector(v, name);
        //TODO: mapping[0] is a workaround, as deal does not support interpolate
        // boundary_values with a mapping collection at this point.
        data_out.build_patches(
            this->GetProblem()->GetSpaceTimeHandler()->GetMapping()[0]);

        std::ofstream output(outfile.c_str());

        if (filetype == ".vtk")
        {
          data_out.write_vtk(output);
        }
        else if (filetype == ".gpl")
        {
          data_out.write_gnuplot(output);
        }
        else
        {
          throw DOpEException(
              "Don't know how to write filetype `" + filetype + "'!",
              "StatPDEProblem::WriteToFile");
        }
        data_out.clear();
      }
      else
      {
        throw DOpEException("No such DoFHandler `" + dof_type + "'!",
            "StatPDEProblem::WriteToFile");
      }
    }

  /******************************************************/

  template<typename NONLINEARSOLVER, typename INTEGRATOR, typename PROBLEM,
      typename VECTOR, int dealdim>
    void
    StatPDEProblem<NONLINEARSOLVER, INTEGRATOR, PROBLEM, VECTOR, dealdim>::WriteToFile(
      const ControlVector<VECTOR> &/*v*/, std::string /*name*/, std::string /*dof_type*/)
    {
      throw DOpEException("This Problem does not support ControlVectors","StatPDEProblem::WriteToFile");
    }
}
#endif