Select Git revision
split.script
Forked from
gmsh / gmsh
Source project has a limited visibility.
main.cpp 12.26 KiB
#include <gmsh.h>
#include <gmshddm/Formulation.h>
#include <gmshddm/GmshDdm.h>
#include <gmshddm/MPIInterface.h>
#include <gmshfem/Message.h>
#include <gmshfem/AnalyticalFunction.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::analytics;
using namespace gmshfem::post;
void sphere(const double Rin, const double Rext, const double lc)
{
gmsh::model::add("sphere");
int s_in = gmsh::model::occ::addSphere(0, 0, 0, Rin);
int s_ext = gmsh::model::occ::addSphere(0, 0, 0, Rext);
std::vector< std::pair< int, int > > ov_phy;
std::vector< std::vector< std::pair< int, int > > > ovv_phy;
gmsh::model::occ::cut({{3, s_ext}}, {{3, s_in}}, ov_phy, ovv_phy, 5);
gmsh::model::occ::synchronize();
gmsh::model::addPhysicalGroup(2, {s_ext}, 1);
gmsh::model::setPhysicalName(2, 1, "gamma_scat");
gmsh::model::addPhysicalGroup(2, {s_in}, 2);
gmsh::model::setPhysicalName(2, 2, "gamma_ext");
gmsh::model::addPhysicalGroup(3, {5}, 1);
gmsh::model::setPhysicalName(3, 1, "omega");
gmsh::option::setNumber("Mesh.MeshSizeMax", lc);
gmsh::model::occ::synchronize();
gmsh::model::mesh::generate(3);
}
int main(int argc, char **argv)
{
GmshDdm gmshDdm(argc, argv);
const double pi = 3.14159265358979323846264338327950288;
const std::complex< double > im(0., 1.);
std::string mesh = "";
gmshDdm.userDefinedParameter(mesh, "mesh");
int nDom = 2;
gmshDdm.userDefinedParameter(nDom, "nDom");
double lc = 0.1;
gmshDdm.userDefinedParameter(lc, "lc");
double Rin = 0.5;
gmshDdm.userDefinedParameter(Rin, "Rin");
double Rext = 1.5;
gmshDdm.userDefinedParameter(Rin, "Rext");
double k = 1;
gmshDdm.userDefinedParameter(k, "k");
double theta = pi/3.; // orientation of the incident wave - azimuthal angle [0, 2*pi]
gmshDdm.userDefinedParameter(theta, "theta");
double phi = pi/2.; // orientation of the incident wave - polar angle [0, pi]
gmshDdm.userDefinedParameter(phi, "phi");
bool ABC2 = true; // false: Sommerfeld ABC, true: 2nd order ABC
gmshDdm.userDefinedParameter(ABC2, "ABC2");
if(ABC2 && !getMPIRank())
msg::info << " - use 2nd order exterior absorbing boundary condition" << msg::endl;
int FEMorder = 2;
gmshDdm.userDefinedParameter(FEMorder, "FEMorder");
std::string gauss = "Gauss" + std::to_string(2 * FEMorder + 1);
gmsh::option::setNumber("Mesh.ElementOrder", FEMorder);
gmsh::option::setNumber("Mesh.SecondOrderLinear", 1); // straight-sided elements
std::string Transmission = "Taylor2";
gmshDdm.userDefinedParameter(Transmission, "TC");
double alpha = -pi / 2.;
gmshDdm.userDefinedParameter(alpha, "alpha");
// if no mesh is given, create a partitioned mesh of a sphere on rank 0
if(mesh.empty()) {
if(!getMPIRank()) {
float pointsByWl = 2 * pi * FEMorder / k / lc;
msg::info << " - dofs density by wavelength = " << pointsByWl << "" << msg::endl;
if(pointsByWl < 6)
msg::warning << " - less than 6 dofs per wavelength !" << msg::endl;
// build geometry and mesh
sphere(Rin, Rext, lc);
// partition
gmsh::model::mesh::partition(nDom);
// save all, as partitioning will create some entities without physical groups
gmsh::option::setNumber("Mesh.SaveAll", 1);
// save partitioned mesh in single file for sequential runs
gmsh::option::setNumber("Mesh.PartitionSplitMeshFiles", 0);
gmsh::write("sphere.msh");
// save partitioned mesh in separate file for distributed runs
gmsh::option::setNumber("Mesh.PartitionSplitMeshFiles", 1);
gmsh::write("sphere.msh");
gmsh::clear();
}
mesh = "sphere";
}
barrier();
// read partitioned mesh
if(getMPISize() == 1) {
gmsh::merge(mesh + ".msh");
}
else {
for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) {
if (mpi::isItMySubdomain(i)) {
gmsh::merge(mesh + "_" + std::to_string(i + 1) + ".msh");
}
}
}
if(gmsh::model::getNumberOfPartitions() != nDom)
msg::error << "Wrong number of partitions in mesh file " << gmsh::model::getNumberOfPartitions() << " vs. " << nDom << msg::endl;
Domain omegaMono = gmshfem::domain::Domain("omega");
Domain gammaScatMono = gmshfem::domain::Domain("gamma_scat");
Domain gammaExtMono = gmshfem::domain::Domain("gamma_ext");
Subdomain omega = Subdomain::buildSubdomainOf(omegaMono);
Subdomain gammaScat = Subdomain::buildSubdomainOf(gammaScatMono);
Subdomain gammaExt = Subdomain::buildSubdomainOf(gammaExtMono);
std::vector< std::vector< unsigned int > > topology;
Interface sigma = Interface::buildInterface(topology);
gmshddm::problem::Formulation< std::complex< double > > formulation("Helmholtz", topology);
SubdomainField< std::complex< double >, Form::Form0 > u("u", omega | gammaScat | gammaExt | sigma, FunctionSpaceTypeForm0::Lagrange);
InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::Lagrange);
formulation.addInterfaceField(g);
// Define Analytical solution
Function< std::complex< double >, Degree::Degree0 > *u_exact = nullptr;
if(mesh == "sphere") {
u_exact = new AnalyticalFunction< helmholtz3D::ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta, phi); // -k because we use e^iwt convention
u_exact->activateMemorization();
}
// Define incident field
Function< std::complex< double >, Degree::Degree0 > uinc, dn_uinc;
if(mesh == "sphere") {
uinc = exp< std::complex< double > >(im * k * ( sin(phi)*cos(theta)*x< std::complex< double > >() + sin(phi)*sin(theta)*y< std::complex< double > >() + cos(phi)*z< std::complex< double > >() ) );
dn_uinc = AnalyticalFunction< helmholtz3D::dr_ScatteringByAHardSphere< std::complex< double > > >(-k, Rin, 0., 0., 0., 2 * k, theta, phi);
}
else {
uinc = exp< std::complex< double > >(im * k * z< std::complex< double > >());
dn_uinc = normal< std::complex< double > >() * vector< std::complex< double > > (0, 0, im * k * uinc);
}
ScalarFunction< std::complex< double > > Si = 0., dS = 0.;
if(Transmission == "Taylor0") {
if(!getMPIRank())
msg::info << " - Use Sommerfeld transmission condition with rotation angle = " << alpha << "" << msg::endl;
Si = im * k * exp(im * alpha / 2.);
dS = 0.;
}
else if(Transmission == "Taylor2") {
if(!getMPIRank())
msg::info << " - Use second order transmission condition with rotation angle = " << alpha << "" << msg::endl;
Si = im * k * cos(alpha / 2.);
dS = im * exp(-im * alpha / 2.) / (2.0 * k);
}
for(unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) {
// Volume weak form
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);
// Sommerfeld BC
formulation(i).integral(im * k * dof(u(i)), tf(u(i)), gammaExt(i), gauss);
if (ABC2) {
formulation(i).integral( ( 1. / Rext ) * dof(u(i)), tf(u(i)), gammaExt(i), gauss);
formulation(i).integral( ( -im / (2. * k) + 1. / (2. * k * k *Rext) ) * grad(dof(u(i))), grad(tf(u(i))), gammaExt(i), gauss);
}
// Incident field
formulation(i).integral( formulation.physicalSource(dn_uinc) , tf(u(i)), gammaScat(i), gauss);
// Artificial source
for(unsigned int jj = 0; jj < topology[i].size(); ++jj) {
const unsigned int j = topology[i][jj];
if(!sigma(i, j).isEmpty()) {
formulation(i).integral(formulation.artificialSource(-g(j, i)), tf(u(i)), sigma(i, j), gauss);
// Transmission conditions terms
formulation(i).integral(Si * dof(u(i)), tf(u(i)), sigma(i, j), gauss);
formulation(i).integral(-dS * grad(dof(u(i))), grad(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];
if(!sigma(i, j).isEmpty()) {
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);
// Transmission conditions terms
formulation(i, j).integral(-2. * Si * u(i), tf(g(i, j)), sigma(i, j), gauss);
formulation(i, j).integral(2. * dS * grad(u(i)), grad(tf(g(i, j))), sigma(i, j), gauss);
}
}
}
// Solve DDM
double tolerence = 1e-6;
int maxIter = 1000;
formulation.pre();
formulation.solve("gmres", tolerence, maxIter, true);
if(u_exact) {
// Compute analytic L2-error
for(auto i = 0; i < nDom; ++i) {
if(gmshddm::mpi::isItMySubdomain(i)) {
std::complex< double > denLocal = integrate(pow(abs(*u_exact), 2), omega(i), gauss);
std::complex< double > numLocal = integrate(pow(abs(*u_exact - u(i)), 2), omega(i), gauss);
double num = numLocal.real();
double den = denLocal.real();
double local_err = sqrt(num / den);
msg::info << "Analytic-ddm L2 error (%) = " << 100.*local_err << " on subdomain " << i << msg::endl;
}
}
}
// save solutions
gmsh::option::setNumber("Mesh.Binary", 1);
gmsh::option::setNumber("PostProcessing.Binary", 1);
gmsh::option::setNumber("PostProcessing.SaveMesh", 0);
for(auto i = 0; i < nDom; ++i) {
if(gmshddm::mpi::isItMySubdomain(i)) {
int tag = save(u(i), omega(i), "u", "", "", true, 0, 0., i + 1);
gmsh::view::write(tag, "u_" + std::to_string(i + 1) + ".msh");
// also save a cut
for(int step = 0; step < 2; step++) {
gmsh::view::option::setNumber(tag, "TimeStep", step);
gmsh::view::option::setNumber(tag, "AdaptVisualizationGrid", 1);
gmsh::plugin::setNumber("CutPlane", "View", gmsh::view::getIndex(tag));
gmsh::plugin::setNumber("CutPlane", "A", 0);
gmsh::plugin::setNumber("CutPlane", "B", 0);
gmsh::plugin::setNumber("CutPlane", "C", 1);
gmsh::plugin::setNumber("CutPlane", "D", 0);
gmsh::plugin::setNumber("CutPlane", "RecurLevel", FEMorder);
gmsh::plugin::setNumber("CutPlane", "TargetError", -1); // 1e-4 to get decent AMR
int tag2 = gmsh::plugin::run("CutPlane");
std::string name2 = std::string("u_cut_step") + std::to_string(step) + "_" + std::to_string(i + 1);
gmsh::view::option::setString(tag2, "Name", name2);
gmsh::view::write(tag2, name2 + ".pos");
}
}
}
bool ComputeMono = false; // monodomain solution
gmshDdm.userDefinedParameter(ComputeMono, "ComputeMono");
if (ComputeMono) {
// monodomain solution
gmshfem::field::Field< std::complex< double >, Form::Form0 > uMono("uMono", omegaMono | gammaScatMono | gammaExtMono, FunctionSpaceTypeForm0::Lagrange);
gmshfem::problem::Formulation< std::complex< double > > monodomain("Helmholtz");
monodomain.integral(grad(dof(uMono)), grad(tf(uMono)), omegaMono, gauss);
monodomain.integral(-k * k * dof(uMono), tf(uMono), omegaMono, gauss);
// Sommerfeld BC
monodomain.integral(im * k * dof(uMono), tf(uMono), gammaExtMono, gauss);
if (ABC2) { // 2nd order ABC
monodomain.integral( ( 1. / Rext ) * dof(uMono), tf(uMono), gammaExtMono, gauss);
monodomain.integral( ( -im / (2. * k) + 1. / (2. * k * k *Rext) ) * grad(dof(uMono)), grad(tf(uMono)), gammaExtMono, gauss);
}
// Incident field
monodomain.integral(dn_uinc, tf(uMono), gammaScatMono, gauss);
// solve
monodomain.pre();
monodomain.assemble();
monodomain.solve();
// Compute analytic L2-error
double numM = 0., denM = 0.;
double mono_err; // L2 local error
for (unsigned int i = 0; i < static_cast< unsigned int >(nDom); ++i) {
if(gmshddm::mpi::isItMySubdomain(i)) {
std::complex< double > denMono = integrate(pow(abs(uMono), 2), omega(i), gauss);
std::complex< double > numMono = integrate(pow(abs(uMono - u(i)), 2), omega(i), gauss);
numM += numMono.real();
denM += denMono.real();
mono_err = sqrt(numM / denM);
msg::info << "Monodomain-ddm L2 error (%) = " << 100.*mono_err << " on subdomain " << i << msg::endl;
}
}
} // end Monodomain
return 0;
}