Boundary conditions on SimpleWave

common/configuration.h

@@ -6,6 +6,7 @@
 #include "Domain.h"
 //Standard Library
 #include <array>
+#include <map>
 //GmshFWI Library
 #include "model/element.h"
 #include "model/parametrization.h"
@@ -21,6 +22,7 @@ protected:
     std::array<gmshfem::domain::Domain,2> _wave_omega;
     std::array<gmshfem::domain::Domain,2> _model_known;
     std::array<gmshfem::domain::Domain,2> _model_unknown;
+    std::map<std::string, gmshfem::domain::Domain> _named_domains;
     gmshfem::domain::Domain _data_omega;
 
     unsigned int _np;
@@ -74,6 +76,7 @@ public:
     gmshfem::domain::Domain model_evaldomain() const {return _model_unknown[Support::BLK] | _model_unknown[Support::BND];};
     gmshfem::domain::Domain wave_evaldomain() const {return _wave_omega[Support::BLK]  | _wave_omega[Support::BND] | _points;};
     gmshfem::domain::Domain data_evaldomain() const {return _data_omega;};
+    const std::map<std::string, gmshfem::domain::Domain>& named_domains() const {return _named_domains;}
 
     gmshfem::domain::Domain emitter(unsigned int shot, unsigned int em) const;
     gmshfem::domain::Domain receiver(unsigned int shot, unsigned int rec) const;

specific/configuration/surface_acquisition.cpp

@@ -113,6 +113,8 @@ namespace surface_acquisition
         _supersurface[Support::BLK] = Domain(2, 22);
         _supersurface[Support::BND] = Domain(1, 12);
 
+        _named_domains["top_bnd"] = Domain("top_bnd");
+
         _wave_omega[Support::BLK] = _subsurface[Support::BLK] | _supersurface[Support::BLK];
         _wave_omega[Support::BND] = _subsurface[Support::BND] | _supersurface[Support::BND];
 
@@ -341,6 +343,7 @@ namespace surface_acquisition
                 lymin_reverse[n] = -lymin[n];
             }
 
+            // Nomenclature: line from one y to another, left or right. All lines go up
             int lHyminl = factory::addLine(pHl,pyminl);
             int lymin0l = factory::addLine(pyminl,p0l);
             int l0ymaxl = factory::addLine(p0l,pymaxl);
@@ -430,6 +433,10 @@ namespace surface_acquisition
             gmodel::addPhysicalGroup(1, lsuper, 12);
             gmodel::setPhysicalName(1, 12, "supersurface_bnd");
 
+            // Add another physical groupe for the top boundary, to represent e.g. reflection at the air-water interface
+            gmodel::addPhysicalGroup(1, lymax, 13);
+            gmodel::setPhysicalName(1, 13, "top_bnd");
+
             gmodel::addPhysicalGroup(2, {s_sub}, 21);
             gmodel::setPhysicalName(2, 21, "subsurface_vol");
             gmodel::addPhysicalGroup(2, {s_supermin,s_supermax}, 22);
@@ -507,6 +514,10 @@ namespace surface_acquisition
             gmodel::addPhysicalGroup(1, lsuper, 12);
             gmodel::setPhysicalName(1, 12, "supersurface_bnd");
 
+            // Add another physical groupe for the top boundary, to represent e.g. reflection at the air-water interface
+            gmodel::addPhysicalGroup(1, lymax, 13);
+            gmodel::setPhysicalName(1, 13, "top_bnd");
+
             gmodel::addPhysicalGroup(2, {s_sub}, 21);
             gmodel::setPhysicalName(2, 21, "subsurface_vol");
             gmodel::addPhysicalGroup(2, {s_supermax}, 22);

specific/wave/equation/simplewave.cpp

@@ -72,6 +72,15 @@ namespace simplewave
         {
                 throw Exception("Impossible to compute boundary terms with real impedance: sensitivities are wrong.");
         }
+
+        if (gmshFem.userDefinedParameter(_topBC, "top_bc")) {
+            if (_topBC != "Absorbing" && _topBC != "Dirichlet" && _topBC != "Neumann") {
+                throw Exception("Invalid value of parameter top_bc: " + _topBC);
+            }
+        }
+        else {
+            msg::warning << "SimpleWave without specified top_bc defaults to Absorbing." << msg::endl;
+        }
     }
 
     bool Equation::modelIsObsolete()
@@ -104,16 +113,35 @@ namespace simplewave
     }
     void Equation::setLHS(const ModelStateEvaluator& ms)
     {
+        auto top = _config->named_domains().find("top_bnd");
+        if (_topBC == "Dirichlet" || _topBC == "Neumann")
+        {
+            if (top == _config->named_domains().end())
+                throw Exception("Asked for specific BCs on top but this entity isn't defined.");
+        }
+
+        if (_topBC == "Dirichlet")
+        {
+            _v.addConstraint(top->second, 0.0);
+        }
+
         _formulation.integral(-1.*grad(dof(_v)), grad(tf(_v)), _config->wave_omega(Support::BLK),integrationType(2*_w_degree));
         _formulation.integral(_pulsation * _pulsation * m0 * dof(_v), tf(_v), _config->model_known(Support::BLK), integrationType(_integrationDegreeBlk));
         
+        auto absorbingPartKnown = _config->model_known(Support::BND);
+        auto absorbingPartUnknown = _config->model_unknown(Support::BND);
+        if (_topBC == "Neumann") {
+            absorbingPartKnown &= ~(top->second);
+            absorbingPartUnknown &= ~(top->second);
+        }
+
         if(_realImpedance)
         {
-            _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), _config->model_known(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+            _formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
         }
         else
         {
-            _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), _config->model_known(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+            _formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
         }
         if(!_config->model_unknown(Support::BLK).isEmpty())
         {
@@ -123,11 +151,11 @@ namespace simplewave
         {
             if(_realImpedance)
             {
-                _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) *  dof(_v), tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+                _formulation.integral(-1. * im * _pulsation * sqrt(real(m)) *  dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
             }
             else
             {
-                _formulation.integral(-1. * im * _pulsation * sqrt(m) *  dof(_v), tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
+                _formulation.integral(-1. * im * _pulsation * sqrt(m) *  dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
             }
         }
     }

specific/wave/equation/simplewave.h

@@ -21,6 +21,7 @@ namespace simplewave
         bool _realImpedance;
         double _const_preconditioner;
         using EquationInterface<Physic::acoustic>::_boundary;
+        std::string _topBC;
     public:
         Equation(unsigned int f_idx,double pulsation,const ConfigurationInterface* const config, const wave::Discretization<Physic::acoustic>& w_discret,const gmshfem::common::GmshFem& gmshFem, std::string suffix = "");