diff --git a/contrib/bessel/CMakeLists.txt b/contrib/bessel/CMakeLists.txt index 357acc318fa53a64c11e9532c4f924c9e8f19212..898de89d96f406cd063f9bb217ed9f4795aa1281 100644 --- a/contrib/bessel/CMakeLists.txt +++ b/contrib/bessel/CMakeLists.txt @@ -12,4 +12,4 @@ if(ENABLE_BESSEL) endif(ENABLE_BESSEL) file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h) -append_gmshfem_contrib(Numeric "${SRC};${HDR}") +append_gmshfem_contrib(bessel "${SRC};${HDR}") diff --git a/examples/helmholtz2d/CMakeLists.txt b/examples/helmholtz2d/diskScattering/CMakeLists.txt similarity index 100% rename from examples/helmholtz2d/CMakeLists.txt rename to examples/helmholtz2d/diskScattering/CMakeLists.txt diff --git a/examples/helmholtz2d/README.md b/examples/helmholtz2d/diskScattering/README.md similarity index 100% rename from examples/helmholtz2d/README.md rename to examples/helmholtz2d/diskScattering/README.md diff --git a/examples/helmholtz2d/main.cpp b/examples/helmholtz2d/diskScattering/main.cpp old mode 100755 new mode 100644 similarity index 100% rename from examples/helmholtz2d/main.cpp rename to examples/helmholtz2d/diskScattering/main.cpp diff --git a/examples/helmholtz2d/waveguide/Readme.md b/examples/helmholtz2d/waveguide/Readme.md new file mode 100644 index 0000000000000000000000000000000000000000..67cc7e981ca30d85e3cbf4ffbb31e9ee9ddf7b19 --- /dev/null +++ b/examples/helmholtz2d/waveguide/Readme.md @@ -0,0 +1,99 @@ +# Assessment of Local Absorbing Boundary Conditions for Heterogeneous Helmholtz Problems +> Marchner P. + +Dependency(ies): +* GmshFEM : v 1.0.0 +* Gmsh : v 4.8.2 + +## Problem description + +This code solves a 2d straight waveguide Helmholtz problem +$$\rho_0^{-1}\partial_x(\rho_0\partial_x) u + \rho_0^{-1}\partial_y(\rho_0\partial_y) u + \omega^2 c_0^{-2} u = 0, \quad \text{in} \; \Omega=[0, L]\times[0, H]$$ +with a either +* a longitudinal, linear speed of sound profile $$c_0(x)=ax+b, \, a>0$$ (`longitudinal/`), +* a transverse, Gaussian $$c_0(y) ,\, \rho_0(y)$$ variation of the speed of sound and density (`transverse/`). + +High-order finite elements are used to discretize the weak formulations. +A given mode of the form $$g=\cos(n \pi y / H), n\in\mathbb{N}$$ is imposed as input boundary condition at $$x=0$$. +Homogeneous Dirichlet or Neumann boundary conditions are available for the upper and lower walls. + +Different ABCs are tested on the output boundary. +* In the folder `longitudinal/`, the relative $$L^2$$-error is computed thanks to an analytic solution. +* In the folder `transverse/`, the relative $$L^2$$-error is computed thanks to a PML reference solution. +For the PML reference, Bermudez' function is used with the parameter $$\sigma_0=40$$ and 40 layers in order to ensure the attenuation of all modes including grazing modes. + +## Installation + +For the examples in the folder `longitudinal`, GmshFEM must be compiled with the Complex Bessel library + +``` +cmake -DENABLE_BESSEL=1 .. +``` + +Once GmshFEM is installed, go to the folder and run + +``` + mkdir build && cd build + cmake .. + make +``` + +## Usage + +Simply run: +``` + ./example [PARAM] +``` +with `[PARAM]`: +* `-h [x]` where `[x]` is the mesh size. +* `-L [x]` is the waveguide length and ABC position. +* `-H [x]` is the waveguide height. +* `-FEMorder [x]` is the order of the finite element hierarchical basis. +* `-problem [x]` could be `hard` or `soft` and corresponds respectively to homogeneous Neumann and Dirichlet boundary conditions for the upper and lower duct walls. +* `-omega [x]` is the angular frequency. +* `-mode [x]` is the mode number for the input boundary condition. +* `-abcName [x]` is the ABC imposed on the right boundary. +* `-padeOrder [x]` is the Pade order if `abcName = Pade`. +* `-angle [x]` is the rotation angle of the branch-cut if `abcName = Pade`. +* `-runcase [x]` selects a specific configuration for the results reproductibility. + +Each subfolder contains +* the main script (`main.cpp`) +* a loop for testing the ABCs over a given frequency range (`main_omega.cpp`) +* a loop for varying the position of the ABC (`main_length.cpp`) +To run one of the loops just change the name in the `CMakeLists.txt` and recompile + +## Results reproductibility +##### Folder `longitudinal/` reproduce the results from Section 3 +Pressure maps from Figure 2 can be obtained by compiling `main.cpp` and running for example +``` + ./example -omega 20 +``` +The solution `v_exact.msh` can be open in Gmsh. + +For Figure 3 use `main_omega.cpp` and run +``` + ./example -runcase [i] +``` +with i=1 for Taylor-based ABCs or i=2 for Padé-based ABCs. + +For Figure 4 use `main_length.cpp` +``` + ./example -runcase [i] +``` +with i={1,2,3} to obtain the data from Figures 4a, 4b and 4c respectively. + +##### Folder `transverse/` reproduce the results from Section 5 +For Figures 13a, 13b and 13c, use respectively `main.cpp`, `main_omega.cpp` and `main_length.cpp` with +``` +./example -runcase 1 +``` + +Figure 15a and 15b, use respectively `main.cpp` and `main_omega.cpp` with +``` +./example -runcase 2 +``` + +## Reference + +// Preprint diff --git a/examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt b/examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6f8caaab62b9209c4cd7f98d220eaa43d6841dbe --- /dev/null +++ b/examples/helmholtz2d/waveguide/longitudinal/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(example CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake) + +add_executable(example main.cpp ${EXTRA_INCS}) +target_link_libraries(example ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/main.cpp b/examples/helmholtz2d/waveguide/longitudinal/main.cpp similarity index 99% rename from examples/helmholtzFlow/helmholtzDuctHeteroX/main.cpp rename to examples/helmholtz2d/waveguide/longitudinal/main.cpp index fda2673ce34a6eb1afbc715c2206fe8ec13f081a..0661c47df21936eaa36a6930efb3b4e2c1113755 100644 --- a/examples/helmholtzFlow/helmholtzDuctHeteroX/main.cpp +++ b/examples/helmholtz2d/waveguide/longitudinal/main.cpp @@ -76,7 +76,6 @@ int main(int argc, char **argv) int FEMorder = 6; gmshFem.userDefinedParameter(FEMorder, "FEMorder"); std::string gauss = "Gauss"+std::to_string(2*FEMorder+1); - gmshFem.userDefinedParameter(gauss, "gauss"); // meshing double L = 1; // duct length diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp b/examples/helmholtz2d/waveguide/longitudinal/main_length.cpp similarity index 99% rename from examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp rename to examples/helmholtz2d/waveguide/longitudinal/main_length.cpp index c953106e8f04dbda837e98d3cb20b99072f024dd..2bce8fb842a742516d9893fc834f706e5cb63ba8 100644 --- a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_length.cpp +++ b/examples/helmholtz2d/waveguide/longitudinal/main_length.cpp @@ -80,7 +80,6 @@ int main(int argc, char **argv) int FEMorder = 6; gmshFem.userDefinedParameter(FEMorder, "FEMorder"); std::string gauss = "Gauss"+std::to_string(2*FEMorder+1); - gmshFem.userDefinedParameter(gauss, "gauss"); // meshing double H = 0.5; // duct height diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp b/examples/helmholtz2d/waveguide/longitudinal/main_omega.cpp similarity index 99% rename from examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp rename to examples/helmholtz2d/waveguide/longitudinal/main_omega.cpp index a2f302ac9e6d2d4565589a990527f842bfcbda51..73421a4362a6fbaa118003dc057b6528bb500e5d 100644 --- a/examples/helmholtzFlow/helmholtzDuctHeteroX/main_loop_omega.cpp +++ b/examples/helmholtz2d/waveguide/longitudinal/main_omega.cpp @@ -79,7 +79,6 @@ int main(int argc, char **argv) int FEMorder = 6; gmshFem.userDefinedParameter(FEMorder, "FEMorder"); std::string gauss = "Gauss"+std::to_string(2*FEMorder+1); - gmshFem.userDefinedParameter(gauss, "gauss"); // meshing double L = 1; // duct length diff --git a/examples/helmholtz2d/waveguide/transverse/CMakeLists.txt b/examples/helmholtz2d/waveguide/transverse/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6f8caaab62b9209c4cd7f98d220eaa43d6841dbe --- /dev/null +++ b/examples/helmholtz2d/waveguide/transverse/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(example CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake) + +add_executable(example main.cpp ${EXTRA_INCS}) +target_link_libraries(example ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/main.cpp b/examples/helmholtz2d/waveguide/transverse/main.cpp similarity index 100% rename from examples/helmholtzFlow/helmholtzDuctHeteroY/main.cpp rename to examples/helmholtz2d/waveguide/transverse/main.cpp diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_length.cpp b/examples/helmholtz2d/waveguide/transverse/main_length.cpp similarity index 100% rename from examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_length.cpp rename to examples/helmholtz2d/waveguide/transverse/main_length.cpp diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_omega.cpp b/examples/helmholtz2d/waveguide/transverse/main_omega.cpp similarity index 100% rename from examples/helmholtzFlow/helmholtzDuctHeteroY/main_loop_omega.cpp rename to examples/helmholtz2d/waveguide/transverse/main_omega.cpp diff --git a/examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt b/examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6f8caaab62b9209c4cd7f98d220eaa43d6841dbe --- /dev/null +++ b/examples/helmholtzFlow/flowPastCylinder/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(example CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake) + +add_executable(example main.cpp ${EXTRA_INCS}) +target_link_libraries(example ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/main.cpp b/examples/helmholtzFlow/flowPastCylinder/main.cpp similarity index 100% rename from examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/main.cpp rename to examples/helmholtzFlow/flowPastCylinder/main.cpp diff --git a/examples/helmholtzFlow/freefield/CMakeLists.txt b/examples/helmholtzFlow/freefield/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..93e63bab9656fc6f30083e44df4fbb07266fe3aa --- /dev/null +++ b/examples/helmholtzFlow/freefield/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(example CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../example.cmake) + +add_executable(example main.cpp ${EXTRA_INCS}) +target_link_libraries(example ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/freefield/main.cpp b/examples/helmholtzFlow/freefield/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a63b33ac7d6fb961aceeeca426a2aaa10021a250 --- /dev/null +++ b/examples/helmholtzFlow/freefield/main.cpp @@ -0,0 +1,564 @@ +#include "GmshFem.h" +#include "Formulation.h" +#include "gmsh.h" +#include "AnalyticalFunction.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::analytics; +using namespace gmshfem::post; +using namespace gmshfem::equation; + +struct Edge +{ + Domain gamma; + Domain corner[2]; +}; + +struct HelmholtzDomain +{ + Domain omega; + Domain gammaPhy; + Domain gammaPml; + Domain omegaPml; + Domain source; + Domain corners; + std::vector< Domain > corner; + std::vector< Edge > edge; +}; + +//*************************************** +// GEOMETRY +//*************************************** + +void meshCircle(GmshFem &gmshFem, HelmholtzDomain &domains, const double R, const double Rpml, const double xs, const double ys, const double lc, bool withPml) +{ + gmsh::model::add("geometry"); + + gmsh::model::geo::addPoint(0., 0., 0., lc, 1); // Center + gmsh::model::geo::addPoint(xs, ys, 0., lc, 100); // Source point + // physical circle + gmsh::model::geo::addPoint(R, 0., 0., lc, 2); + gmsh::model::geo::addPoint(0, R, 0., lc, 3); + gmsh::model::geo::addPoint(-R, 0., 0., lc, 4); + gmsh::model::geo::addPoint(0, -R, 0., lc, 5); + + gmsh::model::geo::addCircleArc(2, 1, 3, 1); + gmsh::model::geo::addCircleArc(3, 1, 4, 2); + gmsh::model::geo::addCircleArc(4, 1, 5, 3); + gmsh::model::geo::addCircleArc(5, 1, 2, 4); + gmsh::model::geo::addCurveLoop({1,2,3,4}, 1); + // surface + gmsh::model::geo::addPlaneSurface({1}, 1); + + if(withPml) { + gmsh::model::geo::addPoint(Rpml, 0., 0., lc, 6); + gmsh::model::geo::addPoint(0, Rpml, 0., lc, 7); + gmsh::model::geo::addPoint(-Rpml, 0., 0., lc, 8); + gmsh::model::geo::addPoint(0, -Rpml, 0., lc, 9); + + gmsh::model::geo::addCircleArc(6, 1, 7, 5); + gmsh::model::geo::addCircleArc(7, 1, 8, 6); + gmsh::model::geo::addCircleArc(8, 1, 9, 7); + gmsh::model::geo::addCircleArc(9, 1, 6, 8); + gmsh::model::geo::addCurveLoop({5,6,7,8}, 2); + gmsh::model::geo::addPlaneSurface({1,2}, 2); + } + gmsh::model::geo::synchronize(); + gmsh::model::mesh::embed(0,{100},2,1); + + // physicals + gmsh::model::addPhysicalGroup(2, {1}, 1); + gmsh::model::setPhysicalName(2, 1, "omega"); + gmsh::model::addPhysicalGroup(1, {1,2,3,4}, 3); + gmsh::model::setPhysicalName(1, 3, "gammaPhy"); + gmsh::model::addPhysicalGroup(0, {100}, 100); + gmsh::model::setPhysicalName(0, 100, "source"); + if(withPml) { + gmsh::model::addPhysicalGroup(2, {2}, 2); + gmsh::model::setPhysicalName(2, 2, "omegaPml"); + gmsh::model::addPhysicalGroup(1, {5,6,7,8}, 4); + gmsh::model::setPhysicalName(1, 4, "gammaPml"); + } + + // generate mesh + gmsh::model::geo::synchronize(); + gmsh::model::mesh::generate(2); + // gmsh::model::mesh::setOrder(2); // T6 elem + // gmsh::write("m.msh"); + + // Domains + domains.omega = Domain(2, 1); + domains.gammaPhy = Domain(1, 3); + domains.source = Domain(0, 100); + + if(withPml) { + domains.omegaPml = Domain(2, 2); + domains.gammaPml = Domain(1, 4); + } +} + +void meshSquare(GmshFem &gmshFem, HelmholtzDomain &domains, const double L, const double Lpml, const double xs, const double ys, const double lc, bool withPml) +{ + gmsh::model::add("geometry"); + + gmsh::model::geo::addPoint(0., 0., 0., lc, 1); // Center + gmsh::model::geo::addPoint(xs, ys, 0., lc, 100); // Source point + + // polygone + gmsh::model::geo::addPoint(L, L, 0., lc, 2); + gmsh::model::geo::addPoint(-L, L, 0., lc, 3); + gmsh::model::geo::addPoint(-L, -L, 0., lc, 4); + gmsh::model::geo::addPoint(L, -L, 0., lc, 5); + + gmsh::model::geo::addLine(2, 3, 1); + gmsh::model::geo::addLine(3, 4, 2); + gmsh::model::geo::addLine(4, 5, 3); + gmsh::model::geo::addLine(5, 2, 4); + + gmsh::model::geo::addCurveLoop({1,2,3,4}, 1); + gmsh::model::geo::addPlaneSurface({1}, 1); + + if(withPml) { + gmsh::model::geo::addPoint(Lpml, Lpml, 0., lc, 6); + gmsh::model::geo::addPoint(-Lpml, Lpml, 0., lc, 7); + gmsh::model::geo::addPoint(-Lpml, -Lpml, 0., lc, 8); + gmsh::model::geo::addPoint(Lpml, -Lpml, 0., lc, 9); + + gmsh::model::geo::addLine(6, 7, 5); + gmsh::model::geo::addLine(7, 8, 6); + gmsh::model::geo::addLine(8, 9, 7); + gmsh::model::geo::addLine(9, 6, 8); + + gmsh::model::geo::addCurveLoop({5,6,7,8}, 2); + gmsh::model::geo::addPlaneSurface({1,2}, 2); + } + gmsh::model::geo::synchronize(); + gmsh::model::mesh::embed(0,{100},2,1); + // physicals + gmsh::model::addPhysicalGroup(2, {1}, 1); + gmsh::model::setPhysicalName(2, 1, "omega"); + gmsh::model::addPhysicalGroup(1, {1,2,3,4}, 3); + gmsh::model::setPhysicalName(1, 3, "gammaPhy"); + gmsh::model::addPhysicalGroup(0, {100}, 100); + gmsh::model::setPhysicalName(0, 100, "source"); + // Physical edges + gmsh::model::addPhysicalGroup(1, {1}, 20); + gmsh::model::setPhysicalName(1, 20, "gammaPhy_top"); + gmsh::model::addPhysicalGroup(1, {2}, 30); + gmsh::model::setPhysicalName(1, 30, "gammaPhy_left"); + gmsh::model::addPhysicalGroup(1, {3}, 40); + gmsh::model::setPhysicalName(1, 40, "gammaPhy_bottom"); + gmsh::model::addPhysicalGroup(1, {4}, 50); + gmsh::model::setPhysicalName(1, 50, "gammaPhy_right"); + // Physical domain corners + gmsh::model::addPhysicalGroup(0, {2}, 2); + gmsh::model::setPhysicalName(0, 2, "UpRightCorner"); + gmsh::model::addPhysicalGroup(0, {3}, 3); + gmsh::model::setPhysicalName(0, 3, "UpLeftCorner"); + gmsh::model::addPhysicalGroup(0, {4}, 4); + gmsh::model::setPhysicalName(0, 4, "BottomLeftCorner"); + gmsh::model::addPhysicalGroup(0, {5}, 5); + gmsh::model::setPhysicalName(0, 5, "BottomRightCorner"); + if(withPml) { + gmsh::model::addPhysicalGroup(2, {2}, 2); + gmsh::model::setPhysicalName(2, 2, "omegaPml"); + gmsh::model::addPhysicalGroup(1, {5,6,7,8}, 4); + gmsh::model::setPhysicalName(1, 4, "gammaPml"); + } + + // generate mesh + gmsh::model::geo::synchronize(); + gmsh::model::mesh::generate(2); + //gmsh::model::mesh::setOrder(2); // T6 elem + // gmsh::write("m.msh"); + + // Domains + domains.omega = Domain(2, 1); + domains.gammaPhy = Domain(1, 3); + domains.source = Domain(0, 100); + + if(withPml) { + domains.omegaPml = Domain(2, 2); + domains.gammaPml = Domain(1, 4); + } + domains.edge.resize(4); + domains.corner.resize(4); + for(int i = 0; i < 4; ++i) { + domains.edge[i].gamma = Domain(1, 10*(i + 2) ); + domains.edge[i].corner[0] = Domain(0, i + 2 ); + domains.edge[i].corner[1] = Domain(0, (i + 1) % 4 + 2 ); + + domains.corner[i] = Domain(0, i + 2); + + domains.corners |= domains.edge[i].corner[0]; + domains.gammaPhy |= domains.edge[i].gamma; + } +} + +//*************************************** +// FORMULATION +//*************************************** + +int main(int argc, char **argv) +{ + GmshFem gmshFem(argc, argv); + + // geometry parameters + int geometry = 0; // 0-circle, 1-square + gmshFem.userDefinedParameter(geometry, "geometry"); + double R = 1.0; // Circle radius or square of size [-R,R] x [-R,R] + gmshFem.userDefinedParameter(R, "R"); + // Point source location + double xs = 0.3; + double ys = 0.2; + gmshFem.userDefinedParameter(xs, "xs"); + gmshFem.userDefinedParameter(ys, "ys"); + msg::info << "Point source located at xs=" << xs << ", ys=" << ys << msg::endl; + + // physical parameters + const double pi = 3.14159265358979323846264338327950288; + const double k = 3 * pi; // free field wavenumber + const std::complex< double > im(0., 1.); + double M = 0.; // Mach number + gmshFem.userDefinedParameter(M, "M"); + double theta = pi/4; // mean flow orientation + gmshFem.userDefinedParameter(theta, "theta"); + msg::info << "Mach number " << M << " oriented at theta=" << theta << " rad" << msg::endl; + + // numerical parameters + int order = 4; // FEM shape function order + gmshFem.userDefinedParameter(order, "order"); + std::string gauss = "Gauss10"; // increase Nr of Gauss points when order > 6 + double lc = 0.03; // mesh size + gmshFem.userDefinedParameter(lc, "lc"); + + // Choose ABC - exterior boundary condition + std::string abcName = "Pade"; // available choices: ABC-0, ABC-2, Pade, pml + gmshFem.userDefinedParameter(abcName, "abcName"); + int Npml; + double Rpml; + if (abcName == "pml") { // PML parameters + Npml = 4; // number of layers + gmshFem.userDefinedParameter(Npml, "Npml"); + Rpml = R+Npml*lc; // extended domain + } + + HelmholtzDomain domains; + + switch (geometry) { + case 0: + meshCircle(gmshFem, domains, R, Rpml, xs, ys, lc, abcName == "pml"); + break; + case 1: + meshSquare(gmshFem, domains, R, Rpml, xs, ys, lc, abcName == "pml"); + break; + default: + msg::error << "geometry not available ! " << msg::endl; + exit(0); + break; + } + + float pointsByWl = (2*pi*order*(1-abs(M)))/(lc*k); + msg::info << " - FEM basis order = " << order << "" << msg::endl; + msg::info << " - Smallest number of dofs by wavelength = " << pointsByWl << "" << msg::endl; + if (pointsByWl < 6) { + msg::warning << " - Less than 6 points per wavelength ! " << msg::endl; + } + + // Mach number projections + double Mx = M*std::cos(theta); + double My = M*std::sin(theta); + // Mach-velocity vector + VectorFunction< std::complex< double > > MM = vector< std::complex< double > > (Mx,My,0.); + // Normal and tangential projections + ScalarFunction< std::complex< double > > Mn = vector< std::complex< double > > (Mx,My,0.)*normal< std::complex< double > >(); + ScalarFunction< std::complex< double > > Mt = vector< std::complex< double > > (Mx,My,0.)*tangent< std::complex< double > >(); + // Parameters for Inverse Lorentz transformation + double beta = sqrt(1-M*M); // Jacobian + double Alphax = 1 + Mx*Mx/(beta*(1+beta)); + double Alphay = 1 + My*My/(beta*(1+beta)); + double K = Mx*My/(beta*(1+beta)); + // Tensor of the Inverse Lorentz transformation + TensorFunction< std::complex< double > > Linv = tensor< std::complex <double > > (beta*Alphay, -beta*K, 0., -beta*K, beta*Alphax,0.,0.,0.,0.); + + std::vector< FieldInterface< std::complex< double > > * > fieldBucket; + Formulation< std::complex< double > > formulation("helmholtzflow"); + Field< std::complex< double >, Form::Form0 > u("u", domains.omega | domains.omegaPml | domains.gammaPhy | domains.gammaPml | domains.corners | domains.source, FunctionSpaceTypeForm0::HierarchicalH1, order); + + // Define analytic solution + Function< std::complex< double >, Degree::Degree0 > *solution = nullptr; + solution = new AnalyticalFunction< helmholtz2D::MonopoleFreeField< std::complex< double > > >(k, M, theta, 1., xs, ys, 0., 0.); + solution->activateMemorization(); + + // convected Helmholz weak form + formulation.galerkin(vector< std::complex< double > >(1-Mx*Mx, -Mx*My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), domains.omega, gauss); + formulation.galerkin(vector< std::complex< double > >(-Mx*My, 1-My*My, 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omega, gauss); + formulation.galerkin(- k * k * dof(u), tf(u), domains.omega, gauss); + formulation.galerkin(dof(u), vector< std::complex< double > >(-im*k*Mx,-im*k*My, 0.) * grad(tf(u)), domains.omega, gauss); + formulation.galerkin(vector< std::complex< double > >(im*k*Mx, im*k*My, 0.) * grad(dof(u)), tf(u), domains.omega, gauss); + + // Source point + formulation.galerkin(-1.,tf(u),domains.source,gauss); + + // Exact DtN: beta_n**2 * kx = - Mn * (k0 - Mt*ky) + sqrt( k0**2 -2*k0*Mt*ky - (1-abs(M))^2 ky**2 ) + if(abcName == "ABC-0") { // zeroth order Taylor approximation + // boundary contributions + formulation.galerkin(im*k*Mn * dof(u), tf(u), domains.gammaPhy, gauss); + formulation.galerkin(Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss); + // ABC + msg::info << "Use zeroth order ABC" << msg::endl; + formulation.galerkin(im*k*(1-Mn) * dof(u), tf(u), domains.gammaPhy, gauss); + formulation.galerkin(-Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss); + } + else if (abcName == "ABC-2") { + // boundary contributions + formulation.galerkin(im*k*Mn * dof(u), tf(u), domains.gammaPhy, gauss); + formulation.galerkin(Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss); + // ABC + msg::info << "Use second order ABC" << msg::endl; + // zeroth order contribution + formulation.galerkin(im*k*(1-Mn) * dof(u), tf(u), domains.gammaPhy, gauss); + // first order contribution + formulation.galerkin( (1-Mn) * Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss); + // second order contribution + formulation.galerkin(- im * (beta*beta) /(2*k) * grad(dof(u)), grad(tf(u)), domains.gammaPhy, gauss); + } + else if (abcName == "Pade") { + int padeOrder = 2; + gmshFem.userDefinedParameter(padeOrder, "padeOrder"); + double angle = -pi/4.; + gmshFem.userDefinedParameter(angle, "angle"); + msg::info << "Use Pade ABC of order " << padeOrder << " with angle " << angle << " rad"<< msg::endl; + + const double Np = 2. * padeOrder + 1.; + const std::complex< double > exp1 = std::complex<double>(std::cos(angle),std::sin(angle)); + const std::complex< double > exp2 = std::complex<double>(std::cos(angle/2.),std::sin(angle/2.)); + const std::complex< double > coef = 2./Np; + std::vector< std::complex< double > > c(padeOrder, 0.); + for(int i = 0; i < padeOrder; ++i) { + c[i] = std::tan((i + 1) * pi / Np); + c[i] *= c[i]; + } + + if(geometry == 0) { + // define the auxiliary fields + std::vector< Field< std::complex< double >, Form::Form0 >* > phi; + for(int i = 0; i < padeOrder; ++i) { + phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i), domains.gammaPhy, FunctionSpaceTypeForm0::HierarchicalH1, order)); + fieldBucket.push_back(phi.back()); + } + + // write the augmented weak form - approximation of the square-root + formulation.galerkin( im * k * exp2 * dof(u), tf(u), domains.gammaPhy, gauss); + for(int i = 0; i < padeOrder; ++i) { + // boundary integral terms relating the auxiliary fields + formulation.galerkin( im * k * exp2 * coef * c[i] * dof(*phi[i]), tf(u), domains.gammaPhy, gauss); + formulation.galerkin( im * k * exp2 * coef * c[i] * dof(u), tf(u), domains.gammaPhy, gauss); + + // coupling of the auxiliary equations + formulation.galerkin(- (beta*beta) * grad(dof(*phi[i])), grad(tf(*phi[i])), domains.gammaPhy, gauss); + formulation.galerkin(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[i])), tf(*phi[i]), domains.gammaPhy, gauss); + formulation.galerkin( (k*k) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), domains.gammaPhy, gauss); + formulation.galerkin( (k*k) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[i]), domains.gammaPhy, gauss); + } + } + else if(geometry == 1) { + // define the auxilary fields + std::vector< Field< std::complex< double >, Form::Form0 > * > phi[4]; + for(int x = 0; x < 4; ++x) { + for(int i = 0; i < padeOrder; ++i) { + phi[x].push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i) + "^" + std::to_string(x), domains.edge[x].gamma | domains.edge[x].corner[0] | domains.edge[x].corner[1], FunctionSpaceTypeForm0::HierarchicalH1,order)); + fieldBucket.push_back(phi[x].back()); + } + } + + formulation.galerkin( im * k * exp2 * dof(u), tf(u), domains.gammaPhy, gauss); + for(int x = 0; x < 4; ++x) { + for(int i = 0; i < padeOrder; ++i) { + // Boundary conditions + formulation.galerkin( im * k * exp2 * coef * c[i] * dof(*phi[x][i]), tf(u), domains.edge[x].gamma, gauss); + formulation.galerkin( im * k * exp2 * coef * c[i] * dof(u), tf(u), domains.edge[x].gamma, gauss); + + // Auxiliary equations + formulation.galerkin(- (beta*beta) * grad(dof(*phi[x][i])), grad(tf(*phi[x][i])), domains.edge[x].gamma, gauss); + formulation.galerkin(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[x][i])), tf(*phi[x][i]), domains.edge[x].gamma, gauss); + formulation.galerkin( (k*k) * (exp1 * c[i] + 1.) * dof(*phi[x][i]), tf(*phi[x][i]), domains.edge[x].gamma, gauss); + formulation.galerkin( (k*k) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[x][i]), domains.edge[x].gamma, gauss); + } + } + + bool withCorner = true; + gmshFem.userDefinedParameter(withCorner, "withCorner"); + std::string CornerStrategy = "Sommerfeld"; // Sommerfeld or HABC + gmshFem.userDefinedParameter(CornerStrategy, "CornerStrategy"); + + if(withCorner) { + msg::info << "with corner treatment" << msg::endl; + msg::info << "use " << CornerStrategy << " condition at corners" << msg::endl; + // define the corner fields + std::vector< std::vector< std::vector< Field< std::complex< double >, Form::Form0 > * > > > corner_phi(4); + for(int k = 0; k < 4; ++k) { + corner_phi[k].resize(padeOrder); + for(int i = 0; i < padeOrder; ++i) { + for(int j = 0; j < padeOrder; ++j) { + corner_phi[k][i].push_back(new Field< std::complex< double >, Form::Form0 >("corner_phi_(" + std::to_string(i) + ", " + std::to_string(j) + ")^(" + std::to_string(k) + ")", domains.corner[k], FunctionSpaceTypeForm0::HierarchicalH1,order)); + fieldBucket.push_back(corner_phi[k][i].back()); + } + } + } + + for(int ci = 0; ci < 4; ++ci) { + int x = ci; + int y = (ci - 1 < 0 ? 4 - 1 : ci - 1); + for(int i = 0; i < padeOrder; ++i) { + if (CornerStrategy == "Sommerfeld"){ + // strategy 1 - Sommerfeld condition at corners + formulation.galerkin(- im * k * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss); + formulation.galerkin(- im * k * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss); + } + + else if (CornerStrategy == "HABC") { + // strategy 2 - HABC at corners + formulation.galerkin(- im * k * exp2 * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss); + formulation.galerkin(- im * k * exp2 * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss); + for(int j = 0; j < padeOrder; ++j) { + // Auxiliary equations for corner + formulation.galerkin( (exp1*c[i] + exp1*c[j] + 1.) * dof(*corner_phi[ci][i][j]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss); + formulation.galerkin( exp1 * (c[j] + 1.) * dof(*phi[x][i]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss); + formulation.galerkin( exp1 * (c[i] + 1.) * dof(*phi[y][j]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss); + + // Corner conditions + formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss); + formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss); + formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*corner_phi[ci][i][j]), tf(*phi[x][i]), domains.corner[ci], gauss); + formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*corner_phi[ci][j][i]), tf(*phi[y][i]), domains.corner[ci], gauss); + } + } + else { + msg::error << "Corner strategy not available !" << msg::endl; + exit(0); + } + } + } + } + } + } + else if (abcName == "pml") { + msg::info << "Use a PML with " << Npml << " layers" << msg::endl; + double Sigma0 = beta; + ScalarFunction< std::complex< double > > det_J; + TensorFunction< std::complex< double > > J_PML_inv_T; + const double Wpml = Rpml - R; + if (geometry == 0) { + msg::info << "stabilized PML - circular domain" << msg::endl; + msg::info << "use unbounded absorbing function with sigma_0 = " << Sigma0 << msg::endl; + ScalarFunction< std::complex< double > > cosT = x< std::complex< double > >() / r2d< std::complex< double > >(); + ScalarFunction< std::complex< double > > sinT = y< std::complex< double > >() / r2d< std::complex< double > >(); + ScalarFunction< std::complex< double > > dampingProfileR = Sigma0 / (Wpml - (r2d< std::complex< double > >() - R)); + ScalarFunction< std::complex< double > > dampingProfileInt = -Sigma0*ln((Wpml - (r2d< std::complex< double > >() - R)) / Wpml); + ScalarFunction< std::complex< double > > gamma = 1. - im * dampingProfileR / k; + ScalarFunction< std::complex< double > > gamma_hat = 1. - im * (1. / r2d< std::complex< double > >()) * dampingProfileInt / k; + det_J = gamma * gamma_hat; + J_PML_inv_T = tensor< std::complex <double > >(cosT/gamma,sinT/gamma,0.,-sinT/gamma_hat,cosT/gamma_hat,0.,0.,0.,0.); + } + else if (geometry == 1) { + msg::info << "stabilized PML - square domain" << msg::endl; + msg::info << "use unbounded absorbing function with sigma_0 = " << Sigma0 << msg::endl; + ScalarFunction< std::complex< double > > SigmaX = heaviside(abs(x< std::complex< double > >()) - R ) * Sigma0/(Rpml - abs(x< std::complex< double > >()) ); + ScalarFunction< std::complex< double > > gammaX = 1-(im/k)*SigmaX; + ScalarFunction< std::complex< double > > SigmaY = heaviside(abs(y< std::complex< double > >()) - R ) * Sigma0/(Rpml - abs(y< std::complex< double > >()) ); + ScalarFunction< std::complex< double > > gammaY = 1-(im/k)*SigmaY; + det_J = gammaX * gammaY; + J_PML_inv_T = tensor< std::complex <double > >(1./gammaX,0.,0.,0.,1./gammaY,0.,0.,0.,0.); + + /* + msg::info << "classical PML - square domain - unstable ! " << msg::endl; + formulation.galerkin(vector< std::complex< double > >((gammaY/gammaX)*(1-Mx*Mx), -Mx*My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin(vector< std::complex< double > >(-Mx*My, (gammaX/gammaY)*(1-My*My), 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin(dof(u), vector< std::complex< double > >(-gammaY*im*k*Mx,-gammaX*im*k*My, 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin(vector< std::complex< double > >(gammaY*im*k*Mx, gammaX*im*k*My, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss); + formulation.galerkin(- k*k* gammaX*gammaY * dof(u), tf(u), domains.omegaPml, gauss); + */ + + /* + msg::info << "Alternative stable PML - square domain only " << msg::endl; + ScalarFunction< std::complex< double > > Ax = (gammaY/gammaX)*(1-Mx*Mx); + std::complex< double > Cx =-im*k*Mx*( 1 - (My*My)/(1-Mx*Mx) ) / (beta*beta); + double Ka = -Mx*My/(1-Mx*Mx); + ScalarFunction< std::complex< double > > Ay = (gammaX/gammaY)*( 1 - (My*My)/(1-Mx*Mx) ); + std::complex< double > Cy = -im*k*My/(beta*beta); + + // modified bilinear term in dx + formulation.galerkin( Ax * vector< std::complex< double > >(1. , Ka, 0.) * grad(dof(u)), vector< std::complex< double > >(1., Ka, 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin( Ax * Cx * dof(u), vector< std::complex< double > >(1., Ka, 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin( Ax * (-Cx) * vector< std::complex< double > >(1. , Ka, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss); + formulation.galerkin( Ax * Cx * (-Cx) *dof(u), tf(u), domains.omegaPml, gauss); + // modified bilinear term in dy + formulation.galerkin( Ay * vector< std::complex< double > >(0. , 1., 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin( Ay * Cy * dof(u), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin( Ay * vector< std::complex< double > >(0. , -Cy, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss); + formulation.galerkin( Ay * Cy * (-Cy) *dof(u), tf(u), domains.omegaPml, gauss); + // + formulation.galerkin(- k*k/(beta*beta) * gammaX*gammaY * dof(u), tf(u), domains.omegaPml, gauss); + */ + } + // build vector and matrix for the weak form + VectorFunction< std::complex< double > > J_PML_inv_T_M = J_PML_inv_T * MM; + TensorFunction< std::complex< double > > J_PML_Linv = J_PML_inv_T * Linv; + // stabilized PML weak form - valid for a general domain + formulation.galerkin( det_J * J_PML_Linv*grad(dof(u)) , J_PML_Linv*grad(tf(u)) , domains.omegaPml, gauss); + + formulation.galerkin(+ k*k/(beta*beta) * det_J * J_PML_inv_T_M * dof(u) , J_PML_inv_T_M * tf(u), domains.omegaPml, gauss); + formulation.galerkin(+ im*k/beta* det_J * J_PML_Linv * grad(dof(u)), J_PML_inv_T_M * tf(u), domains.omegaPml, gauss); + + formulation.galerkin(- im*k/beta* det_J * J_PML_inv_T_M * dof(u), J_PML_Linv * grad(tf(u)), domains.omegaPml, gauss); + formulation.galerkin(- k*k/(beta*beta) * det_J * dof(u), tf(u), domains.omegaPml, gauss); + } + else { + msg::error << "ABC Type not available !" << msg::endl; + exit(0); + } + + // solve + formulation.pre(); + formulation.assemble(); + formulation.solve(); + + // compute the projection of the analytic solution on the FEM basis - best FEM approximation + Formulation< std::complex< double > > projection("helmholtzflow"); + Field< std::complex< double >, Form::Form0 > uP("uP", domains.omega, FunctionSpaceTypeForm0::HierarchicalH1, order); + projection.galerkin( -dof(uP), tf(uP), domains.omega, gauss); + projection.galerkin( *solution, tf(uP), domains.omega, gauss); + + msg::info << "Compute best approximation... " << msg::endl; + projection.pre(); + projection.assemble(); + projection.solve(); + + // save results and compute errors + save(+u, domains.omega, "u"); + save(+uP, domains.omega, "uP"); + if (abcName == "pml"){ + save(+u, domains.omegaPml, "u_pml"); + } + save(*solution, domains.omega, "u_exact"); + save(*solution - u, domains.omega, "error"); + + std::complex< double > num = integrate(pow(abs(*solution - u), 2), domains.omega, gauss); + std::complex< double > numP = integrate(pow(abs(*solution - uP), 2), domains.omega, gauss); + std::complex< double > den = integrate(pow(abs(*solution), 2), domains.omega, gauss); + msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl; + msg::info << "Best L_2 error = " << 100.*sqrt(numP / den) << " %" << msg::endl; + msg::warning << "L_2 error should exclude neighbourhood elements to the point source" << msg::endl; + + for(unsigned int i = 0; i < fieldBucket.size(); ++i) { + delete fieldBucket[i]; + } + + return 0; +} diff --git a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_circle.sh b/examples/helmholtzFlow/freefield/runTests_circle.sh similarity index 72% rename from examples/helmholtzFlow/helmholtzflow/FreeField/runTests_circle.sh rename to examples/helmholtzFlow/freefield/runTests_circle.sh index 071aef80cbfed92bf720caf73de6441a4c40547b..38df525d983669110137250d7a282618aff59345 100644 --- a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_circle.sh +++ b/examples/helmholtzFlow/freefield/runTests_circle.sh @@ -6,18 +6,18 @@ do echo "****************************" echo " Test zeroth order ABC " echo "****************************" - ./demo -abcName ABC-0 -M $i -geometry 0 + ./example -abcName ABC-0 -M $i -geometry 0 echo "****************************" echo " Test second order ABC " echo "****************************" - ./demo -abcName ABC-2 -M $i -geometry 0 + ./example -abcName ABC-2 -M $i -geometry 0 echo "****************************" echo " Test Pade type ABC " echo "****************************" - ./demo -abcName Pade -M $i -geometry 0 + ./example -abcName Pade -M $i -geometry 0 echo "****************************" echo " Test 4 layers pml " echo "****************************" - ./demo -abcName pml -Npml 4 -M $i -geometry 0 + ./example -abcName pml -Npml 4 -M $i -geometry 0 done diff --git a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_square.sh b/examples/helmholtzFlow/freefield/runTests_square.sh similarity index 63% rename from examples/helmholtzFlow/helmholtzflow/FreeField/runTests_square.sh rename to examples/helmholtzFlow/freefield/runTests_square.sh index 4f9d23aeb4fd310309709067db718a473e0392b3..8994a558e3ed19e9238dec613e4e7e63854041de 100644 --- a/examples/helmholtzFlow/helmholtzflow/FreeField/runTests_square.sh +++ b/examples/helmholtzFlow/freefield/runTests_square.sh @@ -6,19 +6,19 @@ do echo "****************************" echo " Test zeroth order ABC " echo "****************************" - ./demo -abcName ABC-0 -M $i -geometry 1 + ./example -abcName ABC-0 -M $i -geometry 1 echo "****************************" echo " Test second order ABC " echo "****************************" - ./demo -abcName ABC-2 -M $i -geometry 1 + ./example -abcName ABC-2 -M $i -geometry 1 echo "****************************" echo " Test Pade type ABC " echo "****************************" - ./demo -abcName Pade -M $i -geometry 1 -CornerStrategy Sommerfeld - ./demo -abcName Pade -M $i -geometry 1 -CornerStrategy HABC + ./example -abcName Pade -M $i -geometry 1 -CornerStrategy Sommerfeld + ./example -abcName Pade -M $i -geometry 1 -CornerStrategy HABC echo "****************************" echo " Test 4 layers pml " echo "****************************" - ./demo -abcName pml -Npml 4 -M $i -geometry 1 + ./example -abcName pml -Npml 4 -M $i -geometry 1 done diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt b/examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt deleted file mode 100644 index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..0000000000000000000000000000000000000000 --- a/examples/helmholtzFlow/helmholtzDuctHeteroX/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -cmake_minimum_required(VERSION 3.0 FATAL_ERROR) - -project(demo CXX) - -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) - -add_executable(demo main.cpp ${EXTRA_INCS}) -target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt b/examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt deleted file mode 100644 index 4d2e4b0c91d2a7e0da5a224e581c56eca24b2f71..0000000000000000000000000000000000000000 --- a/examples/helmholtzFlow/helmholtzDuctHeteroY/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -cmake_minimum_required(VERSION 3.0 FATAL_ERROR) - -project(demo CXX) - -include(${CMAKE_CURRENT_SOURCE_DIR}/../../demos.cmake) - -add_executable(demo main.cpp ${EXTRA_INCS}) -target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt b/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt deleted file mode 100644 index f01d5714177b61745622f922015d7abd78c1cdcf..0000000000000000000000000000000000000000 --- a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -cmake_minimum_required(VERSION 3.0 FATAL_ERROR) - -project(demo CXX) - -include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake) - -add_executable(demo main.cpp ${EXTRA_INCS}) -target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt b/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt deleted file mode 100644 index f01d5714177b61745622f922015d7abd78c1cdcf..0000000000000000000000000000000000000000 --- a/examples/helmholtzFlow/helmholtzflow/FlowPastCylinder/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -cmake_minimum_required(VERSION 3.0 FATAL_ERROR) - -project(demo CXX) - -include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake) - -add_executable(demo main.cpp ${EXTRA_INCS}) -target_link_libraries(demo ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/waveguide/Readme.md b/examples/helmholtzFlow/waveguide/Readme.md new file mode 100644 index 0000000000000000000000000000000000000000..441fbd98b5877dc93be0f92992b9afe2c34c3ed8 --- /dev/null +++ b/examples/helmholtzFlow/waveguide/Readme.md @@ -0,0 +1,78 @@ +# Assessment of Local Absorbing Boundary Conditions for Heterogeneous Convected Helmholtz Problems +> Marchner P. + +Dependency(ies): +* GmshFEM : v 1.0.0 +* Gmsh : v 4.8.2 + +## Problem description + +This code solves a convected Helmholtz problem in a 2d straight waveguide +$$\frac{D_0}{Dt}\left(\frac{1}{c_0^2}\frac{D_0}{Dt}\right) u - \rho_0^{-1}\nabla \cdot (\rho_0 \nabla ) u = 0, \quad \text{in} \; \Omega=[0, L]\times[0, H]$$ +with either a uniform (`uniform/`) or non-uniform flow (`nonUniform/`) + +Different ABCs are tested on the output boundary. +For the non-uniform flow, the relative $$L^2$$-error is computed thanks to a PML reference solution. + +## Installation + +``` + mkdir build && cd build + cmake .. + make +``` + +## Usage + +Simply run: +``` + ./example [PARAM] +``` +with `[PARAM]`: +* `-h [x]` where `[x]` is the mesh size. +* `-L [x]` where `[x]` is the waveguide length and ABC position. +* `-H [x]` where `[x]` is the waveguide height. +* `-FEMorder [x]` where `[x]` is the order of the finite element hierarchical basis. +* `-omega [x]` where `[x]` is the angular frequency for (`waveguideUniform/`) or `-kinf [x]` where `[x]` is the reference wavenumber for (`waveguideNonuniform/`). +* `-mode [x]` where `[x]` is the mode number for the input boundary condition. +* `-abcName [x]` where `[x]` is the ABC imposed on the right boundary. +* `-padeOrder [x]` is the Pade order if `abcName = Pade`. +* `-angle [x]` is the rotation angle of the branch-cut if `abcName = Pade`. +* `-runcase [x]` selects a specific configuration for the results reproductibility. + +Each subfolder contains +* the main script (`main.cpp`) +* a loop for testing the ABCs over a given frequency range (`main_omega.cpp`) +To run one of the loops just change the name in the `CMakeLists.txt` and recompile + +## Results reproductibility +##### Folder `uniform/` + +Figure 6 can be obtained with `main_omega.cpp` and +``` + ./example -runcase [i] +``` +with i=1 or i=2 for respectively Figures 6a and 6b. + +##### Folder `nonUniform/` +Pressure maps from Figure 8 can be obtained by compiling `main.cpp` and running +``` + ./example -kinf 40 +``` +The solution `vpml.msh` can be open in Gmsh. + +For Figure 9, compile `main_omega.cpp` and run the command +``` +./example -runcase 1 -L [i] +``` +with i={0.75,1,1.25,2} to obtain the 4 sub-figures. + +For Figure 10, compile `main_omega.cpp` and run the command +``` +./example -runcase 2 -L [i] +``` +with i={0.75,1,1.25,2} to obtain the 4 sub-figures. + +## Reference + +// Preprint diff --git a/examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt b/examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6f8caaab62b9209c4cd7f98d220eaa43d6841dbe --- /dev/null +++ b/examples/helmholtzFlow/waveguide/nonUniform/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(example CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake) + +add_executable(example main.cpp ${EXTRA_INCS}) +target_link_libraries(example ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main.cpp b/examples/helmholtzFlow/waveguide/nonUniform/main.cpp similarity index 99% rename from examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main.cpp rename to examples/helmholtzFlow/waveguide/nonUniform/main.cpp index b153a8e9d274fa996e87e6830a8151bafedc158d..8543cdf127301afbacfc15a92a1ab50b649b1394 100644 --- a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main.cpp +++ b/examples/helmholtzFlow/waveguide/nonUniform/main.cpp @@ -113,7 +113,6 @@ int main(int argc, char **argv) int FEMorder = 8; gmshFem.userDefinedParameter(FEMorder, "FEMorder"); std::string gauss = "Gauss"+std::to_string(2*FEMorder+1); - gmshFem.userDefinedParameter(gauss, "gauss"); // meshing double L = 2.; // duct length @@ -279,7 +278,7 @@ int main(int argc, char **argv) // Compute reference solution with large PML msg::info << "Use a PML with N = " << Npml << " layers" << msg::endl; const double Lpml = Lext+Npml*h; - ScalarFunction< std::complex< double > > Sigma0 = 1;//4*betax*betax; + ScalarFunction< std::complex< double > > Sigma0 = 4; // Bermudez function ScalarFunction< std::complex< double > > SigmaX = Sigma0 / (Lpml - x< std::complex< double > >() ); ScalarFunction< std::complex< double > > gammaX = 1-(im/k0)*SigmaX; diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main_loop_omega.cpp b/examples/helmholtzFlow/waveguide/nonUniform/main_omega.cpp similarity index 100% rename from examples/helmholtzFlow/helmholtzflow/Duct_HeteroX/main_loop_omega.cpp rename to examples/helmholtzFlow/waveguide/nonUniform/main_omega.cpp diff --git a/examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt b/examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6f8caaab62b9209c4cd7f98d220eaa43d6841dbe --- /dev/null +++ b/examples/helmholtzFlow/waveguide/uniform/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) + +project(example CXX) + +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../example.cmake) + +add_executable(example main.cpp ${EXTRA_INCS}) +target_link_libraries(example ${EXTRA_LIBS}) diff --git a/examples/helmholtzFlow/waveguide/uniform/main.cpp b/examples/helmholtzFlow/waveguide/uniform/main.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7522d00b3b30b3829f1d6740637f8e82d6a3c43e --- /dev/null +++ b/examples/helmholtzFlow/waveguide/uniform/main.cpp @@ -0,0 +1,266 @@ +#include "GmshFem.h" +#include "Formulation.h" +#include "gmsh.h" +#include "AnalyticalFunction.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::analytics; +using namespace gmshfem::post; +using namespace gmshfem::equation; + +void meshDuct(const double L, const double H, const double h, const int MeshOrder) +{ + msg::info << " - Duct of length L = " << L << msg::endl; + msg::info << " - Duct of heigth H = " << H << msg::endl; + gmsh::model::add("geometry"); + + // duct definition + gmsh::model::geo::addPoint(L, H, 0., h, 6); + gmsh::model::geo::addPoint(0, H, 0., h, 7); + gmsh::model::geo::addPoint(0, 0, 0., h, 8); + gmsh::model::geo::addPoint(L, 0, 0., h, 9); + + gmsh::model::geo::addLine(6, 7, 5); + gmsh::model::geo::addLine(7, 8, 6); + gmsh::model::geo::addLine(8, 9, 7); + gmsh::model::geo::addLine(9, 6, 8); + gmsh::model::geo::addCurveLoop({5,6,7,8}, 1); + gmsh::model::geo::addPlaneSurface({1}, 1); + gmsh::model::geo::mesh::setTransfiniteSurface(1); + gmsh::model::geo::mesh::setRecombine(1,1); + + // physical groups + gmsh::model::addPhysicalGroup(2, {1}, 1); + gmsh::model::setPhysicalName(2, 3, "omega"); + gmsh::model::addPhysicalGroup(1, {5}, 2); + gmsh::model::setPhysicalName(1, 2, "gammaTop"); + gmsh::model::addPhysicalGroup(1, {6}, 3); + gmsh::model::setPhysicalName(1, 3, "gammaLeft"); + gmsh::model::addPhysicalGroup(1, {7}, 4); + gmsh::model::setPhysicalName(1, 4, "gammaBottom"); + gmsh::model::addPhysicalGroup(1, {8}, 5); + gmsh::model::setPhysicalName(1, 5, "gammaRight"); + + gmsh::model::addPhysicalGroup(0, {6}, 6); + gmsh::model::setPhysicalName(0, 6, "gammaTopCorner"); + gmsh::model::addPhysicalGroup(0, {9}, 7); + gmsh::model::setPhysicalName(0, 7, "gammaBottomCorner"); + + // generate mesh + gmsh::model::geo::synchronize(); + gmsh::model::mesh::generate(2); + gmsh::model::mesh::setOrder(MeshOrder); + gmsh::option::setNumber("Mesh.ElementOrder",MeshOrder); + gmsh::model::mesh::setOrder(MeshOrder); +} + +int main(int argc, char **argv) +{ + //***** + // Problem declaration + //***** + GmshFem gmshFem(argc, argv); + const double pi = 3.14159265358979323846264338327950288; + const std::complex< double > im(0., 1.); + + // numerical parameters + int FEMorder = 8; + gmshFem.userDefinedParameter(FEMorder, "FEMorder"); + std::string gauss = "Gauss"+std::to_string(2*FEMorder+1); + + // meshing + double L = 0.5; // duct length + gmshFem.userDefinedParameter(L, "L"); + double H = 0.25; // duct height + gmshFem.userDefinedParameter(H, "H"); + double h = 1./40; // meshsize + gmshFem.userDefinedParameter(h, "h"); + meshDuct(L, H, h, 1); + + // duct problem specifications + std::string problem = "hard"; + gmshFem.userDefinedParameter(problem, "problem"); + double w = 30; // frequency + gmshFem.userDefinedParameter(w, "omega"); + std::vector< unsigned int > n = {3}; // single mode + // std::vector< unsigned int > n = {0,2,4,6,8}; // multiple modes modes + std::vector< double > A(n.size(),1.); // select amplitudes + int N_modes = n.size(); + msg::info << " There are " << N_modes << " input modes " << msg::endl; + + // Choose ABC type + std::string abcName = "Pade"; + gmshFem.userDefinedParameter(abcName, "abcName"); + + // mean flow + const double M = 0.8; + double beta = sqrt(1-M*M); + + std::vector< std::complex< double > > KX(N_modes); + std::vector< double > KY(N_modes), SIGMA(N_modes); + for(unsigned int j = 0; j < n.size(); ++j) { + KY[j] = n[j]*pi/H; + SIGMA[j] = w*w-beta*beta*KY[j]*KY[j]; + if ( SIGMA[j] >=0 ) { + KX[j] = (1/(beta*beta))*(-M*w + sqrt(SIGMA[j])); + msg::info << " Mode " << n[j] << " is propagative " << msg::endl; + if ( w>(beta*KY[j]) && w<KY[j] ) { + msg::info << " Mode " << n[j] << " is inverse upstream ! " << msg::endl; + } + } + else { + KX[j] = (1/(beta*beta))*(-M*w - im* sqrt(abs(SIGMA[j]))); + msg::info << " Mode " << n[j] << " is evanescent " << msg::endl; + } + msg::info << " - propagating wavenumber, kx = " << KX[j] << " of mode " << n[j] << msg::endl; + } + + float PtsByWl = 2*pi*FEMorder*(1+M) / w / h; + msg::info << "********************************" << msg::endl; + msg::info << " - convected Helmholtz duct problem " << msg::endl; + msg::info << " - wave number = " << w << msg::endl; + msg::info << " - FEM basis order = " << FEMorder << "" << msg::endl; + msg::info << " - Approximate dofs by wavelength = " << PtsByWl << "" << msg::endl; + msg::info << "********************************" << msg::endl; + + // declare formulation and FEM domains + std::vector< FieldInterface< std::complex< double > > * > fieldBucket; + Formulation< std::complex< double > > formulation("helmholtzflow"); + + Domain omega(2, 1); + Domain gammaTop(1, 2); + Domain gammaLeft(1, 3); + Domain gammaBottom(1, 4); + Domain gammaRight(1, 5); + Domain gammaTopCorner(0, 6); + Domain gammaBottomCorner(0, 7); + + Field< std::complex< double >, Form::Form0 > v("v", omega | gammaLeft | gammaRight | gammaTop | gammaBottom, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder); + // define analytical solution + Function< std::complex< double >, Degree::Degree0 > *solution = nullptr; + + if(problem == "soft") { + v.addConstraint(gammaTop, 0.); + v.addConstraint(gammaBottom, 0.); + solution = new AnalyticalFunction< helmholtz2D::DuctModeSolutionMulti< std::complex< double > > >(w, M, H, A, n, 0., 0., 1); + } + else { + solution = new AnalyticalFunction< helmholtz2D::DuctModeSolutionMulti< std::complex< double > > >(w, M, H, A, n, 0., 0., 0); + } + solution->activateMemorization(); + + formulation.galerkin(vector< std::complex< double > >(beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss); + formulation.galerkin(vector< std::complex< double > >(0., 1., 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omega, gauss); + formulation.galerkin(vector< std::complex< double > >(im*w*M, 0., 0.) * grad(dof(v)), tf(v), omega, gauss); + formulation.galerkin(-im*w*M*dof(v), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss); + formulation.galerkin(- w * w * dof(v), tf(v), omega, gauss); + + // input BC + ScalarFunction< std::complex< double > > Input_modal_sum=0.; + for(unsigned int j = 0; j < n.size(); ++j) { + if(problem == "soft"){ + Input_modal_sum = Input_modal_sum + A[j] * KX[j] * sin( KY[j] * y< std::complex< double > >() ); + } + else if(problem == "hard"){ + Input_modal_sum = Input_modal_sum + A[j] * KX[j] * cos( KY[j] * y< std::complex< double > >() ); + } + else { + msg::error << " this problem is not defined " << msg::endl; + } + } + formulation.galerkin(-beta*beta * im * Input_modal_sum , tf(v), gammaLeft, gauss); + formulation.galerkin(-im * w * M * dof(v), tf(v), gammaLeft, gauss); + + // output BC + if(abcName == "ABC-0") { + msg::info << "Use Sommerfeld ABC." << msg::endl; + const double angle = 0.; + std::complex< double > wM = w*(M-exp(im*angle/2.)) / (beta*beta); //angle=0 -> -k/(1+M); + formulation.galerkin(-beta*beta * im * wM * dof(v), tf(v), gammaRight, gauss); + formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss); + } + else if(abcName == "ABC-2") { + msg::info << "Use second order ABC." << msg::endl; + const double angle = 0.; + // contribution to the mass and rigidity matrices + std::complex< double > wM = w*(M-cos(angle/2.)) / (beta*beta); + std::complex< double > wK = -exp(-im*angle/2.)/(2*w); + formulation.galerkin(-beta*beta * im * wM * dof(v), tf(v), gammaRight, gauss); + formulation.galerkin( beta*beta * im * wK * grad(dof(v)), grad(tf(v)), gammaRight, gauss); + formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss); + + /* use DtN eigenvalues instead of Laplace Beltrami operator + std::complex< double > kT = -k + KY[0]*KY[0]/(2*k); + formulation.galerkin(-beta*beta * im * kT * dof(v), tf(v), gammaRight, gauss); + */ + } + else if(abcName == "Pade") { + int padeOrder = 4; + gmshFem.userDefinedParameter(padeOrder, "padeOrder"); + double angle = -pi/4.; + gmshFem.userDefinedParameter(angle, "angle"); + msg::info << "Use Pade ABC of order " << padeOrder << " with angle " << angle << " rad"<< msg::endl; + + const double Np = 2. * padeOrder + 1.; + const std::complex< double > exp1 = std::complex<double>(std::cos(angle),std::sin(angle)); + const std::complex< double > exp2 = std::complex<double>(std::cos(angle/2.),std::sin(angle/2.)); + const std::complex< double > coef = 2./Np; + std::vector< std::complex< double > > c(padeOrder, 0.); + for(int i = 0; i < padeOrder; ++i) { + c[i] = std::tan((i + 1) * pi / Np); + c[i] *= c[i]; + } + + // define the auxiliary fields + std::vector< Field< std::complex< double >, Form::Form0 >* > phi; + Field< std::complex< double >, Form::Form0 >* phis; + for(int i = 0; i < padeOrder; ++i) { + phis = new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i), gammaRight | gammaTopCorner | gammaBottomCorner, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder); + if(problem == "soft") { // impose Dirichlet constraints on the auxiliary fields + phis->addConstraint(gammaTopCorner,0.); + phis->addConstraint(gammaBottomCorner,0.); + } + phi.push_back(phis); + fieldBucket.push_back(phi.back()); + } + + // write the augmented weak form + formulation.galerkin(im * w * exp2 * dof(v), tf(v), gammaRight, gauss); + for(int i = 0; i < padeOrder; ++i) { + // boundary integral terms relating the auxiliary fields + formulation.galerkin(im * w * exp2 * coef * c[i] * dof(*phi[i]), tf(v), gammaRight, gauss); + formulation.galerkin(im * w * exp2 * coef * c[i] * dof(v), tf(v), gammaRight, gauss); + + // coupling of the auxiliary equations + formulation.galerkin( grad(dof(*phi[i])), grad(tf(*phi[i])), gammaRight, gauss); + formulation.galerkin(- (w*w)/(beta*beta) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), gammaRight, gauss); + formulation.galerkin(- (w*w)/(beta*beta) * exp1 * (c[i] + 1.) * dof(v), tf(*phi[i]), gammaRight, gauss); + } + } + else if(abcName == "DtN") { + msg::info << "Use analytical DtN." << msg::endl; + formulation.galerkin(beta*beta * im * KX[0] * dof(v), tf(v), gammaRight, gauss); + formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss); + } + + formulation.pre(); + formulation.assemble(); + formulation.solve(); + save(+v, omega, "v"); + save(*solution, omega, "v_exact"); + save(*solution - v, omega, "error"); + std::complex< double > num = integrate(pow(abs(*solution - v), 2), omega, gauss); + std::complex< double > den = integrate(pow(abs(*solution), 2), omega, gauss); + msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl; + + for(unsigned int i = 0; i < fieldBucket.size(); ++i) { + delete fieldBucket[i]; + } + + return 0; +} diff --git a/examples/helmholtzFlow/helmholtzflow/Duct_ABCs/main_omega.cpp b/examples/helmholtzFlow/waveguide/uniform/main_omega.cpp similarity index 100% rename from examples/helmholtzFlow/helmholtzflow/Duct_ABCs/main_omega.cpp rename to examples/helmholtzFlow/waveguide/uniform/main_omega.cpp