Helmholtz 2D tutorial with multi source

tutorials/helmholtz/free_space/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/helmholtz/free_space/main.cpp

+#include "gmsh.h"
+
+#include <gmshddm/Formulation.h>
+#include <gmshddm/GmshDdm.h>
+#include <gmshddm/MPIInterface.h>
+#include <gmshfem/Message.h>
+
+using namespace gmshddm;
+using namespace gmshddm::common;
+using namespace gmshddm::domain;
+using namespace gmshddm::problem;
+using namespace gmshddm::field;
+using namespace gmshddm::mpi;
+
+using namespace gmshfem;
+using namespace gmshfem::domain;
+using namespace gmshfem::equation;
+using namespace gmshfem::problem;
+using namespace gmshfem::field;
+using namespace gmshfem::function;
+using namespace gmshfem::post;
+
+#include <cstring>
+#include <fstream>
+#include <sstream>
+
+static std::string freqToStr(double f) {
+  char buffer[100];
+  sprintf(buffer, "%.2f", f);
+  return std::string(buffer);
+}
+
+struct sourcePosition{
+  double x;
+  double y;
+  sourcePosition(double x, double y) : x(x), y(y) {}
+};
+
+std::vector<sourcePosition> readSources(const std::string& path) {
+    std::ifstream inFile(path);
+    std::vector<sourcePosition> positions;
+
+    if (!inFile.is_open()) {
+        throw std::runtime_error("Unable to open file: " + path);
+    }
+
+    std::string line;
+    int lineNumber = 0;
+    while (std::getline(inFile, line)) {
+        lineNumber++;
+        std::istringstream lineStream(line);
+        double x, y;
+        char delimiter;
+
+        if (lineStream >> x >> delimiter >> y && delimiter == ';') {
+            positions.push_back({x, y});
+        } else {
+            throw std::runtime_error("Invalid line format in file: " + path + " at line " + std::to_string(lineNumber));
+        }
+    }
+
+    inFile.close();
+    return positions;
+}
+
+
+/**
+ * Returns ID of the domains containing the source
+*/
+std::vector<int> mesh(double L, double H, double x0, double y0, const std::vector<sourcePosition>& sourcePositions, unsigned nX = 2, unsigned nY = 2, double lc = 0.1) {
+  gmsh::model::add("freeSpace");
+  // Disable warnings for large meshes
+  gmsh::option::setNumber("General.ExpertMode", 1);
+
+  namespace gmodel = gmsh::model;
+
+  for(unsigned i = 0; i < sourcePositions.size(); ++i) {
+      double x_src = sourcePositions[i].x;
+      double y_src = sourcePositions[i].y;
+        if(x_src < x0 || y_src < y0 || x_src > (x0 + L) || y_src > (y0 + H)) {
+            throw gmshfem::common::Exception("Invalid source location.");
+        }
+  }
+
+  std::vector<std::vector<int>> points(nX + 1);
+  std::vector<std::vector<int>> edgesH(nX);
+  std::vector<std::vector<int>> edgesV(nX+1);
+  std::vector<std::vector<int>> curvedLoops(nX);
+  // Indices from 0 to N (N+1 points)
+  for (unsigned i = 0; i <= nX; ++i)
+    points[i].resize(nY + 1);
+
+  const double dx = L / nX;
+  const double dy = H / nY;
+
+  // Create all points
+  for (unsigned i = 0; i <= nX; ++i) {
+    for (unsigned j = 0; j <= nY; ++j) {
+      points[i][j] = gmodel::geo::addPoint(x0 + i * dx, y0 + j * dy, 0.0, lc);
+    }
+  }
+
+  // Create all horizontal edges
+  for (unsigned i = 0; i < nX; ++i) {
+    edgesH[i].resize(nY+1);
+    for (unsigned j = 0; j <= nY; ++j) {
+      edgesH[i][j] = gmodel::geo::addLine(points.at(i).at(j), points.at(i+1).at(j));
+    }
+  }
+
+  // Create all vertical edges
+  for (unsigned i = 0; i <= nX; ++i) {
+    edgesV[i].resize(nY);
+    for (unsigned j = 0; j < nY; ++j) {
+      edgesV[i][j] = gmodel::geo::addLine(points.at(i).at(j), points.at(i).at(j+1));
+    }
+  }
+
+  // Create all surfaces
+  for (unsigned i = 0; i < nX; ++i) {
+    curvedLoops[i].resize(nY);
+    for (unsigned j = 0; j < nY; ++j) {
+      curvedLoops[i][j] = gmodel::geo::addCurveLoop({edgesH.at(i).at(j), edgesV.at(i+1).at(j), -edgesH.at(i).at(j+1), -edgesV.at(i).at(j)});
+      gmodel::geo::addPlaneSurface({curvedLoops[i][j]}, curvedLoops[i][j]); // Same indices
+    }
+  }
+
+  // Create sources. TODO: optimize
+  std::vector<int> sources;
+  std::vector<int> sourcesDomains;
+  for (const auto& pos : sourcePositions) {
+    int source = gmodel::geo::addPoint(pos.x, pos.y, 0, lc);
+    int i_source = int(pos.x / dx);
+    int j_source = int(pos.y / dy);
+    gmodel::geo::synchronize();
+    gmodel::mesh::embed(0, {source}, 2, curvedLoops[i_source][j_source]);
+    sources.push_back(source);
+    sourcesDomains.push_back(i_source + j_source * nX);
+  }
+
+  gmodel::geo::synchronize();
+
+  /* Define physical entities */
+
+  
+
+  // Domains
+  std::vector<int> omegaAll;
+  for (unsigned i = 0; i < nX; ++i) {
+    for (unsigned j = 0; j < nY; ++j) {
+      int omegaIJ = gmodel::addPhysicalGroup(2, {curvedLoops[i][j]}, -1);
+      omegaAll.push_back(omegaIJ);
+      gmodel::setPhysicalName(2, omegaIJ, "omega_" + std::to_string(i) + "_" + std::to_string(j));
+    }
+  }
+
+  int omegaAllIdx = gmodel::addPhysicalGroup(2, omegaAll, -1);
+  gmodel::setPhysicalName(2, omegaAllIdx, "omega");
+
+  // Outer boundary
+  std::vector<int> gammaOut;
+  for (unsigned i = 0; i < nX; ++i)
+  {
+    gammaOut.push_back(edgesH[i][0]);
+    gammaOut.push_back(edgesH[i][nY]);
+  }
+
+  for (unsigned j = 0; j < nY; ++j)
+  {
+    gammaOut.push_back(edgesV[0][j]);
+    gammaOut.push_back(edgesV[nX][j]);
+  }
+
+  int gammaOutIdx = gmodel::addPhysicalGroup(1, gammaOut, -1);
+  gmodel::setPhysicalName(1, gammaOutIdx, "gammaOut");
+
+  // All eges
+  // H-edges
+  for (unsigned i = 0; i < nX; ++i) {
+    for (unsigned j = 0; j <= nY; ++j) {
+      int HIJ = gmodel::addPhysicalGroup(1, {edgesH[i][j]}, -1);
+      gmodel::setPhysicalName(1, HIJ, "edge_h_" + std::to_string(i) + "_" + std::to_string(j));
+    }
+  }
+
+  // V-edges
+  for (unsigned i = 0; i <= nX; ++i) {
+    for (unsigned j = 0; j < nY; ++j) {
+      int VIJ = gmodel::addPhysicalGroup(1, {edgesV[i][j]}, -1);
+      gmodel::setPhysicalName(1, VIJ, "edge_v_" + std::to_string(i) + "_" + std::to_string(j));
+    }
+  }
+
+  // Points
+  for (unsigned i = 0; i <= nX; ++i) {
+    for (unsigned j = 0; j <= nY; ++j) {
+      int VIJ = gmodel::addPhysicalGroup(1, {points[i][j]}, -1);
+      gmodel::setPhysicalName(1, VIJ, "point_" + std::to_string(i) + "_" + std::to_string(j));
+    }
+  }
+
+  // Src
+  for (unsigned srcId = 0; srcId < sources.size(); ++srcId) {
+    int srcEntity = gmodel::addPhysicalGroup(0, {sources[srcId]});
+    gmodel::setPhysicalName(0, srcEntity, "source_" + std::to_string(srcId));
+
+  }
+
+  gmodel::geo::synchronize();
+  gmodel::mesh::generate();
+
+  return sourcesDomains;
+
+}
+
+int main(int argc, char **argv)
+{
+  GmshDdm gmshDdm(argc, argv);
+  constexpr std::complex<double> im(0, 1);
+  constexpr double pi = 3.14159265359;
+
+  unsigned rank = gmshddm::mpi::getMPIRank();
+
+  int order = 2;
+  gmshDdm.userDefinedParameter(order, "order");
+  double f = 2;
+  gmshDdm.userDefinedParameter(f, "f");
+
+  double lc = (1. / f) / 10;
+
+
+  gmshDdm.userDefinedParameter(lc, "lc");
+
+  std::string gauss = "Gauss10";
+  gmshDdm.userDefinedParameter(gauss, "gauss");
+
+  // Output parameters
+  bool mustSolve{true};
+  gmshDdm.userDefinedParameter(mustSolve, "noSolve");
+
+
+  auto nX = 2;
+  auto nY = 2;
+  gmshDdm.userDefinedParameter(nX, "nX");
+  gmshDdm.userDefinedParameter(nY, "nY");
+
+  double L{1}, H{1};
+  gmshDdm.userDefinedParameter(L, "L");
+  gmshDdm.userDefinedParameter(H, "H");
+
+  
+
+  bool save_matrix = false;
+  gmshDdm.userDefinedParameter(save_matrix, "save_iteration_matrix");
+
+  // TODO : Weird behavior for negative y-axis -> Unsupported x0, y0
+  std::vector<sourcePosition> sources;
+  std::string sourcePath = "../srcs.csv";
+  gmshDdm.userDefinedParameter(sourcePath, "sourcePath");
+  sources = readSources(sourcePath);
+  auto d_source = mesh(L, H, 0.0, 0, sources, nX, nY, lc);
+
+
+  // Configure wavenumber. If no model is given, assume homogenuous velocity of 1
+  // (i.e. k = 2 * pi * f)
+  // If a model is given, it is assumed to be a squared slowness.
+  gmshfem::function::Function<std::complex<double>, gmshfem::Degree::Degree0> k;
+
+  std::string model_path = "";
+  std::string model_name;
+  gmshDdm.userDefinedParameter(model_path, "model_path");
+  if (model_path != "") {
+    model_name = model_path.substr(0, model_path.size() - 4);
+    gmsh::merge(model_path);
+    std::vector<int> tags;
+    gmsh::view::getTags(tags);
+
+    auto view = tags.back();
+    ScalarFunction<std::complex<double>> model = probeScalarView<std::complex<double>>(view);
+    // k = omega * slowness
+    k = 2 * pi * f * sqrt(model);
+  }
+  else {
+    // Unit velocity
+    k = 2 * pi * f;
+    model_name = "free_space";
+  }
+  
+  gmshDdm.userDefinedParameter(model_name, "model_name");
+  msg::info << "Model name is " << model_name << gmshfem::msg::endl;
+  
+
+  auto nDom = nX * nY;
+  // Define domain
+  Subdomain omega(nDom);
+  Subdomain gammaOut(nDom);
+  Interface sigma(nDom);
+
+  #define OMEGA(i,j) "omega_" + std::to_string((i)) + "_" + std::to_string((j))
+  #define POINT(i,j) "point_" + std::to_string((i)) + "_" + std::to_string((j))
+  #define EDGE_H(i,j) "edge_h_" + std::to_string((i)) + "_" + std::to_string((j))
+  #define EDGE_V(i,j) "edge_v_" + std::to_string((i)) + "_" + std::to_string((j))
+
+  for(int i = 0; i < nX; ++i) {
+    for (int j = 0; j < nY; ++j) {
+      int d = i + nX * j;
+      omega(d) = Domain(OMEGA(i,j));
+
+      // Left border / interface
+      if(i == 0) {
+        gammaOut(d) |= Domain(EDGE_V(i, j));
+      }
+      else {
+        sigma(d, d - 1) =  Domain(EDGE_V(i, j));
+      }
+      // Right border / interface
+      if(i == nX - 1) {
+        gammaOut(d) |= Domain(EDGE_V(i + 1, j));
+      }
+      else {
+        sigma(d, d + 1) =  Domain(EDGE_V(i + 1, j));
+      }
+      // Bottom border / interface
+      if(j == 0) {
+        gammaOut(d) |= Domain(EDGE_H(i, j));
+      }
+      else {
+        sigma(d, d - nX) =  Domain(EDGE_H(i, j));
+      }
+      // Top border / interface
+      if(j == nY - 1) {
+        gammaOut(d) |= Domain(EDGE_H(i, j + 1));
+      }
+      else {
+        sigma(d, d + nX) =  Domain(EDGE_H(i, j + 1));
+      }
+
+      
+    }
+  }
+
+  // Define topology
+  std::vector< std::vector< unsigned int > > topology(nDom);
+  for(unsigned int i = 0; i < nX; ++i) {
+    for (unsigned int j = 0; j < nY; ++j) {
+        unsigned int d = i + nX * j;
+        if (i != 0) topology[d].push_back(d-1);
+        if (i != nX-1) topology[d].push_back(d+1);
+        if (j != 0) topology[d].push_back(d-nX);
+        if (j != nY-1) topology[d].push_back(d+nX);
+    }
+  }
+
+
+  // Create DDM formulation
+  gmshddm::problem::Formulation< std::complex< double > > formulation("Helmholtz", topology);
+
+  // Create fields
+  SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaOut, FunctionSpaceTypeForm0::HierarchicalH1, order);
+  InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, order);
+  // Tell to the formulation that g is field that have to be exchanged between subdomains
+  formulation.addInterfaceField(g);
+
+  // Add terms common to all RHS
+  for(unsigned int i = 0; i < nDom; ++i) {
+        // VOLUME TERMS
+        formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss);
+        formulation(i).integral(-k * k * dof(u(i)), tf(u(i)), omega(i), gauss);
+        formulation(i).integral(-im * k * dof(u(i)), tf(u(i)), gammaOut(i), gauss);
+
+        // Artificial source
+        for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
+        const unsigned int j = topology[i][jj];
+        formulation(i).integral(-im * k * dof(u(i)), tf(u(i)), sigma(i, j), gauss);
+        formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss);
+        }
+
+        // INTERFACE TERMS
+        for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
+        const unsigned int j = topology[i][jj];
+        formulation(i, j).integral(dof(g(i, j)), tf(g(i, j)), sigma(i, j), gauss);
+        formulation(i, j).integral(formulation.artificialSource(g(j, i)), tf(g(i, j)), sigma(i, j), gauss);
+        formulation(i, j).integral(2. * im * k * u(i), tf(g(i, j)), sigma(i, j), gauss);
+        }
+    }
+
+  // Add all sources
+  formulation.setShotNumber(sources.size());
+  for(unsigned isrc = 0; isrc < sources.size(); ++isrc) {
+    formulation.numberedPhysicalSourceTerm(isrc, formulation(d_source[isrc]).integral(formulation.physicalSource(-1), tf(u(d_source[isrc])), Domain("source_" + std::to_string(isrc)), gauss));
+  }
+
+      formulation.pre(false);
+
+
+  for(unsigned isrc = 0; isrc < sources.size(); ++isrc) {
+    formulation.disableAllShots();
+    formulation.enableShot(isrc);
+
+    // Solve the DDM formulation
+    formulation.assemble();
+
+    if(mustSolve) {
+        int iterMax = 1000;
+        std::string solver = "gmres";
+        gmshDdm.userDefinedParameter(solver, "solver");
+        gmshDdm.userDefinedParameter(iterMax, "iter");
+        formulation.solve(solver, 1e-6, iterMax);
+         save(u.buildUnifiedField(), u.getLocalDomain(), "u_rhs_" + std::to_string(isrc) + "_rank_" + std::to_string(rank));
+    }
+
+
+    std::string suffix = model_name + "_" + freqToStr(f) + "_" + std::to_string(nX) + "_" + std::to_string(nY);
+    std::string mat_name{"matrix_" + suffix}, vec_name{"rhs_" + std::to_string(isrc) + "_" + suffix};
+    if(save_matrix) {
+        // Export all RHS but only one matrix
+        if(isrc == 0) {
+        auto mat = formulation.computeMatrix();
+        msg::info << "nnz = " << mat.numberOfNonZero() << ", symmetric = " << mat.isSymmetric() << ", hermitian = " << mat.isHermitian() << msg::endl;
+        mat.save(mat_name);
+        }
+
+        formulation.getRHS().save(vec_name);
+    }
+  }
+
+
+  return 0;
+}

tutorials/helmholtz/free_space/srcs.csv

+0.4;0.4
+0.4;0.6
+0.6;0.4
+0.6;0.6