MultiAppProjectionTransfer.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "MultiAppProjectionTransfer.h"

// MOOSE includes
#include "AddVariableAction.h"
#include "FEProblem.h"
#include "MooseMesh.h"
#include "MooseVariableFE.h"
#include "SystemBase.h"

#include "libmesh/dof_map.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/mesh_function.h"
#include "libmesh/mesh_tools.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/parallel_algebra.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/default_coupling.h"

void
assemble_l2(EquationSystems & es, const std::string & system_name)
{
  MultiAppProjectionTransfer * transfer =
      es.parameters.get<MultiAppProjectionTransfer *>("transfer");
  transfer->assembleL2(es, system_name);
}

registerMooseObject("MooseApp", MultiAppProjectionTransfer);

template <>
InputParameters
validParams<MultiAppProjectionTransfer>()
{
  InputParameters params = validParams<MultiAppFieldTransfer>();

  MooseEnum proj_type("l2", "l2");
  params.addParam<MooseEnum>("proj_type", proj_type, "The type of the projection.");

  params.addParam<bool>("fixed_meshes",
                        false,
                        "Set to true when the meshes are not changing (ie, "
                        "no movement or adaptivity).  This will cache some "
                        "information to speed up the transfer.");

  // Need one layer of ghosting
  params.addRelationshipManager("ElementSideNeighborLayers",
                                Moose::RelationshipManagerType::GEOMETRIC |
                                    Moose::RelationshipManagerType::ALGEBRAIC);
  return params;
}

MultiAppProjectionTransfer::MultiAppProjectionTransfer(const InputParameters & parameters)
  : MultiAppFieldTransfer(parameters),
    _proj_type(getParam<MooseEnum>("proj_type")),
    _compute_matrix(true),
    _fixed_meshes(getParam<bool>("fixed_meshes")),
    _qps_cached(false)
{
  if (_to_var_names.size() != 1 || _from_var_names.size() != 1)
    mooseError(" Support single variable only ");
}

void
MultiAppProjectionTransfer::initialSetup()
{
  MultiAppFieldTransfer::initialSetup();

  getAppInfo();

  _proj_sys.resize(_to_problems.size(), NULL);

  for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
  {
    FEProblemBase & to_problem = *_to_problems[i_to];
    EquationSystems & to_es = to_problem.es();

    // Add the projection system.
    FEType fe_type = to_problem
                         .getVariable(0,
                                      _to_var_name,
                                      Moose::VarKindType::VAR_ANY,
                                      Moose::VarFieldType::VAR_FIELD_STANDARD)
                         .feType();

    LinearImplicitSystem & proj_sys = to_es.add_system<LinearImplicitSystem>("proj-sys-" + name());

    proj_sys.get_dof_map().add_coupling_functor(
        proj_sys.get_dof_map().default_coupling(),
        false); // The false keeps it from getting added to the mesh

    _proj_var_num = proj_sys.add_variable("var", fe_type);
    proj_sys.attach_assemble_function(assemble_l2);
    _proj_sys[i_to] = &proj_sys;

    // Prevent the projection system from being written to checkpoint
    // files.  In the event of a recover or restart, we'll read the checkpoint
    // before this initialSetup method is called.  As a result, we'll find
    // systems in the checkpoint (the projection systems) that we don't know
    // what to do with, and there will be a crash.  We could fix this by making
    // the systems in the constructor, except we don't know how many sub apps
    // there are at the time of construction.  So instead, we'll just nuke the
    // projection system and rebuild it from scratch every recover/restart.
    proj_sys.hide_output() = true;

    // Reinitialize EquationSystems since we added a system.
    to_es.reinit();
  }

  if (_fixed_meshes)
  {
    _cached_qps.resize(n_processors());
    _cached_index_map.resize(n_processors());
  }
}

void
MultiAppProjectionTransfer::assembleL2(EquationSystems & es, const std::string & system_name)
{
  // Get the system and mesh from the input arguments.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(system_name);
  MeshBase & to_mesh = es.get_mesh();

  // Get the meshfunction evaluations and the map that was stashed in the es.
  std::vector<Real> & final_evals = *es.parameters.get<std::vector<Real> *>("final_evals");
  std::map<dof_id_type, unsigned int> & element_map =
      *es.parameters.get<std::map<dof_id_type, unsigned int> *>("element_map");

  // Setup system vectors and matrices.
  FEType fe_type = system.variable_type(0);
  std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
  QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
  fe->attach_quadrature_rule(&qrule);
  const DofMap & dof_map = system.get_dof_map();
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;
  std::vector<dof_id_type> dof_indices;
  const std::vector<Real> & JxW = fe->get_JxW();
  const std::vector<std::vector<Real>> & phi = fe->get_phi();

  for (const auto & elem : to_mesh.active_local_element_ptr_range())
  {
    fe->reinit(elem);

    dof_map.dof_indices(elem, dof_indices);
    Ke.resize(dof_indices.size(), dof_indices.size());
    Fe.resize(dof_indices.size());

    for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
    {
      Real meshfun_eval = 0.;
      if (element_map.find(elem->id()) != element_map.end())
      {
        // We have evaluations for this element.
        meshfun_eval = final_evals[element_map[elem->id()] + qp];
      }

      // Now compute the element matrix and RHS contributions.
      for (unsigned int i = 0; i < phi.size(); i++)
      {
        // RHS
        Fe(i) += JxW[qp] * (meshfun_eval * phi[i][qp]);

        if (_compute_matrix)
          for (unsigned int j = 0; j < phi.size(); j++)
          {
            // The matrix contribution
            Ke(i, j) += JxW[qp] * (phi[i][qp] * phi[j][qp]);
          }
      }
      dof_map.constrain_element_matrix_and_vector(Ke, Fe, dof_indices);

      if (_compute_matrix)
        system.matrix->add_matrix(Ke, dof_indices);
      system.rhs->add_vector(Fe, dof_indices);
    }
  }
}

void
MultiAppProjectionTransfer::execute()
{
  _console << "Beginning projection transfer " << name() << std::endl;

  getAppInfo();

  // We are going to project the solutions by solving some linear systems.  In
  // order to assemble the systems, we need to evaluate the "from" domain
  // solutions at quadrature points in the "to" domain.  Some parallel
  // communication is necessary because each processor doesn't necessarily have
  // all the "from" information it needs to set its "to" values.  We don't want
  // to use a bunch of big all-to-all broadcasts, so we'll use bounding boxes to
  // figure out which processors have the information we need and only
  // communicate with those processors.
  //
  // Each processor will
  // 1. Check its local quadrature points in the "to" domains to see which
  //    "from" domains they might be in.
  // 2. Send quadrature points to the processors with "from" domains that might
  //    contain those points.
  // 3. Recieve quadrature points from other processors, evaluate its mesh
  //    functions at those points, and send the values back to the proper
  //    processor
  // 4. Recieve mesh function evaluations from all relevant processors and
  //    decide which one to use at every quadrature point (the lowest global app
  //    index always wins)
  // 5. And use the mesh function evaluations to assemble and solve an L2
  //    projection system on its local elements.

  // For every combination of global "from" problem and local "to" problem, find
  // which "from" bounding boxes overlap with which "to" elements.  Keep track
  // of which processors own bounding boxes that overlap with which elements.
  // Build vectors of quadrature points to send to other processors for mesh
  // function evaluations.

  // Get the bounding boxes for the "from" domains.
  std::vector<BoundingBox> bboxes = getFromBoundingBoxes();

  // Figure out how many "from" domains each processor owns.
  std::vector<unsigned int> froms_per_proc = getFromsPerProc();

  std::vector<std::vector<Point>> outgoing_qps(n_processors());
  std::vector<std::map<std::pair<unsigned int, unsigned int>, unsigned int>> element_index_map(
      n_processors());
  // element_index_map[i_to, element_id] = index
  // outgoing_qps[index] is the first quadrature point in element

  if (!_qps_cached)
  {
    for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
    {
      MeshBase & to_mesh = _to_meshes[i_to]->getMesh();

      LinearImplicitSystem & system = *_proj_sys[i_to];

      FEType fe_type = system.variable_type(0);
      std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
      QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
      fe->attach_quadrature_rule(&qrule);
      const std::vector<Point> & xyz = fe->get_xyz();

      unsigned int from0 = 0;
      for (processor_id_type i_proc = 0; i_proc < n_processors();
           from0 += froms_per_proc[i_proc], i_proc++)
      {
        for (const auto & elem :
             as_range(to_mesh.local_elements_begin(), to_mesh.local_elements_end()))
        {
          fe->reinit(elem);

          bool qp_hit = false;
          for (unsigned int i_from = 0; i_from < froms_per_proc[i_proc] && !qp_hit; i_from++)
          {
            for (unsigned int qp = 0; qp < qrule.n_points() && !qp_hit; qp++)
            {
              Point qpt = xyz[qp];
              if (bboxes[from0 + i_from].contains_point(qpt + _to_positions[i_to]))
                qp_hit = true;
            }
          }

          if (qp_hit)
          {
            // The selected processor's bounding box contains at least one
            // quadrature point from this element.  Send all qps from this element
            // and remember where they are in the array using the map.
            std::pair<unsigned int, unsigned int> key(i_to, elem->id());
            element_index_map[i_proc][key] = outgoing_qps[i_proc].size();
            for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
            {
              Point qpt = xyz[qp];
              outgoing_qps[i_proc].push_back(qpt + _to_positions[i_to]);
            }
          }
        }
      }
    }

    if (_fixed_meshes)
      _cached_index_map = element_index_map;
  }
  else
  {
    element_index_map = _cached_index_map;
  }

  // Request quadrature point evaluations from other processors and handle
  // requests sent to this processor.

  // Non-blocking send quadrature points to other processors.
  std::vector<Parallel::Request> send_qps(n_processors());
  if (!_qps_cached)
    for (processor_id_type i_proc = 0; i_proc < n_processors(); i_proc++)
      if (i_proc != processor_id())
        _communicator.send(i_proc, outgoing_qps[i_proc], send_qps[i_proc]);

  // Get the local bounding boxes.
  std::vector<BoundingBox> local_bboxes(froms_per_proc[processor_id()]);
  {
    // Find the index to the first of this processor's local bounding boxes.
    unsigned int local_start = 0;
    for (processor_id_type i_proc = 0; i_proc < n_processors() && i_proc != processor_id();
         i_proc++)
      local_start += froms_per_proc[i_proc];

    // Extract the local bounding boxes.
    for (unsigned int i_from = 0; i_from < froms_per_proc[processor_id()]; i_from++)
      local_bboxes[i_from] = bboxes[local_start + i_from];
  }

  // Setup the local mesh functions.
  std::vector<MeshFunction *> local_meshfuns(froms_per_proc[processor_id()], NULL);
  for (unsigned int i_from = 0; i_from < _from_problems.size(); i_from++)
  {
    FEProblemBase & from_problem = *_from_problems[i_from];
    MooseVariableFEBase & from_var = from_problem.getVariable(
        0, _from_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
    System & from_sys = from_var.sys().system();
    unsigned int from_var_num = from_sys.variable_number(from_var.name());

    MeshFunction * from_func = new MeshFunction(
        from_problem.es(), *from_sys.current_local_solution, from_sys.get_dof_map(), from_var_num);
    from_func->init(Trees::ELEMENTS);
    from_func->enable_out_of_mesh_mode(OutOfMeshValue);
    local_meshfuns[i_from] = from_func;
  }

  // Recieve quadrature points from other processors, evaluate mesh frunctions
  // at those points, and send the values back.
  std::vector<Parallel::Request> send_evals(n_processors());
  std::vector<Parallel::Request> send_ids(n_processors());
  std::vector<std::vector<Real>> outgoing_evals(n_processors());
  std::vector<std::vector<unsigned int>> outgoing_ids(n_processors());
  std::vector<std::vector<Real>> incoming_evals(n_processors());
  std::vector<std::vector<unsigned int>> incoming_app_ids(n_processors());
  for (processor_id_type i_proc = 0; i_proc < n_processors(); i_proc++)
  {
    // Use the cached qps if they're available.
    std::vector<Point> incoming_qps;
    if (!_qps_cached)
    {
      if (i_proc == processor_id())
        incoming_qps = outgoing_qps[i_proc];
      else
        _communicator.receive(i_proc, incoming_qps);
      // Cache these qps for later if _fixed_meshes
      if (_fixed_meshes)
        _cached_qps[i_proc] = incoming_qps;
    }
    else
    {
      incoming_qps = _cached_qps[i_proc];
    }

    outgoing_evals[i_proc].resize(incoming_qps.size(), OutOfMeshValue);
    if (_direction == FROM_MULTIAPP)
      outgoing_ids[i_proc].resize(incoming_qps.size(), libMesh::invalid_uint);
    for (unsigned int qp = 0; qp < incoming_qps.size(); qp++)
    {
      Point qpt = incoming_qps[qp];

      // Loop until we've found the lowest-ranked app that actually contains
      // the quadrature point.
      for (unsigned int i_from = 0; i_from < _from_problems.size(); i_from++)
      {
        if (local_bboxes[i_from].contains_point(qpt))
        {
          outgoing_evals[i_proc][qp] = (*local_meshfuns[i_from])(qpt - _from_positions[i_from]);
          if (_direction == FROM_MULTIAPP)
            outgoing_ids[i_proc][qp] = _local2global_map[i_from];
        }
      }
    }

    if (i_proc == processor_id())
    {
      incoming_evals[i_proc] = outgoing_evals[i_proc];
      if (_direction == FROM_MULTIAPP)
        incoming_app_ids[i_proc] = outgoing_ids[i_proc];
    }
    else
    {
      _communicator.send(i_proc, outgoing_evals[i_proc], send_evals[i_proc]);
      if (_direction == FROM_MULTIAPP)
        _communicator.send(i_proc, outgoing_ids[i_proc], send_ids[i_proc]);
    }
  }

  // Gather all of the qp evaluations and pick out the best ones for each qp.
  for (processor_id_type i_proc = 0; i_proc < n_processors(); i_proc++)
  {
    if (i_proc == processor_id())
      continue;
    _communicator.receive(i_proc, incoming_evals[i_proc]);
    if (_direction == FROM_MULTIAPP)
      _communicator.receive(i_proc, incoming_app_ids[i_proc]);
  }

  std::vector<std::vector<Real>> final_evals(_to_problems.size());
  std::vector<std::map<dof_id_type, unsigned int>> trimmed_element_maps(_to_problems.size());

  for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
  {
    MeshBase & to_mesh = _to_meshes[i_to]->getMesh();
    LinearImplicitSystem & system = *_proj_sys[i_to];

    FEType fe_type = system.variable_type(0);
    std::unique_ptr<FEBase> fe(FEBase::build(to_mesh.mesh_dimension(), fe_type));
    QGauss qrule(to_mesh.mesh_dimension(), fe_type.default_quadrature_order());
    fe->attach_quadrature_rule(&qrule);
    const std::vector<Point> & xyz = fe->get_xyz();

    for (const auto & elem : to_mesh.active_local_element_ptr_range())
    {
      fe->reinit(elem);

      bool element_is_evaled = false;
      std::vector<Real> evals(qrule.n_points(), 0.);

      for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
      {
        Point qpt = xyz[qp];

        unsigned int lowest_app_rank = libMesh::invalid_uint;
        for (unsigned int i_proc = 0; i_proc < n_processors(); i_proc++)
        {
          // Ignore the selected processor if the element wasn't found in it's
          // bounding box.
          std::map<std::pair<unsigned int, unsigned int>, unsigned int> & map =
              element_index_map[i_proc];
          std::pair<unsigned int, unsigned int> key(i_to, elem->id());
          if (map.find(key) == map.end())
            continue;
          unsigned int qp0 = map[key];

          // Ignore the selected processor if it's app has a higher rank than the
          // previously found lowest app rank.
          if (_direction == FROM_MULTIAPP)
            if (incoming_app_ids[i_proc][qp0 + qp] >= lowest_app_rank)
              continue;

          // Ignore the selected processor if the qp was actually outside the
          // processor's subapp's mesh.
          if (incoming_evals[i_proc][qp0 + qp] == OutOfMeshValue)
            continue;

          // This is the best meshfunction evaluation so far, save it.
          element_is_evaled = true;
          evals[qp] = incoming_evals[i_proc][qp0 + qp];
        }
      }

      // If we found good evaluations for any of the qps in this element, save
      // those evaluations for later.
      if (element_is_evaled)
      {
        trimmed_element_maps[i_to][elem->id()] = final_evals[i_to].size();
        for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
          final_evals[i_to].push_back(evals[qp]);
      }
    }
  }

  // We now have just one or zero mesh function values at all of our local
  // quadrature points.  Stash those values (and a map linking them to element
  // ids) in the equation systems parameters and project the solution.

  for (unsigned int i_to = 0; i_to < _to_problems.size(); i_to++)
  {
    _to_es[i_to]->parameters.set<std::vector<Real> *>("final_evals") = &final_evals[i_to];
    _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>("element_map") =
        &trimmed_element_maps[i_to];
    projectSolution(i_to);
    _to_es[i_to]->parameters.set<std::vector<Real> *>("final_evals") = NULL;
    _to_es[i_to]->parameters.set<std::map<dof_id_type, unsigned int> *>("element_map") = NULL;
  }

  for (unsigned int i = 0; i < _from_problems.size(); i++)
    delete local_meshfuns[i];

  // Make sure all our sends succeeded.
  for (processor_id_type i_proc = 0; i_proc < n_processors(); i_proc++)
  {
    if (i_proc == processor_id())
      continue;
    if (!_qps_cached)
      send_qps[i_proc].wait();
    send_evals[i_proc].wait();
    if (_direction == FROM_MULTIAPP)
      send_ids[i_proc].wait();
  }

  if (_fixed_meshes)
    _qps_cached = true;

  _console << "Finished projection transfer " << name() << std::endl;

  postExecute();
}

void
MultiAppProjectionTransfer::projectSolution(unsigned int i_to)
{
  FEProblemBase & to_problem = *_to_problems[i_to];
  EquationSystems & proj_es = to_problem.es();
  LinearImplicitSystem & ls = *_proj_sys[i_to];
  // activate the current transfer
  proj_es.parameters.set<MultiAppProjectionTransfer *>("transfer") = this;

  // TODO: specify solver params in an input file
  // solver tolerance
  Real tol = proj_es.parameters.get<Real>("linear solver tolerance");
  proj_es.parameters.set<Real>("linear solver tolerance") = 1e-10; // set our tolerance
  // solve it
  ls.solve();
  proj_es.parameters.set<Real>("linear solver tolerance") = tol; // restore the original tolerance

  // copy projected solution into target es
  MeshBase & to_mesh = proj_es.get_mesh();

  MooseVariableFEBase & to_var = to_problem.getVariable(
      0, _to_var_name, Moose::VarKindType::VAR_ANY, Moose::VarFieldType::VAR_FIELD_STANDARD);
  System & to_sys = to_var.sys().system();
  NumericVector<Number> * to_solution = to_sys.solution.get();

  for (const auto & node : to_mesh.local_node_ptr_range())
  {
    for (unsigned int comp = 0; comp < node->n_comp(to_sys.number(), to_var.number()); comp++)
    {
      const dof_id_type proj_index = node->dof_number(ls.number(), _proj_var_num, comp);
      const dof_id_type to_index = node->dof_number(to_sys.number(), to_var.number(), comp);
      to_solution->set(to_index, (*ls.solution)(proj_index));
    }
  }
  for (const auto & elem : to_mesh.active_local_element_ptr_range())
    for (unsigned int comp = 0; comp < elem->n_comp(to_sys.number(), to_var.number()); comp++)
    {
      const dof_id_type proj_index = elem->dof_number(ls.number(), _proj_var_num, comp);
      const dof_id_type to_index = elem->dof_number(to_sys.number(), to_var.number(), comp);
      to_solution->set(to_index, (*ls.solution)(proj_index));
    }

  to_solution->close();
  to_sys.update();
}