ProblemStat.inc.hpp

#pragma once


#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>

#include "AdaptInfo.hpp"
#include "Loops.hpp"
#include "Timer.hpp"

namespace AMDiS
{
    template <class Param>
    void ProblemStat<Param>::initialize(Flag initFlag,
					ProblemStat* adoptProblem,
					Flag adoptFlag) 
    {
	// === create meshes ===
	if (mesh) {
	    AMDIS_WARNING("mesh already created");
	}
	else {
	    if (initFlag.isSet(CREATE_MESH) ||
		(!adoptFlag.isSet(INIT_MESH) &&
		(initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
		createMesh();
	    }

	    if (adoptProblem &&
		(adoptFlag.isSet(INIT_MESH) ||
		adoptFlag.isSet(INIT_SYSTEM) ||
		adoptFlag.isSet(INIT_FE_SPACE))) {
		mesh = adoptProblem->getMesh();
	    }
	    
	    componentMeshes.resize(nComponents, mesh.get());
	}

	if (!mesh) {
	    AMDIS_WARNING("no mesh created");
	}
	
	
	// === create fespace ===
	if (feSpaces) {
	    AMDIS_WARNING("feSpaces already created");
	}
	else {
	    if (initFlag.isSet(INIT_FE_SPACE) ||
		(initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
		createFeSpaces();
	    }

	    if (adoptProblem &&
		(adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
		feSpaces = adoptProblem->getFeSpaces();
	    }
	}
	
	
	int globalRefinements = 0;
	Parameters::get(meshName + "->global refinements", globalRefinements);
	mesh->globalRefine(globalRefinements);
	
	solution->compress();
	rhs->compress();
    }
  
  
    template <class Param>
    void ProblemStat<Param>::addMatrixOperator(OperatorType& op, 
					       int i, int j,
					       double* factor, 
					       double* estFactor)
    {
	// TODO: currently the factors are ignored
	assert( factor == NULL && estFactor == NULL );
	
	matrixOperators[std::make_pair(i,j)].push_back(std::make_shared<OperatorType>(op));
    }

    
    template <class Param>
    void ProblemStat<Param>::addVectorOperator(OperatorType& op, 
					       int i,
					       double* factor, 
					       double* estFactor)
    {
	// TODO: currently the factors are ignored
	assert( factor == NULL && estFactor == NULL );
	
	vectorOperators[i].push_back(std::make_shared<OperatorType>(op));
    }

  
  
  /// Adds a Dirichlet boundary condition
  template <class Param>
    template <class Predicate, class Values>
  void ProblemStat<Param>::addDirichletBC(Predicate const& predicate, 
					  int row, int col, 
					  Values const& values)
  {
    static_assert(Dune::Functions::Concept::isCallable<Predicate, WorldVector>(), 
      "Function passed to addDirichletBC for predicate does not model the Callable<WorldVector> concept");
    
    static_assert(Dune::Functions::Concept::isCallable<Values, WorldVector>(), 
      "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

    AMDIS_TEST_EXIT(row >= 0 && row < nComponents, "row number out of range: " << row);
    AMDIS_TEST_EXIT(col >= 0 && col < nComponents, "col number out of range: " << col);
    
//     boundaryConditionSet = true;

    using BcType = DirichletBC<WorldVector>;
    std::shared_ptr<BcType> dirichlet = std::make_shared<BcType>(predicate, values);
      
    dirichletBc[std::make_pair(row, col)].push_back( dirichlet );
  }
  
  
  template <class Param>
  void ProblemStat<Param>::solve(AdaptInfo& adaptInfo, 
				 bool createMatrixData, 
				 bool storeMatrixData)
  {
//     if (!solver) {
//       WARNING("no solver\n");
//       return;
//     }

    Timer t;
    
    // NOTE: just to show the priciple, with initial solution 0
//     For<0, nComponents>::loop([this](auto const _i) {
// 	solution->getVector(_i) = 0.0;
//     });
    
    using VectorType = std::tuple_element_t<0, typename SystemVectorType::BaseVectors>;
    
    // Technicality : turn the matrix into a linear operator
    Dune::MatrixAdapter<DOFMatrix, VectorType, VectorType> op(systemMatrix[0][0]);

    // Sequential incomplete LU decomposition as the preconditioner
    Dune::SeqILU0<DOFMatrix, VectorType, VectorType> ilu0(systemMatrix[0][0], 1.0);

    // Preconditioned conjugate−gradient solver
    Dune::CGSolver<VectorType> cg(op, ilu0, 1.e-4, 500, 2);

    // storing some statistics
    Dune::InverseOperatorResult statistics;

    // solve the linear system
    constexpr index_<0> _0;
    cg.apply(solution->getVector(_0), rhs->getVector(_0), statistics);
    
    AMDIS_MSG("solution of discrete system needed " << t.elapsed() << " seconds");

    adaptInfo.setSolverIterations(statistics.iterations);
//     adaptInfo.setMaxSolverIterations(solver->getMaxIterations());
//     adaptInfo.setSolverTolerance(solver->getTolerance());
    adaptInfo.setSolverResidual(statistics.reduction);
  }
  
  
  template <class Param>
  void ProblemStat<Param>::buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag,
					     bool asmMatrix, bool asmVector)
  {
    Timer t;
    solution->compress();
      

    For<0, nComponents>::loop([this, asmMatrix, asmVector](auto const _r) 
    {
      if (asmVector) {
	rhs->compress(_r);
        rhs->getVector(_r) = 0.0;
      }

      For<0, nComponents>::loop([this, asmMatrix, asmVector, _r](auto const _c) 
      {
	auto const& rowFeSpace = this->getFeSpace(_r);
        auto const& colFeSpace = this->getFeSpace(_c);
	
	int r = int(_r), c = int(_c);
	auto row_col = std::make_pair(r, c);
	
        // The DOFMatrix which should be assembled
        DOFMatrix& matrix = systemMatrix[_r][_c];
	  
	if (asmMatrix) {
	  // TODO: calc occupationPattern only once for each feSpace
	  Dune::MatrixIndexSet occupationPattern;
	  this->getOccupationPattern(rowFeSpace, colFeSpace, occupationPattern);
	  occupationPattern.exportIdx(matrix);

	  // TODO: only set to zero if not keep the matrix
	  matrix = 0.0;

	  // init boundary condition
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
		for (auto bc : bc_list.second)
		    bc->init(c == int(_c), 
			     rowFeSpace, colFeSpace, &matrix, 
			     &rhs->getVector(_r), 
			     &solution->getVector(_c));
	    }
	  }
	  
	  this->assemble(row_col, rowFeSpace, colFeSpace, &matrix,
		  (_r == _c && asmVector ? &rhs->getVector(_r) : NULL));

	  // finish boundary condition
	  for (auto bc_list : dirichletBc) {
	    std::tie(r, c) = bc_list.first;
	    if (r == int(_r)) {
		for (auto bc : bc_list.second)
		    bc->finish(c == int(_c), 
			       rowFeSpace, colFeSpace, &matrix, 
			       &rhs->getVector(_r), 
			       &solution->getVector(_c));
	    }
	  }
	}
      });
    });
    
    AMDIS_MSG("buildAfterCoarsen needed " << t.elapsed() << " seconds");
  }
  

  template <class Param>
  void ProblemStat<Param>::writeFiles(AdaptInfo& adaptInfo, double time)
  {
    Timer t;
    
    Dune::VTKWriter<MeshView> vtkWriter(*meshView);
    For<0, nComponents>::loop([this, &vtkWriter](const auto _i) {
	vtkWriter.addVertexData(solution->getVector(_i), componentNames[_i]);
    });
    
    std::string filename = "";
    Parameters::get(name + "->output->filename", filename);
    if (time >= 0.0)
	filename += std::to_string(time);
    
    vtkWriter.write(filename);


    AMDIS_MSG("writeFiles needed " << t.elapsed() << " seconds");
  }
  
  
  
  // Get the occupation pattern
  template <class Param>
    template <class RowFeSpace, class ColFeSpace>
  void ProblemStat<Param>::getOccupationPattern(const RowFeSpace& rowFeSpace, 
						const ColFeSpace& colFeSpace, 
						Dune::MatrixIndexSet& nb)
  {
      nb.resize(rowFeSpace.size(), colFeSpace.size());

      // A loop over all elements of the grid
      auto rowView = rowFeSpace.localView();
      auto rowIndexSet = rowFeSpace.localIndexSet();
      
      auto colView = colFeSpace.localView();
      auto colIndexSet = colFeSpace.localIndexSet();

      for (const auto& element : elements(*meshView)) {
	  rowView.bind(element);
	  rowIndexSet.bind(rowView); // TODO: can this be moved out of the loop?
	  
	  colView.bind(element);
	  colIndexSet.bind(colView); // TODO: can this be moved out of the loop?
	  
	  for (size_t i = 0; i < rowIndexSet.size(); ++i) {
	      // The global index of the i-th vertex of the element
	      auto row = rowIndexSet.index(i);
	      for (size_t j = 0; j < colIndexSet.size(); ++j) {
		  // The global index of the j-th vertex of the element
		  auto col = colIndexSet.index(j);
		  nb.add(row, col);
	      }
	  }
      }
  }


  template <class Param>
    template <class RowFeSpace, class ColFeSpace, class Vector>
  void ProblemStat<Param>::assemble(std::pair<int, int> row_col,
				    const RowFeSpace& rowFeSpace, 
				    const ColFeSpace& colFeSpace, 
				    DOFMatrix* matrix,
				    Vector* rhs)
  {
      for (auto op : matrixOperators[row_col])
	op->init(rowFeSpace, colFeSpace);
      for (auto op : vectorOperators[row_col.first])
	op->init(rowFeSpace, colFeSpace);
	      
      auto rowView = rowFeSpace.localView();
      auto rowIndexSet = rowFeSpace.localIndexSet();
      
      auto colView = colFeSpace.localView();
      auto colIndexSet = colFeSpace.localIndexSet();
      
      for (const auto& element : elements(*meshView)) {
	  rowView.bind(element);
	  rowIndexSet.bind(rowView);
	  
	  colView.bind(element);
	  colIndexSet.bind(colView);
	  
	  if (matrix) {	      
	      ElementMatrix elementMatrix;

	      // DOFMatrix::assemble {
	      bool add = getElementMatrix(rowView, colView, elementMatrix, matrixOperators[row_col]);

	      // DOFMatrix::addElementMatrix {
	      if (add) {
		  for (size_t i = 0; i < elementMatrix.N(); ++i) {
		      // The global index of the i−th vertex of the element
		      const auto row = rowIndexSet.index(i);
		      for (size_t j = 0; j < elementMatrix.M(); ++j) {
			  // The global index of the j−th vertex of the element
			  const auto col = colIndexSet.index(j);
			  (*matrix)[row][col] += elementMatrix[i][j];
		      }
		  }
	      }
	      // }
	      // }
	  }
	  
	  if (rhs) {
	      ElementVector elementVector;
	      
	      // DOFVectorBase<double>::assemble {
	      bool add = getElementVector(rowView, elementVector, vectorOperators[row_col.first]);
	      
	      // DOFVectorBase<T>::addElementVector {
	      if (add) {
		  for (size_t i = 0; i < elementVector.size(); ++i) {
		      // The global index of the i-th vertex
		      const auto row = rowIndexSet.index(i);
		      (*rhs)[row] += elementVector[i];
		  }
	      }
	      // }
	      // }
	  }
      }
  }
  
  
  template <class Param>
    template <class RowView, class ColView>
  bool ProblemStat<Param>::getElementMatrix(RowView const& rowView,
					    ColView const& colView,
					    ElementMatrix& elementMatrix,
					    std::list<std::shared_ptr<OperatorType>>& operators)
  {
    const auto nRows = rowView.tree().finiteElement().size();
    const auto nCols = colView.tree().finiteElement().size();

    elementMatrix.setSize(nRows, nCols);
    elementMatrix = 0.0; // fills the entire matrix with zeroes

    for (auto op : operators)
      op->getElementMatrix(rowView, colView, elementMatrix);
    
    return !operators.empty();
  }
  
  
  
  template <class Param>
    template <class RowView>
  bool ProblemStat<Param>::getElementVector(RowView const& rowView,
					    ElementVector& elementVector,
					    std::list<std::shared_ptr<OperatorType>>& operators)
  {
    const auto nRows = rowView.tree().finiteElement().size();

    // Set all entries to zero
    elementVector.resize(nRows);
    elementVector = 0.0;

    for (auto op : operators)
      op->getElementVector(rowView, elementVector);
    
    return !operators.empty();
  }

} // end namespace AMDiS