tutorial for checking sum of RHS

tutorials/mpi_rhs/CMakeLists.txt

+cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
+
+project(tuto CXX)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../tutorial.cmake)
+
+add_executable(tuto main.cpp ${EXTRA_INCS})
+target_link_libraries(tuto ${EXTRA_LIBS})

tutorials/mpi_rhs/main.cpp

+
+#include <gmsh.h>
+#include <gmshfem/Formulation.h>
+#include <gmshfem/GmshFem.h>
+#include <gmshfem/DistributedField.h>
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::problem;
+using namespace gmshfem::domain;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+using namespace gmshfem::equation;
+
+int main(int argc, char **argv)
+{
+  GmshFem gmshFem(argc, argv);
+
+  double lc = 0.05;
+  gmshFem.userDefinedParameter(lc, "lc");
+  /****
+   * Gmsh part
+   ****/
+
+  gmsh::model::add("poissonDirichlet");
+  gmshfem::function::ScalarFunction<std::complex<double>> sol = x<double>() + 2*y<double>() + 0.5 * pow(x<double>(), 2) + x<double>()*y<double>() + 0.5 * pow(y<double>(), 2);
+
+
+  // small circle
+  gmsh::model::geo::addPoint(0., 0., 0., lc, 1);
+  gmsh::model::geo::addPoint(1, 0., 0., lc, 2);
+  gmsh::model::geo::addPoint(1, 1., 0., lc, 3);
+  gmsh::model::geo::addPoint(0, 1., 0., lc, 4);
+
+  gmsh::model::geo::addLine(1, 2, 1);
+  gmsh::model::geo::addLine(2, 3, 2);
+  gmsh::model::geo::addLine(3, 4, 3);
+  gmsh::model::geo::addLine(4, 1, 4);
+
+  gmsh::model::geo::addCurveLoop({1, 2, 3, 4}, 2);
+
+  // surface
+  gmsh::model::geo::addPlaneSurface({2}, 2);
+  gmsh::model::geo::synchronize();
+
+  // physicals
+  gmsh::model::addPhysicalGroup(2, {2}, 1);
+  gmsh::model::setPhysicalName(2, 1, "omega");
+  gmsh::model::addPhysicalGroup(1, {1, 2, 3, 4}, 1);
+  gmsh::model::setPhysicalName(1, 1, "gammaSmall");
+
+  gmsh::model::geo::synchronize();
+  gmsh::model::mesh::generate();
+  gmsh::model::mesh::createEdges();
+
+  /****
+   * GmshFem part
+   *****/
+
+  if(gmshFem.getMPISize() != 2) {
+    gmshfem::msg::error << "This example needs 2 MPI ranks" << gmshfem::msg::endl;
+    return 1;
+  }
+
+  gmsh::model::mesh::partition(2);
+  gmsh::model::mesh::generateOverlapForEntity(2, 2);
+  std::vector<int> tagsOvlp, tagsOvlpBnd, tagsBndOvlp;
+  auto rank = gmshFem.getMPIRank();
+  gmsh::model::mesh::findOverlapEntities(2, 2, rank + 1, tagsOvlp);
+
+  gmsh::vectorpair entities;
+  gmshfem::domain::Domain omegaOvlp, omegaLoc;
+  gmsh::model::getEntities(entities, 2);
+  for (auto entity: entities) {
+    std::vector<int> parts;
+    gmsh::model::getPartitions(entity.first, entity.second, parts);
+    if (std::find(parts.begin(), parts.end(), rank + 1) != parts.end()) {
+      omegaLoc.addEntity(entity.first, entity.second);
+    } 
+  }
+  
+  omegaOvlp.addEntity(2, tagsOvlp[0]);
+
+  auto computeError = [&omegaOvlp, &omegaLoc, &sol](bool useSum, bool addTermOnOverlap) {
+    Formulation< std::complex< double > > formulation("poissonDirichlet");
+
+    DistributedField< std::complex< double >, Form::Form0 > v("v", omegaLoc, omegaOvlp, FunctionSpaceTypeForm0::HierarchicalH1, 2);
+
+
+    formulation.integral(-dof(v), (tf(v)), omegaLoc | omegaOvlp, "Gauss6");
+    formulation.integral(sol, tf(v), omegaLoc, "Gauss6");
+    if (addTermOnOverlap)
+      formulation.integral(sol, tf(v), omegaOvlp, "Gauss6");
+    formulation.pre();
+    formulation.preDistributed();
+
+    // Pro
+    formulation.assemble();
+    formulation.setSumDistributedRHS(useSum);
+    formulation.solveDistributed(gmshfem::system::DISTRIBUTED_SOLVE_TYPE::DIRECT);
+
+    auto err = sol - v;
+    auto errL2 = gmshfem::post::integrate(err * conj(err), omegaLoc, "Gauss6");
+    auto sumErr = gmshfem::common::GmshFem::currentInstance()->reductionSum(real(errL2));
+    return std::sqrt(sumErr);
+  };
+
+  auto classic = computeError(false, true);
+  auto distributed = computeError(true, false);
+  auto tooMuch = computeError(true, true);
+  auto notEnough = computeError(false, false);
+  if(rank == 0) {
+    gmshfem::msg::info << "Classic: " << classic << gmshfem::msg::endl;
+    gmshfem::msg::info << "Distributed: " << distributed << gmshfem::msg::endl;
+    gmshfem::msg::info << "Sum and term on ovlp: " << tooMuch << gmshfem::msg::endl;
+    gmshfem::msg::info << "No sum and no term on ovlp: " << notEnough << gmshfem::msg::endl;
+  }
+
+
+  return 0;
+}