diff --git a/inputFiles/wavePropagation/acous3D_Q3_small_base.xml b/inputFiles/wavePropagation/acous3D_Q3_small_base.xml
index 8529d2ce19a..00d67d72bba 100644
--- a/inputFiles/wavePropagation/acous3D_Q3_small_base.xml
+++ b/inputFiles/wavePropagation/acous3D_Q3_small_base.xml
@@ -67,6 +67,14 @@
scale="1500"
setNames="{ all }"/>
+
+
+
+
-
+
-
+
@@ -168,7 +168,7 @@
right,
front,
back}"/>
-
+
+
+
+
+
+
-
+
-
+
diff --git a/inputFiles/wavePropagation/acous3D_small_base.xml b/inputFiles/wavePropagation/acous3D_small_base.xml
index d917f09b737..7c2cec72128 100644
--- a/inputFiles/wavePropagation/acous3D_small_base.xml
+++ b/inputFiles/wavePropagation/acous3D_small_base.xml
@@ -67,6 +67,14 @@
scale="1500"
setNames="{ all }"/>
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/wavePropagation/benchmarks/elas3D_benchmark_base.xml b/inputFiles/wavePropagation/benchmarks/elas3D_benchmark_base.xml
index a1759a7b731..a7cceeadb85 100644
--- a/inputFiles/wavePropagation/benchmarks/elas3D_benchmark_base.xml
+++ b/inputFiles/wavePropagation/benchmarks/elas3D_benchmark_base.xml
@@ -184,7 +184,7 @@
name="cellDensity"
initialCondition="1"
objectPath="ElementRegions/Region/cb"
- fieldName="mediumDensity"
+ fieldName="mediumDensityE"
scale="1"
setNames="{ all }"/>
diff --git a/inputFiles/wavePropagation/elas3D_DAS_smoke.xml b/inputFiles/wavePropagation/elas3D_DAS_smoke.xml
index 78b0210ad08..030dcc475d5 100644
--- a/inputFiles/wavePropagation/elas3D_DAS_smoke.xml
+++ b/inputFiles/wavePropagation/elas3D_DAS_smoke.xml
@@ -149,7 +149,7 @@
name="cellDensity"
initialCondition="1"
objectPath="ElementRegions/Region/cb"
- fieldName="mediumDensity"
+ fieldName="mediumDensityE"
scale="2000"
setNames="{ all }"/>
diff --git a/inputFiles/wavePropagation/elas3D_Q3_small_base.xml b/inputFiles/wavePropagation/elas3D_Q3_small_base.xml
index 04155600218..bdc3fc93480 100644
--- a/inputFiles/wavePropagation/elas3D_Q3_small_base.xml
+++ b/inputFiles/wavePropagation/elas3D_Q3_small_base.xml
@@ -163,7 +163,7 @@
name="cellDensity"
initialCondition="1"
objectPath="mesh/FE1/ElementRegions/Region/cb"
- fieldName="mediumDensity"
+ fieldName="mediumDensityE"
scale="1"
setNames="{ all }"/>
diff --git a/inputFiles/wavePropagation/elas3D_Q5_small_base.xml b/inputFiles/wavePropagation/elas3D_Q5_small_base.xml
index b811a1c3a3e..58991fa7487 100644
--- a/inputFiles/wavePropagation/elas3D_Q5_small_base.xml
+++ b/inputFiles/wavePropagation/elas3D_Q5_small_base.xml
@@ -163,7 +163,7 @@
name="cellDensity"
initialCondition="1"
objectPath="mesh/FE1/ElementRegions/Region/cb"
- fieldName="mediumDensity"
+ fieldName="mediumDensityA"
scale="1"
setNames="{ all }"/>
diff --git a/inputFiles/wavePropagation/elas3D_small_base.xml b/inputFiles/wavePropagation/elas3D_small_base.xml
index 0ed1a918810..f95c3491a3e 100644
--- a/inputFiles/wavePropagation/elas3D_small_base.xml
+++ b/inputFiles/wavePropagation/elas3D_small_base.xml
@@ -170,7 +170,7 @@
name="cellDensity"
initialCondition="1"
objectPath="ElementRegions/Region/cb"
- fieldName="mediumDensity"
+ fieldName="mediumDensityE"
scale="1"
setNames="{ all }"/>
diff --git a/src/coreComponents/codingUtilities/StringUtilities.cpp b/src/coreComponents/codingUtilities/StringUtilities.cpp
index 5ac99424e64..30976346a53 100644
--- a/src/coreComponents/codingUtilities/StringUtilities.cpp
+++ b/src/coreComponents/codingUtilities/StringUtilities.cpp
@@ -73,6 +73,9 @@ string removeStringAndFollowingContent( string const & str,
template< typename T >
string toMetricPrefixString( T const & value )
{
+ if( std::fpclassify( value ) == FP_ZERO )
+ return " 0.0 ";
+
// These are the metric prefixes corrosponding to kilo, mega, giga...etc.
char const prefixes[12] = { 'f', 'p', 'n', 'u', 'm', ' ', 'K', 'M', 'G', 'T', 'P', 'E'};
string rval;
diff --git a/src/coreComponents/fileIO/Outputs/VTKOutput.cpp b/src/coreComponents/fileIO/Outputs/VTKOutput.cpp
index aca02acc171..68d7ae80fff 100644
--- a/src/coreComponents/fileIO/Outputs/VTKOutput.cpp
+++ b/src/coreComponents/fileIO/Outputs/VTKOutput.cpp
@@ -34,8 +34,10 @@ VTKOutput::VTKOutput( string const & name,
m_plotFileRoot( name ),
m_writeFaceMesh(),
m_plotLevel(),
+ m_highOrder(),
m_onlyPlotSpecifiedFieldNames(),
m_fieldNames(),
+ m_levelNames(),
m_writer( getOutputDirectory() + '/' + m_plotFileRoot )
{
registerWrapper( viewKeysStruct::plotFileRoot, &m_plotFileRoot ).
@@ -63,10 +65,19 @@ VTKOutput::VTKOutput( string const & name,
setDescription(
"If this flag is equal to 1, then we only plot the fields listed in `fieldNames`. Otherwise, we plot all the fields with the required `plotLevel`, plus the fields listed in `fieldNames`" );
+ registerWrapper( viewKeysStruct::highOrder, &m_highOrder ).
+ setApplyDefaultValue( 1 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Should the vtk files contain the high order solution or not." );
+
registerWrapper( viewKeysStruct::fieldNames, &m_fieldNames ).
setInputFlag( InputFlags::OPTIONAL ).
setDescription( "Names of the fields to output. If this attribute is specified, GEOSX outputs all the fields specified by the user, regardless of their `plotLevel`" );
+ registerWrapper( viewKeysStruct::levelNames, &m_levelNames ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Names of mesh levels to output." );
+
registerWrapper( viewKeysStruct::binaryString, &m_writeBinaryData ).
setApplyDefaultValue( m_writeBinaryData ).
setInputFlag( InputFlags::OPTIONAL ).
@@ -85,6 +96,7 @@ void VTKOutput::postProcessInput()
{
m_writer.setOutputLocation( getOutputDirectory(), m_plotFileRoot );
m_writer.setFieldNames( m_fieldNames.toViewConst() );
+ m_writer.setLevelNames( m_levelNames.toViewConst() );
m_writer.setOnlyPlotSpecifiedFieldNamesFlag( m_onlyPlotSpecifiedFieldNames );
string const fieldNamesString = viewKeysStruct::fieldNames;
@@ -129,6 +141,7 @@ bool VTKOutput::execute( real64 const time_n,
m_writer.setOutputMode( m_writeBinaryData );
m_writer.setOutputRegionType( m_outputRegionType );
m_writer.setPlotLevel( m_plotLevel );
+ m_writer.setHighOrder( m_highOrder );
m_writer.write( time_n, cycleNumber, domain );
return false;
diff --git a/src/coreComponents/fileIO/Outputs/VTKOutput.hpp b/src/coreComponents/fileIO/Outputs/VTKOutput.hpp
index 93abd3e865d..215dc68a4fd 100644
--- a/src/coreComponents/fileIO/Outputs/VTKOutput.hpp
+++ b/src/coreComponents/fileIO/Outputs/VTKOutput.hpp
@@ -94,6 +94,8 @@ class VTKOutput : public OutputBase
static constexpr auto outputRegionTypeString = "outputRegionType";
static constexpr auto onlyPlotSpecifiedFieldNames = "onlyPlotSpecifiedFieldNames";
static constexpr auto fieldNames = "fieldNames";
+ static constexpr auto levelNames = "levelNames";
+ static constexpr auto highOrder = "highOrder";
} vtkOutputViewKeys;
/// @endcond
@@ -110,6 +112,7 @@ class VTKOutput : public OutputBase
string m_plotFileRoot;
integer m_writeFaceMesh;
integer m_plotLevel;
+ integer m_highOrder;
/// Should the vtk files contain the ghost cells or not.
integer m_writeGhostCells;
@@ -120,6 +123,9 @@ class VTKOutput : public OutputBase
/// array of names of the fields to output
array1d< string > m_fieldNames;
+ /// array of names of the mesh levels to output
+ array1d< string > m_levelNames;
+
/// VTK output mode
vtk::VTKOutputMode m_writeBinaryData = vtk::VTKOutputMode::BINARY;
diff --git a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp
index d01bf903b57..2672ea45983 100644
--- a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp
+++ b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.cpp
@@ -50,6 +50,7 @@ VTKPolyDataWriterInterface::VTKPolyDataWriterInterface( string name ):
m_pvd( m_outputName + ".pvd" ),
m_writeGhostCells( false ),
m_plotLevel( PlotLevel::LEVEL_1 ),
+ m_highOrder( true ),
m_requireFieldRegistrationCheck( true ),
m_previousCycle( -1 ),
m_outputMode( VTKOutputMode::BINARY ),
@@ -57,7 +58,7 @@ VTKPolyDataWriterInterface::VTKPolyDataWriterInterface( string name ):
{}
static int
-toVTKCellType( ElementType const elementType )
+toVTKCellType( ElementType const elementType, localIndex numNodes )
{
switch( elementType )
{
@@ -69,7 +70,16 @@ toVTKCellType( ElementType const elementType )
case ElementType::Tetrahedron: return VTK_TETRA;
case ElementType::Pyramid: return VTK_PYRAMID;
case ElementType::Wedge: return VTK_WEDGE;
- case ElementType::Hexahedron: return VTK_HEXAHEDRON;
+ case ElementType::Hexahedron:
+ switch( numNodes )
+ {
+ case 8:
+ return VTK_HEXAHEDRON;
+ case 27:
+ return VTK_QUADRATIC_HEXAHEDRON;
+ default:
+ return VTK_HEXAHEDRON;
+ }
case ElementType::Prism5: return VTK_PENTAGONAL_PRISM;
case ElementType::Prism6: return VTK_HEXAGONAL_PRISM;
case ElementType::Prism7: return VTK_POLYHEDRON;
@@ -100,7 +110,7 @@ toVTKCellType( ElementType const elementType )
* nodes for mapping purpose while keeping a face streams data structure. The negative values are
* converted to positives when generating the VTK_POLYHEDRON. Check getVtkCells() for more details.
*/
-static std::vector< int > getVtkConnectivity( ElementType const elementType )
+static std::vector< int > getVtkConnectivity( ElementType const elementType, localIndex const numNodes )
{
switch( elementType )
{
@@ -112,7 +122,20 @@ static std::vector< int > getVtkConnectivity( ElementType const elementType )
case ElementType::Tetrahedron: return { 0, 1, 2, 3 };
case ElementType::Pyramid: return { 0, 1, 3, 2, 4 };
case ElementType::Wedge: return { 0, 4, 2, 1, 5, 3 };
- case ElementType::Hexahedron: return { 0, 1, 3, 2, 4, 5, 7, 6 };
+ case ElementType::Hexahedron:
+ switch( numNodes )
+ {
+ case 8:
+ return { 0, 1, 3, 2, 4, 5, 7, 6 };
+ break;
+ case 27:
+ // WARNING: numbering convention changed between VTK writer 1.0 and 2.2
+ // (https://discourse.julialang.org/t/writevtk-node-numbering-for-27-node-lagrange-hexahedron/93698/8) Geos uses 1.0
+ return { 0, 2, 8, 6, 18, 20, 26, 24, 1, 5, 7, 3, 19, 23, 25, 21, 9, 11, 17, 15, 12, 14, 10, 16, 4, 22, 17 };
+ break;
+ default:
+ return { }; // TODO
+ }
case ElementType::Prism5: return { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
case ElementType::Prism6: return { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
case ElementType::Prism7: return {-9,
@@ -276,7 +299,8 @@ getSurface( FaceElementSubRegion const & subRegion,
geosx2VTKIndexing.reserve( subRegion.size() * subRegion.numNodesPerElement() );
localIndex nodeIndexInVTK = 0;
std::vector< vtkIdType > connectivity( subRegion.numNodesPerElement() );
- std::vector< int > const vtkOrdering = getVtkConnectivity( subRegion.getElementType() );
+ auto numNodes = nodeListPerElement.toViewConst().size();
+ std::vector< int > const vtkOrdering = getVtkConnectivity( subRegion.getElementType(), numNodes );
for( localIndex ei = 0; ei < subRegion.size(); ei++ )
{
@@ -308,10 +332,11 @@ getSurface( FaceElementSubRegion const & subRegion,
VTKCellType const type = [&]()
{
- switch( subRegion.numNodesPerElement() )
+ switch( numNodes )
{
case 6: return VTK_WEDGE;
case 8: return VTK_HEXAHEDRON;
+ case 27: return VTK_QUADRATIC_HEXAHEDRON;
default:
{
GEOS_ERROR( GEOS_FMT( "Elements with {} nodes can't be output in the subregion {}",
@@ -426,7 +451,8 @@ getVtkCells( CellElementRegion const & region,
localIndex numConn = 0;
region.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion const & subRegion )
{
- numConn += subRegion.size() * getVtkConnectivity( subRegion.getElementType() ).size();
+ auto const nodeList = subRegion.nodeList().toViewConst();
+ numConn += subRegion.size() * getVtkConnectivity( subRegion.getElementType(), nodeList.size( 1 ) ).size();
} );
return numConn;
}();
@@ -445,8 +471,9 @@ getVtkCells( CellElementRegion const & region,
localIndex connOffset = 0;
region.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion const & subRegion )
{
- cellTypes.insert( cellTypes.end(), subRegion.size(), toVTKCellType( subRegion.getElementType() ) );
- std::vector< int > const vtkOrdering = getVtkConnectivity( subRegion.getElementType() );
+ auto const numNodes = subRegion.nodeList().toViewConst().size( 1 );
+ cellTypes.insert( cellTypes.end(), subRegion.size(), toVTKCellType( subRegion.getElementType(), numNodes ) );
+ std::vector< int > const vtkOrdering = getVtkConnectivity( subRegion.getElementType(), numNodes );
localIndex const numVtkData = vtkOrdering.size();
auto const nodeList = subRegion.nodeList().toViewConst();
@@ -956,19 +983,25 @@ void VTKPolyDataWriterInterface::writeVtmFile( integer const cycle,
{
meshBody.forMeshLevels( [&]( MeshLevel const & meshLevel )
{
-
if( meshLevel.isShallowCopy() )
- {
return;
+
+ string const & meshLevelName = meshLevel.getName();
+ string const & meshBodyName = meshBody.getName();
+
+ if( !m_levelNames.empty())
+ {
+ if( m_levelNames.find( meshLevelName ) == m_levelNames.end())
+ return;
}
ElementRegionManager const & elemManager = meshLevel.getElemManager();
- string const meshPath = joinPath( getCycleSubFolder( cycle ), meshBody.getName(), meshLevel.getName() );
+ string const meshPath = joinPath( getCycleSubFolder( cycle ), meshBodyName, meshLevelName );
int const mpiSize = MpiWrapper::commSize();
auto addRegion = [&]( ElementRegionBase const & region )
{
- std::vector< string > const blockPath{ meshBody.getName(), meshLevel.getName(), region.getCatalogName(), region.getName() };
+ std::vector< string > const blockPath{ meshBodyName, meshLevelName, region.getCatalogName(), region.getName() };
string const regionPath = joinPath( meshPath, region.getName() );
for( int i = 0; i < mpiSize; i++ )
{
@@ -1052,11 +1085,9 @@ void VTKPolyDataWriterInterface::write( real64 const time,
integer const cycle,
DomainPartition const & domain )
{
- // This guard prevents crashes observed on MacOS due to a floating point exception
+ // This guard prevents crashes observed due to a floating point exception
// triggered inside VTK by a progress indicator
-#if defined(__APPLE__) && defined(__MACH__)
LvArray::system::FloatingPointExceptionGuard guard;
-#endif
string const stepSubDir = joinPath( m_outputName, getCycleSubFolder( cycle ) );
string const stepSubDirFull = joinPath( m_outputDir, stepSubDir );
@@ -1073,11 +1104,8 @@ void VTKPolyDataWriterInterface::write( real64 const time,
{
meshBody.forMeshLevels( [&]( MeshLevel const & meshLevel )
{
-
if( meshLevel.isShallowCopy() )
- {
return;
- }
ElementRegionManager const & elemManager = meshLevel.getElemManager();
NodeManager const & nodeManager = meshLevel.getNodeManager();
@@ -1085,6 +1113,12 @@ void VTKPolyDataWriterInterface::write( real64 const time,
string const & meshLevelName = meshLevel.getName();
string const & meshBodyName = meshBody.getName();
+ if( !m_levelNames.empty())
+ {
+ if( m_levelNames.find( meshLevelName ) == m_levelNames.end())
+ return;
+ }
+
if( m_requireFieldRegistrationCheck && !m_fieldNames.empty() )
{
outputUtilities::checkFieldRegistration( elemManager,
diff --git a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.hpp b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.hpp
index 6990cfec159..1ea250aafd5 100644
--- a/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.hpp
+++ b/src/coreComponents/fileIO/vtk/VTKPolyDataWriterInterface.hpp
@@ -96,6 +96,15 @@ class VTKPolyDataWriterInterface
m_plotLevel = dataRepository::toPlotLevel( plotLevel );
}
+ /**
+ * @brief Defines if the vtk outputs should the vtk files contain the high order solution.
+ * @param highOrder The boolean flag.
+ */
+ void setHighOrder( bool highOrder )
+ {
+ m_highOrder = highOrder;
+ }
+
/**
* @brief Set the binary mode
* @param[in] mode output mode to be set
@@ -144,6 +153,14 @@ class VTKPolyDataWriterInterface
m_fieldNames.insert( fieldNames.begin(), fieldNames.end() );
}
+ /**
+ * @brief Set the names of the mesh levels to output
+ * @param[in] levelNames the mesh levels to output
+ */
+ void setLevelNames( arrayView1d< string const > const & levelNames )
+ {
+ m_levelNames.insert( levelNames.begin(), levelNames.end() );
+ }
/**
* @brief Main method of this class. Write all the files for one time step.
@@ -298,12 +315,18 @@ class VTKPolyDataWriterInterface
/// Maximum plot level to be written.
dataRepository::PlotLevel m_plotLevel;
+ /// Should the vtk files contain the high order solution or not.
+ bool m_highOrder;
+
/// Flag to decide whether we only plot the fields specified by fieldNames, or if we also plot fields based on plotLevel
integer m_onlyPlotSpecifiedFieldNames;
/// Flag to decide whether we check that the specified fieldNames are actually registered
bool m_requireFieldRegistrationCheck;
+ /// Names of the mesh levels to output
+ std::set< string > m_levelNames;
+
/// Names of the fields to output
std::set< string > m_fieldNames;
diff --git a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis2.hpp b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis2.hpp
index 62e267b30cb..fa84fe6d6ab 100644
--- a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis2.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis2.hpp
@@ -172,7 +172,7 @@ class LagrangeBasis2
}
/**
- * @brief The gradient of the basis function for support point 1 evaluated at
+ * @brief The gradient of the basis function for support point 2 evaluated at
* a point along the axes.
* @param xi The coordinate at which to evaluate the gradient.
* @return The gradient of basis function
diff --git a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis3GL.hpp b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis3GL.hpp
index 1e07db8f13d..8a290894350 100644
--- a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis3GL.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis3GL.hpp
@@ -250,7 +250,7 @@ class LagrangeBasis3GL
}
/**
- * @brief The gradient of the basis function for support point 1 evaluated at
+ * @brief The gradient of the basis function for support point 2 evaluated at
* a point along the axes.
* @param xi The coordinate at which to evaluate the gradient.
* @return The gradient of basis function
diff --git a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis4GL.hpp b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis4GL.hpp
index 0b3cb0f0742..4f357b418e8 100644
--- a/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis4GL.hpp
+++ b/src/coreComponents/finiteElement/elementFormulations/LagrangeBasis4GL.hpp
@@ -271,7 +271,7 @@ class LagrangeBasis4GL
}
/**
- * @brief The gradient of the basis function for support point 1 evaluated at
+ * @brief The gradient of the basis function for support point 2 evaluated at
* a point along the axes.
* @param xi The coordinate at which to evaluate the gradient.
* @return The gradient of basis function
diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt
index e834480cd08..a302f8d85d3 100644
--- a/src/coreComponents/physicsSolvers/CMakeLists.txt
+++ b/src/coreComponents/physicsSolvers/CMakeLists.txt
@@ -135,6 +135,7 @@ set( physicsSolvers_headers
wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp
wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp
wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp
+ wavePropagation/AcousticElasticWaveEquationSEM.hpp
)
#
@@ -198,7 +199,7 @@ set( physicsSolvers_sources
wavePropagation/ElasticWaveEquationSEM.cpp
wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp
wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp
-
+ wavePropagation/AcousticElasticWaveEquationSEM.cpp
)
include( solidMechanics/kernels/SolidMechanicsKernels.cmake)
diff --git a/src/coreComponents/physicsSolvers/SolverBase.hpp b/src/coreComponents/physicsSolvers/SolverBase.hpp
index 8b649643546..87cc5aa7dbd 100644
--- a/src/coreComponents/physicsSolvers/SolverBase.hpp
+++ b/src/coreComponents/physicsSolvers/SolverBase.hpp
@@ -673,10 +673,22 @@ class SolverBase : public ExecutableGroup
MeshBody const & meshBody = meshBodies.getGroup< MeshBody >( meshBodyName );
MeshLevel const * meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( meshLevelName );
- if( meshLevelPtr==nullptr )
+ bool nok = meshLevelPtr==nullptr;
+ if( nok )
{
meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( MeshBody::groupStructKeys::baseDiscretizationString() );
}
+ /*
+ std::cout << "\t[SolverBase::forDiscretizationOnMeshTargets] ok=" << (nok ? 'F' : 'T')
+ << " meshBodyName=" << meshBodyName
+ << " meshLevelName=" << meshLevelName
+ << " discretizationName=" << m_discretizationName
+ << std::endl;
+ for (auto const & val: regionNames)
+ {
+ std::cout << "\t[SolverBase::forDiscretizationOnMeshTargets] regionNames=" << val << std::endl;
+ }
+ */
lambda( meshBodyName, *meshLevelPtr, regionNames );
}
}
@@ -699,10 +711,22 @@ class SolverBase : public ExecutableGroup
MeshBody & meshBody = meshBodies.getGroup< MeshBody >( meshBodyName );
MeshLevel * meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( meshLevelName );
- if( meshLevelPtr==nullptr )
+ bool nok = meshLevelPtr==nullptr;
+ if( nok )
{
meshLevelPtr = meshBody.getMeshLevels().getGroupPointer< MeshLevel >( MeshBody::groupStructKeys::baseDiscretizationString() );
}
+ /*
+ std::cout << "\t[SolverBase::forDiscretizationOnMeshTargets] ok=" << (nok ? 'F' : 'T')
+ << " meshBodyName=" << meshBodyName
+ << " meshLevelName=" << meshLevelName
+ << " discretizationName=" << m_discretizationName
+ << std::endl;
+ for (auto const & val: regionNames)
+ {
+ std::cout << "\t[SolverBase::forDiscretizationOnMeshTargets] regionNames=" << val << std::endl;
+ }
+ */
lambda( meshBodyName, *meshLevelPtr, regionNames );
}
}
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEM.cpp
new file mode 100644
index 00000000000..8b3482f7cb2
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEM.cpp
@@ -0,0 +1,221 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019- GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file AcousticElasticWaveEquationSEM.cpp
+ */
+
+#include "AcousticElasticWaveEquationSEM.hpp"
+#include "AcousticElasticWaveEquationSEMKernel.hpp"
+#include "dataRepository/Group.hpp"
+#include
+#include
+
+namespace geos
+{
+using namespace dataRepository;
+
+void AcousticElasticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
+{
+ SolverBase::registerDataOnMesh( meshBodies );
+
+ forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & )
+ {
+ NodeManager & nodeManager = mesh.getNodeManager();
+ nodeManager.registerField< fields::CouplingVectorx >( getName() );
+ nodeManager.registerField< fields::CouplingVectory >( getName() );
+ nodeManager.registerField< fields::CouplingVectorz >( getName() );
+ } );
+}
+
+void AcousticElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
+{
+ SolverBase::initializePostInitialConditionsPreSubGroups();
+
+ auto acousSolver = acousticSolver();
+ auto elasSolver = elasticSolver();
+
+ auto acousNodesSet = acousSolver->getSolverNodesSet();
+ auto elasNodesSet = elasSolver->getSolverNodesSet();
+
+ for( auto val : acousNodesSet )
+ {
+ if( elasNodesSet.contains( val ) )
+ m_interfaceNodesSet.insert( val );
+ }
+ localIndex const numInterfaceNodes = MpiWrapper::sum( m_interfaceNodesSet.size() );
+ GEOS_THROW_IF( numInterfaceNodes == 0, "Failed to compute interface: check xml input (solver order)", std::runtime_error );
+
+ m_acousRegions = acousSolver->getReference< array1d< string > >( SolverBase::viewKeyStruct::targetRegionsString() );
+ m_elasRegions = elasSolver->getReference< array1d< string > >( SolverBase::viewKeyStruct::targetRegionsString() );
+
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & GEOS_UNUSED_PARAM( regionNames ) )
+ {
+ NodeManager & nodeManager = mesh.getNodeManager();
+ FaceManager & faceManager = mesh.getFaceManager();
+ ElementRegionManager & elemManager = mesh.getElemManager();
+
+ arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst();
+
+ arrayView2d< real64 const > const faceNormals = faceManager.faceNormal().toViewConst();
+ ArrayOfArraysView< localIndex const > const faceToNode = faceManager.nodeList().toViewConst();
+ arrayView2d< localIndex const > const faceToSubRegion = faceManager.elementSubRegionList();
+ arrayView2d< localIndex const > const faceToRegion = faceManager.elementRegionList();
+ arrayView2d< localIndex const > const faceToElement = faceManager.elementList();
+
+ arrayView1d< real32 > const couplingVectorx = nodeManager.getField< fields::CouplingVectorx >();
+ couplingVectorx.zero();
+
+ arrayView1d< real32 > const couplingVectory = nodeManager.getField< fields::CouplingVectory >();
+ couplingVectory.zero();
+
+ arrayView1d< real32 > const couplingVectorz = nodeManager.getField< fields::CouplingVectorz >();
+ couplingVectorz.zero();
+
+ elemManager.forElementRegions( m_acousRegions, [&] ( localIndex const regionIndex, ElementRegionBase const & elemRegion )
+ {
+ elemRegion.forElementSubRegionsIndex( [&]( localIndex const subRegionIndex, ElementSubRegionBase const & elementSubRegion )
+ {
+ finiteElement::FiniteElementBase const &
+ fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
+
+ finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
+ {
+ using FE_TYPE = TYPEOFREF( finiteElement );
+
+ acousticElasticWaveEquationSEMKernels::CouplingKernel< FE_TYPE > kernelC;
+ kernelC.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
+ nodeCoords,
+ regionIndex,
+ subRegionIndex,
+ faceToSubRegion,
+ faceToRegion,
+ faceToElement,
+ faceToNode,
+ faceNormals,
+ couplingVectorx,
+ couplingVectory,
+ couplingVectorz );
+ } );
+ } );
+ } );
+ } );
+}
+
+real64 AcousticElasticWaveEquationSEM::solverStep( real64 const & time_n,
+ real64 const & dt,
+ int const cycleNumber,
+ DomainPartition & domain )
+{
+ GEOS_MARK_FUNCTION;
+
+ auto acousSolver = acousticSolver();
+ auto elasSolver = elasticSolver();
+
+ SortedArrayView< localIndex const > const interfaceNodesSet = m_interfaceNodesSet.toViewConst();
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & GEOS_UNUSED_PARAM( regionNames ) )
+ {
+ NodeManager & nodeManager = mesh.getNodeManager();
+
+ arrayView1d< real32 const > const massA = nodeManager.getField< fields::MassVectorA >();
+ arrayView1d< real32 const > const massE = nodeManager.getField< fields::MassVectorE >();
+ arrayView1d< localIndex > const freeSurfaceNodeIndicatorA = nodeManager.getField< fields::FreeSurfaceNodeIndicatorA >();
+ arrayView1d< localIndex > const freeSurfaceNodeIndicatorE = nodeManager.getField< fields::FreeSurfaceNodeIndicatorE >();
+
+ arrayView1d< real32 const > const p_n = nodeManager.getField< fields::Pressure_n >();
+ arrayView1d< real32 const > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >();
+ arrayView1d< real32 const > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >();
+ arrayView1d< real32 const > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >();
+ arrayView1d< real32 const > const ux_n = nodeManager.getField< fields::Displacementx_n >();
+ arrayView1d< real32 const > const uy_n = nodeManager.getField< fields::Displacementy_n >();
+ arrayView1d< real32 const > const uz_n = nodeManager.getField< fields::Displacementz_n >();
+ arrayView1d< real32 const > const atoex = nodeManager.getField< fields::CouplingVectorx >();
+ arrayView1d< real32 const > const atoey = nodeManager.getField< fields::CouplingVectory >();
+ arrayView1d< real32 const > const atoez = nodeManager.getField< fields::CouplingVectorz >();
+
+ arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
+ arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
+ arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
+ arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
+
+ real32 const dt2 = pow( dt, 2 );
+
+ elasSolver->computeUnknowns( time_n, dt, cycleNumber, domain, mesh, m_elasRegions );
+
+ forAll< EXEC_POLICY >( interfaceNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
+ {
+ localIndex const a = interfaceNodesSet[n];
+ if( freeSurfaceNodeIndicatorE[a] == 1 )
+ return;
+
+ real32 const aux = -p_n[a] / massE[a];
+ real32 const localIncrementx = dt2 * atoex[a] * aux;
+ real32 const localIncrementy = dt2 * atoey[a] * aux;
+ real32 const localIncrementz = dt2 * atoez[a] * aux;
+
+ RAJA::atomicAdd< ATOMIC_POLICY >( &ux_np1[a], localIncrementx );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &uy_np1[a], localIncrementy );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &uz_np1[a], localIncrementz );
+ } );
+
+ elasSolver->synchronizeUnknowns( time_n, dt, cycleNumber, domain, mesh, m_elasRegions );
+
+ acousSolver->computeUnknowns( time_n, dt, cycleNumber, domain, mesh, m_acousRegions );
+
+ forAll< EXEC_POLICY >( interfaceNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
+ {
+ localIndex const a = interfaceNodesSet[n];
+ if( freeSurfaceNodeIndicatorA[a] == 1 )
+ return;
+
+ real32 const localIncrement = (
+ atoex[a] * ( ux_np1[a] - 2.0 * ux_n[a] + ux_nm1[a] ) +
+ atoey[a] * ( uy_np1[a] - 2.0 * uy_n[a] + uy_nm1[a] ) +
+ atoez[a] * ( uz_np1[a] - 2.0 * uz_n[a] + uz_nm1[a] )
+ ) / massA[a];
+
+ RAJA::atomicAdd< ATOMIC_POLICY >( &p_np1[a], localIncrement );
+ } );
+
+ acousSolver->synchronizeUnknowns( time_n, dt, cycleNumber, domain, mesh, m_acousRegions );
+
+ acousSolver->prepareNextTimestep( mesh );
+ elasSolver->prepareNextTimestep( mesh );
+ } );
+
+ return dt;
+}
+
+void AcousticElasticWaveEquationSEM::cleanup( real64 const time_n,
+ integer const cycleNumber,
+ integer const eventCounter,
+ real64 const eventProgress,
+ DomainPartition & domain )
+{
+ elasticSolver()->cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
+ acousticSolver()->cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain );
+}
+
+REGISTER_CATALOG_ENTRY( SolverBase, AcousticElasticWaveEquationSEM, string const &, Group * const )
+
+} /* namespace geos */
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEM.hpp
new file mode 100644
index 00000000000..63cf89f193c
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEM.hpp
@@ -0,0 +1,206 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019- GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+
+/**
+ * @file AcousticElasticWaveEquationSEM.hpp
+ */
+
+#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICELASTICWAVEEQUATIONSEM_HPP_
+#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICELASTICWAVEEQUATIONSEM_HPP_
+
+#include "physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp"
+#include "physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp"
+#include "physicsSolvers/SolverBase.hpp"
+#include
+
+namespace geos
+{
+
+template< typename ... SOLVERS >
+class CoupledWaveSolver : public SolverBase
+{
+
+public:
+
+ /**
+ * @brief main constructor for CoupledWaveSolver Objects
+ * @param name the name of this instantiation of CoupledWaveSolver in the repository
+ * @param parent the parent group of this instantiation of CoupledWaveSolver
+ */
+ CoupledWaveSolver( const string & name,
+ Group * const parent )
+ : SolverBase( name, parent )
+ {
+ forEachArgInTuple( m_solvers, [&]( auto solver, auto idx )
+ {
+ using SolverType = TYPEOFPTR( solver );
+ string const key = SolverType::coupledSolverAttributePrefix() + "SolverName";
+ registerWrapper( key, &m_names[idx()] ).
+ setInputFlag( dataRepository::InputFlags::REQUIRED ).
+ setDescription( "Name of the " + SolverType::coupledSolverAttributePrefix() + " solver used by the coupled solver" );
+ } );
+ }
+
+ /// deleted copy constructor
+ CoupledWaveSolver( CoupledWaveSolver const & ) = delete;
+
+ /// default move constructor
+ CoupledWaveSolver( CoupledWaveSolver && ) = default;
+
+ /// deleted assignment operator
+ CoupledWaveSolver & operator=( CoupledWaveSolver const & ) = delete;
+
+ /// deleted move operator
+ CoupledWaveSolver & operator=( CoupledWaveSolver && ) = delete;
+
+ virtual void
+ postProcessInput() override
+ {
+ SolverBase::postProcessInput();
+
+ forEachArgInTuple( m_solvers, [&]( auto & solver, auto idx )
+ {
+ using SolverPtr = TYPEOFREF( solver );
+ using SolverType = TYPEOFPTR( SolverPtr {} );
+ solver = getParent().template getGroupPointer< SolverType >( m_names[idx()] );
+ GEOS_THROW_IF( solver == nullptr,
+ GEOS_FMT( "Could not find solver '{}' of type {}",
+ m_names[idx()], LvArray::system::demangleType< SolverType >() ),
+ InputError );
+ } );
+ }
+
+protected:
+
+ /// Pointers of the single-physics solvers
+ std::tuple< SOLVERS *... > m_solvers;
+
+ /// Names of the single-physics solvers
+ std::array< string, sizeof...( SOLVERS ) > m_names;
+};
+
+
+class AcousticElasticWaveEquationSEM : public CoupledWaveSolver< AcousticWaveEquationSEM, ElasticWaveEquationSEM >
+{
+public:
+ using Base = CoupledWaveSolver< AcousticWaveEquationSEM, ElasticWaveEquationSEM >;
+ using Base::m_solvers;
+ using wsCoordType = AcousticWaveEquationSEM::wsCoordType;
+
+ enum class SolverType : integer
+ {
+ AcousticWaveEquationSEM = 0,
+ ElasticWaveEquationSEM = 1
+ };
+
+ /// String used to form the solverName used to register solvers in CoupledWaveSolver
+ static string coupledSolverAttributePrefix() { return "acousticelastic"; }
+
+ using EXEC_POLICY = parallelDevicePolicy< >;
+ using ATOMIC_POLICY = AtomicPolicy< EXEC_POLICY >;
+
+ virtual void registerDataOnMesh( Group & meshBodies ) override final;
+
+ /**
+ * @brief main constructor for AcousticElasticWaveEquationSEM objects
+ * @param name the name of this instantiation of AcousticElasticWaveEquationSEM in the repository
+ * @param parent the parent group of this instantiation of AcousticElasticWaveEquationSEM
+ */
+ AcousticElasticWaveEquationSEM( const string & name,
+ Group * const parent )
+ : Base( name, parent )
+ { }
+
+ /// Destructor for the class
+ ~AcousticElasticWaveEquationSEM() override {}
+
+ /**
+ * @brief name of the node manager in the object catalog
+ * @return string that contains the catalog name to generate a new AcousticElasticWaveEquationSEM object through the object catalog.
+ */
+ static string catalogName() { return "AcousticElasticSEM"; }
+
+ /**
+ * @brief accessor for the pointer to the solid mechanics solver
+ * @return a pointer to the solid mechanics solver
+ */
+ AcousticWaveEquationSEM * acousticSolver() const
+ {
+ return std::get< toUnderlying( SolverType::AcousticWaveEquationSEM ) >( m_solvers );
+ }
+
+ /**
+ * @brief accessor for the pointer to the flow solver
+ * @return a pointer to the flow solver
+ */
+ ElasticWaveEquationSEM * elasticSolver() const
+ {
+ return std::get< toUnderlying( SolverType::ElasticWaveEquationSEM ) >( m_solvers );
+ }
+
+ // (requires not to be private because it is called from GEOS_HOST_DEVICE method)
+ virtual real64
+ solverStep( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain ) override;
+
+ virtual void
+ cleanup( real64 const time_n,
+ integer const cycleNumber,
+ integer const eventCounter,
+ real64 const eventProgress,
+ DomainPartition & domain ) override;
+
+protected:
+
+ virtual void initializePostInitialConditionsPreSubGroups() override;
+
+ SortedArray< localIndex > m_interfaceNodesSet;
+ arrayView1d< string const > m_acousRegions;
+ arrayView1d< string const > m_elasRegions;
+};
+
+namespace fields
+{
+
+DECLARE_FIELD( CouplingVectorx,
+ "couplingVectorx",
+ array1d< real32 >,
+ 0,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Coupling term on x." );
+
+DECLARE_FIELD( CouplingVectory,
+ "couplingVectory",
+ array1d< real32 >,
+ 0,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Coupling term on y." );
+
+DECLARE_FIELD( CouplingVectorz,
+ "couplingVectorz",
+ array1d< real32 >,
+ 0,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Coupling term on z." );
+}
+
+} /* namespace geos */
+
+#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICELASTICWAVEEQUATIONSEM_HPP_ */
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEMKernel.hpp
new file mode 100644
index 00000000000..aa3a56b1c0b
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticElasticWaveEquationSEMKernel.hpp
@@ -0,0 +1,102 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
+ * Copyright (c) 2018-2020 TotalEnergies
+ * Copyright (c) 2019- GEOSX Contributors
+ * All rights reserved
+ *
+ * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+/**
+ * @file AcousticElasticWaveEquationSEMKernel.hpp
+ */
+
+#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICELASTICWAVEEQUATIONSEMKERNEL_HPP_
+#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICELASTICWAVEEQUATIONSEMKERNEL_HPP_
+
+#include "finiteElement/kernelInterface/KernelBase.hpp"
+#if !defined( GEOS_USE_HIP )
+#include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp"
+#endif
+
+#include
+
+namespace geos
+{
+
+namespace acousticElasticWaveEquationSEMKernels
+{
+
+template< typename FE_TYPE >
+struct CouplingKernel
+{
+ static constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
+
+ template< typename EXEC_POLICY, typename ATOMIC_POLICY >
+ void
+ launch( localIndex const size,
+ arrayView2d< WaveSolverBase::wsCoordType const,
+ nodes::REFERENCE_POSITION_USD > const nodeCoords,
+ localIndex const regionIndex,
+ localIndex const subRegionIndex,
+ arrayView2d< localIndex const > const faceToSubRegion,
+ arrayView2d< localIndex const > const faceToRegion,
+ arrayView2d< localIndex const > const faceToElement,
+ ArrayOfArraysView< localIndex const > const facesToNodes,
+ arrayView2d< real64 const > const faceNormals,
+ arrayView1d< real32 > const couplingVectorx,
+ arrayView1d< real32 > const couplingVectory,
+ arrayView1d< real32 > const couplingVectorz )
+ {
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
+ {
+ localIndex e0 = faceToElement( f, 0 ), e1 = faceToElement( f, 1 );
+ localIndex er0 = faceToRegion( f, 0 ), er1 = faceToRegion( f, 1 );
+ localIndex esr0 = faceToSubRegion( f, 0 ), esr1 = faceToSubRegion( f, 1 );
+
+ if( e0 != -1 && e1 != -1 && er0 != er1 ) // an interface is defined as a transition between regions
+ {
+ // check that one of the region is the fluid subregion for the fluid -> solid coupling term
+ if((er0 == regionIndex && esr0 == subRegionIndex) || (er1 == regionIndex && esr1 == subRegionIndex))
+ {
+ real64 xLocal[ numNodesPerFace ][ 3 ];
+ for( localIndex a = 0; a < numNodesPerFace; ++a )
+ {
+ for( localIndex i = 0; i < 3; ++i )
+ {
+ xLocal[a][i] = nodeCoords( facesToNodes( f, a ), i );
+ }
+ }
+
+ // determine normal sign for fluid -> solid coupling
+ localIndex sgn = er0 == regionIndex ? 1 : (er1 == regionIndex ? -1 : 0);
+
+ for( localIndex q = 0; q < numNodesPerFace; ++q )
+ {
+ real64 const aux = FE_TYPE::computeDampingTerm( q, xLocal );
+
+ real32 const localIncrementx = aux * (sgn * faceNormals( f, 0 ));
+ real32 const localIncrementy = aux * (sgn * faceNormals( f, 1 ));
+ real32 const localIncrementz = aux * (sgn * faceNormals( f, 2 ));
+
+ RAJA::atomicAdd< ATOMIC_POLICY >( &couplingVectorx[facesToNodes( f, q )], localIncrementx );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &couplingVectory[facesToNodes( f, q )], localIncrementy );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &couplingVectorz[facesToNodes( f, q )], localIncrementz );
+ }
+ }
+ }
+ } );
+
+ }
+};
+
+} /* namespace acousticElasticWaveEquationSEMKernels */
+
+} /* namespace geos */
+
+#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICELASTICWAVEEQUATIONSEMKERNEL_HPP_ */
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp
index ad01f33626a..40f83983dfb 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.cpp
@@ -65,21 +65,11 @@ AcousticFirstOrderWaveEquationSEM::AcousticFirstOrderWaveEquationSEM( const std:
setSizedFromParent( 0 ).
setDescription( "Element containing the sources" );
- registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ).
- setInputFlag( InputFlags::FALSE ).
- setSizedFromParent( 0 ).
- setDescription( "Element containing the receivers" );
-
registerWrapper( viewKeyStruct::sourceRegionString(), &m_sourceRegion ).
setInputFlag( InputFlags::FALSE ).
setSizedFromParent( 0 ).
setDescription( "Region containing the sources" );
- registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ).
- setInputFlag( InputFlags::FALSE ).
- setSizedFromParent( 0 ).
- setDescription( "Region containing the receivers" );
-
}
AcousticFirstOrderWaveEquationSEM::~AcousticFirstOrderWaveEquationSEM()
@@ -109,21 +99,21 @@ void AcousticFirstOrderWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
wavesolverfields::ForcingRHS,
wavesolverfields::MassVector,
wavesolverfields::DampingVector,
- wavesolverfields::FreeSurfaceNodeIndicator >( this->getName() );
+ wavesolverfields::FreeSurfaceNodeIndicator >( getName() );
FaceManager & faceManager = mesh.getFaceManager();
- faceManager.registerField< wavesolverfields::FreeSurfaceFaceIndicator >( this->getName() );
+ faceManager.registerField< wavesolverfields::FreeSurfaceFaceIndicator >( getName() );
ElementRegionManager & elemManager = mesh.getElemManager();
elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion )
{
- subRegion.registerField< wavesolverfields::MediumVelocity >( this->getName() );
- subRegion.registerField< wavesolverfields::MediumDensity >( this->getName() );
+ subRegion.registerField< wavesolverfields::MediumVelocity >( getName() );
+ subRegion.registerField< wavesolverfields::MediumDensity >( getName() );
- subRegion.registerField< wavesolverfields::Velocity_x >( this->getName() );
- subRegion.registerField< wavesolverfields::Velocity_y >( this->getName() );
- subRegion.registerField< wavesolverfields::Velocity_z >( this->getName() );
+ subRegion.registerField< wavesolverfields::Velocity_x >( getName() );
+ subRegion.registerField< wavesolverfields::Velocity_y >( getName() );
+ subRegion.registerField< wavesolverfields::Velocity_z >( getName() );
finiteElement::FiniteElementBase const & fe = subRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
@@ -151,10 +141,10 @@ void AcousticFirstOrderWaveEquationSEM::postProcessInput()
localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 );
localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 );
- m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_uxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
+ m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_uxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_uyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_uzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
m_sourceElem.resize( numSourcesGlobal );
m_sourceRegion.resize( numSourcesGlobal );
m_rcvElem.resize( numReceiversGlobal );
@@ -191,15 +181,13 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
receiverConstants.setValues< EXEC_POLICY >( -1 );
receiverIsLocal.zero();
- real32 const timeSourceFrequency = this->m_timeSourceFrequency;
- localIndex const rickerOrder = this->m_rickerOrder;
arrayView2d< real32 > const sourceValue = m_sourceValue.toView();
real64 dt = 0;
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
{
EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
- if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+ if( subEvent->getEventName() == "/Solvers/" + getName() )
{
dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
}
@@ -254,8 +242,9 @@ void AcousticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLev
receiverRegion,
sourceValue,
dt,
- timeSourceFrequency,
- rickerOrder );
+ m_timeSourceFrequency,
+ m_timeSourceDelay,
+ m_rickerOrder );
} );
} );
@@ -287,10 +276,9 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
{
WaveSolverBase::initializePostInitialConditionsPreSubGroups();
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
- real64 const time = 0.0;
- applyFreeSurfaceBC( time, domain );
+ applyFreeSurfaceBC( 0.0, domain );
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
@@ -323,7 +311,7 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
CellElementSubRegion & elementSubRegion )
{
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
- arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
+ arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();
arrayView1d< real32 const > const velocity = elementSubRegion.getField< wavesolverfields::MediumVelocity >();
arrayView1d< real32 const > const density = elementSubRegion.getField< wavesolverfields::MediumDensity >();
@@ -344,9 +332,9 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
acousticFirstOrderWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );
- kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
+ kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
X,
- facesToElements,
+ elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
@@ -355,6 +343,8 @@ void AcousticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGro
} );
} );
} );
+
+ WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
}
@@ -525,9 +515,9 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
cycleNumber,
p_np1 );
} );
- arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView();
- arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView();
- arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView();
+ arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView();
+ arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView();
+ arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView();
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_x, velocity_x, uxReceivers );
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, velocity_y, velocity_y, uyReceivers );
@@ -546,16 +536,12 @@ real64 AcousticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & t
true );
// compute the seismic traces since last step.
- arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
+ arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
computeAllSeismoTraces( time_n, dt, p_np1, p_np1, pReceivers );
+ } );
- // increment m_indexSeismoTrace
- while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
- {
- m_indexSeismoTrace++;
- }
+ incrementIndexSeismoTrace( time_n );
- } );
return dt;
}
@@ -575,62 +561,26 @@ void AcousticFirstOrderWaveEquationSEM::cleanup( real64 const time_n, integer co
arrayView2d< real32 > const velocity_y = elementSubRegion.getField< wavesolverfields::Velocity_y >();
arrayView2d< real32 > const velocity_z = elementSubRegion.getField< wavesolverfields::Velocity_z >();
- arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView();
- arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView();
- arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView();
+ arrayView2d< real32 > const uxReceivers = m_uxNp1AtReceivers.toView();
+ arrayView2d< real32 > const uyReceivers = m_uyNp1AtReceivers.toView();
+ arrayView2d< real32 > const uzReceivers = m_uzNp1AtReceivers.toView();
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_x, velocity_x, uxReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_y, velocity_y, uyReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, velocity_z, velocity_z, uzReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_x, velocity_x, uxReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_y, velocity_y, uyReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, velocity_z, velocity_z, uzReceivers );
- } );
- arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
- computeAllSeismoTraces( time_n, 0, p_np1, p_np1, pReceivers );
- } );
+ WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, uxReceivers, uyReceivers, uzReceivers );
- // increment m_indexSeismoTrace
- while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
- {
- m_indexSeismoTrace++;
- }
-}
-
-void AcousticFirstOrderWaveEquationSEM::computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers )
-{
- localIndex indexSeismoTrace = m_indexSeismoTrace;
- for( real64 timeSeismo;
- (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
- indexSeismoTrace++ )
- {
- WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace,
- m_outputSeismoTrace,
- var_np1, var_n, varAtReceivers );
- }
-}
+ } );
+ arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
+ computeAllSeismoTraces( time_n, 0.0, p_np1, p_np1, pReceivers );
+ WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers );
-void AcousticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( localIndex const regionIndex,
- real64 const time_n,
- real64 const dt,
- arrayView2d< real32 const > const var_np1,
- arrayView2d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers )
-{
- localIndex indexSeismoTrace = m_indexSeismoTrace;
- for( real64 timeSeismo;
- (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
- indexSeismoTrace++ )
- {
- WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dt, regionIndex, m_receiverRegion, timeSeismo, indexSeismoTrace, m_rcvElem, m_receiverConstants, m_receiverIsLocal, m_nsamplesSeismoTrace,
- m_outputSeismoTrace,
- var_np1, var_n, varAtReceivers );
- }
+ } );
}
-
void AcousticFirstOrderWaveEquationSEM::initializePML()
{
GEOS_ERROR( "PML for the first order acoustic wave propagator not yet implemented" );
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp
index 97bcfe470ce..d534765d3f7 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEM.hpp
@@ -21,8 +21,8 @@
#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_ACOUSTICFIRSTORDERWAVEEQUATIONSEM_HPP_
#include "mesh/MeshFields.hpp"
-#include "WaveSolverUtils.hpp"
#include "WaveSolverBaseFields.hpp"
+#include "WaveSolverBase.hpp"
namespace geos
{
@@ -78,38 +78,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
*/
virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs );
- /**
- * TODO: move implementation into WaveSolverUtils
- * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
- * @param time_n the time corresponding to the field values pressure_n
- * @param dt the simulation timestep
- * @param var_np1 the field values at time_n + dt
- * @param var_n the field values at time_n
- * @param var_receivers the array holding the trace values, where the output is written
- */
- virtual void computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers );
-
- /**
- * TODO: move implementation into WaveSolverUtils
- * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt for a 2d variable
- * @param time_n the time corresponding to the field values pressure_n
- * @param dt the simulation timestep
- * @param var_np1 the field values at time_n + dt
- * @param var_n the field values at time_n
- * @param var_receivers the array holding the trace values, where the output is written
- */
- virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex,
- real64 const time_n,
- real64 const dt,
- arrayView2d< real32 const > const var_np1,
- arrayView2d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers );
-
-
/**
* @brief Initialize Perfectly Matched Layer (PML) information
*/
@@ -133,8 +101,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
static constexpr char const * sourceElemString() { return "sourceElem"; }
static constexpr char const * sourceRegionString() { return "sourceRegion"; }
- static constexpr char const * receiverElemString() { return "rcvElem"; }
- static constexpr char const * receiverRegionString() { return "receiverRegion"; }
} waveEquationViewKeys;
@@ -198,12 +164,6 @@ class AcousticFirstOrderWaveEquationSEM : public WaveSolverBase
/// Array containing the elements which contain the region which the source belongs
array1d< localIndex > m_sourceRegion;
- /// Array containing the elements which contain the region which the receiver belongs
- array1d< localIndex > m_receiverRegion;
-
- /// Array containing the elements which contain a receiver
- array1d< localIndex > m_rcvElem;
-
};
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp
index 712e909fdb5..3c0d6ab1c1d 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticFirstOrderWaveEquationSEMKernel.hpp
@@ -61,6 +61,7 @@ struct PrecomputeSourceAndReceiverKernel
* @param[out] sourceValue value of the temporal source (eg. Ricker)
* @param[in] dt time-step
* @param[in] timeSourceFrequency the central frequency of the source
+ * @param[in] timeSourceDelay the time delay of the source
* @param[in] rickerOrder order of the Ricker wavelet
*/
template< typename EXEC_POLICY, typename FE_TYPE >
@@ -91,6 +92,7 @@ struct PrecomputeSourceAndReceiverKernel
arrayView2d< real32 > const sourceValue,
real64 const dt,
real32 const timeSourceFrequency,
+ real32 const timeSourceDelay,
localIndex const rickerOrder )
{
@@ -136,14 +138,13 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- sourceNodeIds[isrc][a] = elemsToNodes[k][a];
+ sourceNodeIds[isrc][a] = elemsToNodes( k, a );
sourceConstants[isrc][a] = Ntest[a];
}
for( localIndex cycle = 0; cycle < sourceValue.size( 0 ); ++cycle )
{
- real64 const time = cycle*dt;
- sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( time, timeSourceFrequency, rickerOrder );
+ sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( cycle * dt, timeSourceFrequency, timeSourceDelay, rickerOrder );
}
}
}
@@ -185,7 +186,7 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- receiverNodeIds[ircv][a] = elemsToNodes[k][a];
+ receiverNodeIds[ircv][a] = elemsToNodes( k, a );
receiverConstants[ircv][a] = Ntest[a];
}
}
@@ -245,7 +246,7 @@ struct MassMatrixKernel
for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
{
real32 const localIncrement = invC2 * m_finiteElement.computeMassTerm( q, xLocal );
- RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( k, q )], localIncrement );
}
} ); // end loop over element
}
@@ -268,8 +269,8 @@ struct DampingMatrixKernel
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
- * @param[in] X coordinates of the nodes
- * @param[in] facesToElems map from faces to elements
+ * @param[in] nodeCoords coordinates of the nodes
+ * @param[in] elemsToFaces map from elements to faces
* @param[in] facesToNodes map from face to nodes
* @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise
* @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise
@@ -279,42 +280,40 @@ struct DampingMatrixKernel
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
void
launch( localIndex const size,
- arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
- arrayView2d< localIndex const > const facesToElems,
+ arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
+ arrayView2d< localIndex const > const elemsToFaces,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView1d< integer const > const facesDomainBoundaryIndicator,
arrayView1d< localIndex const > const freeSurfaceFaceIndicator,
arrayView1d< real32 const > const velocity,
arrayView1d< real32 > const damping )
{
- forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
- // face on the domain boundary and not on free surface
- if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
+ for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i )
{
- localIndex k = facesToElems( f, 0 );
- if( k < 0 )
+ localIndex const f = elemsToFaces( e, i );
+ // face on the domain boundary and not on free surface
+ if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
{
- k = facesToElems( f, 1 );
- }
-
- real32 const alpha = 1.0 / velocity[k];
+ real32 const alpha = 1.0 / velocity[e];
- constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
+ constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
- real64 xLocal[ numNodesPerFace ][ 3 ];
- for( localIndex a = 0; a < numNodesPerFace; ++a )
- {
- for( localIndex i = 0; i < 3; ++i )
+ real64 xLocal[ numNodesPerFace ][ 3 ];
+ for( localIndex a = 0; a < numNodesPerFace; ++a )
{
- xLocal[a][i] = X( facesToNodes( f, a ), i );
+ for( localIndex d = 0; d < 3; ++d )
+ {
+ xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d );
+ }
}
- }
- for( localIndex q = 0; q < numNodesPerFace; ++q )
- {
- real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
- RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes[f][q]], localIncrement );
+ for( localIndex q = 0; q < numNodesPerFace; ++q )
+ {
+ real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes( f, q )], localIncrement );
+ }
}
}
} ); // end loop over element
@@ -396,23 +395,23 @@ struct VelocityComputation
m_finiteElement.template computeFirstOrderStiffnessTermX( q, xLocal, [&] ( int i, int j, real32 dfx1, real32 dfx2, real32 dfx3 )
{
- flowx[j] += dfx1*p_np1[elemsToNodes[k][i]];
- flowy[j] += dfx2*p_np1[elemsToNodes[k][i]];
- flowz[j] += dfx3*p_np1[elemsToNodes[k][i]];
+ flowx[j] += dfx1*p_np1[elemsToNodes( k, i )];
+ flowy[j] += dfx2*p_np1[elemsToNodes( k, i )];
+ flowz[j] += dfx3*p_np1[elemsToNodes( k, i )];
} );
m_finiteElement.template computeFirstOrderStiffnessTermY( q, xLocal, [&] ( int i, int j, real32 dfy1, real32 dfy2, real32 dfy3 )
{
- flowx[j] += dfy1*p_np1[elemsToNodes[k][i]];
- flowy[j] += dfy2*p_np1[elemsToNodes[k][i]];
- flowz[j] += dfy3*p_np1[elemsToNodes[k][i]];
+ flowx[j] += dfy1*p_np1[elemsToNodes( k, i )];
+ flowy[j] += dfy2*p_np1[elemsToNodes( k, i )];
+ flowz[j] += dfy3*p_np1[elemsToNodes( k, i )];
} );
m_finiteElement.template computeFirstOrderStiffnessTermZ( q, xLocal, [&] ( int i, int j, real32 dfz1, real32 dfz2, real32 dfz3 )
{
- flowx[j] += dfz1*p_np1[elemsToNodes[k][i]];
- flowy[j] += dfz2*p_np1[elemsToNodes[k][i]];
- flowz[j] += dfz3*p_np1[elemsToNodes[k][i]];
+ flowx[j] += dfz1*p_np1[elemsToNodes( k, i )];
+ flowy[j] += dfz2*p_np1[elemsToNodes( k, i )];
+ flowz[j] += dfz3*p_np1[elemsToNodes( k, i )];
} );
@@ -551,8 +550,8 @@ struct PressureComputation
real32 diag=(auxx[i]+auyy[i]+auzz[i]);
uelemx[i]+=dt*diag;
- real32 const localIncrement = uelemx[i]/mass[elemsToNodes[k][i]];
- RAJA::atomicAdd< ATOMIC_POLICY >( &p_np1[elemsToNodes[k][i]], localIncrement );
+ real32 const localIncrement = uelemx[i]/mass[elemsToNodes( k, i )];
+ RAJA::atomicAdd< ATOMIC_POLICY >( &p_np1[elemsToNodes( k, i )], localIncrement );
}
//Source Injection
@@ -564,8 +563,8 @@ struct PressureComputation
{
for( localIndex i = 0; i < numNodesPerElem; ++i )
{
- real32 const localIncrement2 = dt*(sourceConstants[isrc][i]*sourceValue[cycleNumber][isrc])/(mass[elemsToNodes[k][i]]);
- RAJA::atomicAdd< ATOMIC_POLICY >( &p_np1[elemsToNodes[k][i]], localIncrement2 );
+ real32 const localIncrement2 = dt*(sourceConstants[isrc][i]*sourceValue[cycleNumber][isrc])/(mass[elemsToNodes( k, i )]);
+ RAJA::atomicAdd< ATOMIC_POLICY >( &p_np1[elemsToNodes( k, i )], localIncrement2 );
}
}
}
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp
index 98ae308427d..4f3ad22d174 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.cpp
@@ -60,6 +60,7 @@ void AcousticWaveEquationSEM::initializePreSubGroups()
void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
{
WaveSolverBase::registerDataOnMesh( meshBodies );
+
forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
MeshLevel & mesh,
arrayView1d< string const > const & )
@@ -71,10 +72,10 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
fields::Pressure_np1,
fields::PressureDoubleDerivative,
fields::ForcingRHS,
- fields::MassVector,
+ fields::MassVectorA,
fields::DampingVector,
fields::StiffnessVector,
- fields::FreeSurfaceNodeIndicator >( this->getName() );
+ fields::FreeSurfaceNodeIndicatorA >( getName() );
/// register PML auxiliary variables only when a PML is specified in the xml
if( m_usePML )
@@ -82,36 +83,31 @@ void AcousticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
nodeManager.registerField< fields::AuxiliaryVar1PML,
fields::AuxiliaryVar2PML,
fields::AuxiliaryVar3PML,
- fields::AuxiliaryVar4PML >( this->getName() );
+ fields::AuxiliaryVar4PML >( getName() );
nodeManager.getField< fields::AuxiliaryVar1PML >().resizeDimension< 1 >( 3 );
nodeManager.getField< fields::AuxiliaryVar2PML >().resizeDimension< 1 >( 3 );
}
FaceManager & faceManager = mesh.getFaceManager();
- faceManager.registerField< fields::FreeSurfaceFaceIndicator >( this->getName() );
+ faceManager.registerField< fields::FreeSurfaceFaceIndicatorA >( getName() );
ElementRegionManager & elemManager = mesh.getElemManager();
elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion )
{
- subRegion.registerField< fields::MediumVelocity >( this->getName() );
- subRegion.registerField< fields::PartialGradient >( this->getName() );
+ subRegion.registerField< fields::MediumVelocity >( getName() );
+ subRegion.registerField< fields::MediumDensityA >( getName() );
+ subRegion.registerField< fields::PartialGradient >( getName() );
} );
} );
}
-
void AcousticWaveEquationSEM::postProcessInput()
{
-
WaveSolverBase::postProcessInput();
-
- localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 );
-
- m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
-
+ m_pressureNp1AtReceivers.resize( m_nsamplesSeismoTrace, m_receiverCoordinates.size( 0 ) + 1 );
}
void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
@@ -144,15 +140,13 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
receiverConstants.setValues< EXEC_POLICY >( -1 );
receiverIsLocal.zero();
- real32 const timeSourceFrequency = this->m_timeSourceFrequency;
- localIndex const rickerOrder = this->m_rickerOrder;
arrayView2d< real32 > const sourceValue = m_sourceValue.toView();
real64 dt = 0;
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
{
EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
- if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+ if( subEvent->getEventName() == "/Solvers/" + getName() )
{
dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
}
@@ -176,16 +170,13 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
{
using FE_TYPE = TYPEOFREF( finiteElement );
- constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
localIndex const numFacesPerElem = elementSubRegion.numFacesPerElement();
-
{
GEOS_MARK_SCOPE( acousticWaveEquationSEMKernels::PrecomputeSourceAndReceiverKernel );
acousticWaveEquationSEMKernels::
PrecomputeSourceAndReceiverKernel::
launch< EXEC_POLICY, FE_TYPE >
( elementSubRegion.size(),
- numNodesPerElem,
numFacesPerElem,
X32,
elemGhostRank,
@@ -204,8 +195,9 @@ void AcousticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
receiverConstants,
sourceValue,
dt,
- timeSourceFrequency,
- rickerOrder );
+ m_timeSourceFrequency,
+ m_timeSourceDelay,
+ m_rickerOrder );
}
} );
} );
@@ -240,14 +232,11 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
WaveSolverBase::initializePostInitialConditionsPreSubGroups();
}
if( m_usePML )
- {
AcousticWaveEquationSEM::initializePML();
- }
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
- real64 const time = 0.0;
- applyFreeSurfaceBC( time, domain );
+ applyFreeSurfaceBC( 0.0, domain );
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
@@ -257,16 +246,14 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
+ ElementRegionManager & elemManager = mesh.getElemManager();
/// get the array of indicators: 1 if the face is on the boundary; 0 otherwise
arrayView1d< integer > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator();
- arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X32 = nodeManager.getField< fields::referencePosition32 >().toViewConst();
-
- /// get face to node map
- ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
+ arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst();
// mass matrix to be computed in this function
- arrayView1d< real32 > const mass = nodeManager.getField< fields::MassVector >();
+ arrayView1d< real32 > const mass = nodeManager.getField< fields::MassVectorA >();
{
GEOS_MARK_SCOPE( mass_zero );
mass.zero();
@@ -277,53 +264,57 @@ void AcousticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
GEOS_MARK_SCOPE( damping_zero );
damping.zero();
}
- /// get array of indicators: 1 if face is on the free surface; 0 otherwise
- arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >();
- mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
- CellElementSubRegion & elementSubRegion )
+ arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicatorA >();
+ ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
+
+ elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
+ CellElementSubRegion & elementSubRegion )
{
+ finiteElement::FiniteElementBase const &
+ fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
+
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
- arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
+ arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();
+
+ computeTargetNodeSet( elemsToNodes, elementSubRegion.size(), fe.getNumQuadraturePoints() );
+
arrayView1d< real32 const > const velocity = elementSubRegion.getField< fields::MediumVelocity >();
+ arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensityA >();
/// Partial gradient if gradient as to be computed
arrayView1d< real32 > grad = elementSubRegion.getField< fields::PartialGradient >();
- {
- grad.zero();
- }
- finiteElement::FiniteElementBase const &
- fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
+ grad.zero();
+
finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );
- {
- GEOS_MARK_SCOPE( MassMatrixKernel );
- acousticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement );
- kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
- X32,
- elemsToNodes,
- velocity,
- mass );
- }
+ acousticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement );
+ kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
+ nodeCoords,
+ elemsToNodes,
+ velocity,
+ density,
+ mass );
{
GEOS_MARK_SCOPE( DampingMatrixKernel );
acousticWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );
-
- kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
- X32,
- facesToElements,
+ kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
+ nodeCoords,
+ elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
velocity,
+ density,
damping );
}
} );
} );
} );
+ WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
}
@@ -343,10 +334,10 @@ void AcousticWaveEquationSEM::applyFreeSurfaceBC( real64 time, DomainPartition &
ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst();
/// array of indicators: 1 if a face is on on free surface; 0 otherwise
- arrayView1d< localIndex > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >();
+ arrayView1d< localIndex > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicatorA >();
/// array of indicators: 1 if a node is on on free surface; 0 otherwise
- arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicator >();
+ arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicatorA >();
//freeSurfaceFaceIndicator.zero();
//freeSurfaceNodeIndicator.zero();
@@ -423,7 +414,7 @@ void AcousticWaveEquationSEM::initializePML()
} );
/// Now compute the PML parameters above internally
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
arrayView1d< string const > const & )
@@ -827,20 +818,14 @@ real64 AcousticWaveEquationSEM::explicitStepForward( real64 const & time_n,
"An error occured while writting "<< fileName,
InputError );
}
-
}
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
- {
- p_nm1[a] = p_n[a];
- p_n[a] = p_np1[a];
- } );
+ prepareNextTimestep( mesh );
} );
return dtOut;
}
-
real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
real64 const & dt,
integer cycleNumber,
@@ -848,6 +833,7 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
bool computeGradient )
{
real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain );
+
forDiscretizationOnMeshTargets( domain.getMeshBodies(),
[&] ( string const &,
MeshLevel & mesh,
@@ -855,13 +841,13 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
{
NodeManager & nodeManager = mesh.getNodeManager();
- arrayView1d< real32 const > const mass = nodeManager.getField< fields::MassVector >();
+ arrayView1d< real32 const > const mass = nodeManager.getField< fields::MassVectorA >();
arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::Pressure_nm1 >();
arrayView1d< real32 > const p_n = nodeManager.getField< fields::Pressure_n >();
arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() );
int const maxCycle = int(round( maxTime/dt ));
@@ -920,180 +906,218 @@ real64 AcousticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
} );
}
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
- {
- p_nm1[a] = p_n[a];
- p_n[a] = p_np1[a];
- } );
+ prepareNextTimestep( mesh );
} );
return dtOut;
}
-real64 AcousticWaveEquationSEM::explicitStepInternal( real64 const & time_n,
- real64 const & dt,
- integer cycleNumber,
- DomainPartition & domain )
+void AcousticWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh )
{
- GEOS_MARK_FUNCTION;
+ NodeManager & nodeManager = mesh.getNodeManager();
- GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc );
+ arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::Pressure_nm1 >();
+ arrayView1d< real32 > const p_n = nodeManager.getField< fields::Pressure_n >();
+ arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
- forDiscretizationOnMeshTargets( domain.getMeshBodies(),
- [&] ( string const &,
- MeshLevel & mesh,
- arrayView1d< string const > const & regionNames )
+ arrayView1d< real32 > const stiffnessVector = nodeManager.getField< fields::StiffnessVector >();
+ arrayView1d< real32 > const rhs = nodeManager.getField< fields::ForcingRHS >();
+
+ SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst();
+
+ forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
{
- NodeManager & nodeManager = mesh.getNodeManager();
+ localIndex const a = solverTargetNodesSet[n];
- arrayView1d< real32 const > const mass = nodeManager.getField< fields::MassVector >();
- arrayView1d< real32 const > const damping = nodeManager.getField< fields::DampingVector >();
+ p_nm1[a] = p_n[a];
+ p_n[a] = p_np1[a];
- arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::Pressure_nm1 >();
- arrayView1d< real32 > const p_n = nodeManager.getField< fields::Pressure_n >();
- arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
+ stiffnessVector[a] = rhs[a] = 0.0;
+ } );
+}
- arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicator >();
- arrayView1d< real32 > const stiffnessVector = nodeManager.getField< fields::StiffnessVector >();
- arrayView1d< real32 > const rhs = nodeManager.getField< fields::ForcingRHS >();
+void AcousticWaveEquationSEM::computeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer cycleNumber,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+{
+ NodeManager & nodeManager = mesh.getNodeManager();
- bool const usePML = m_usePML;
+ arrayView1d< real32 const > const mass = nodeManager.getField< fields::MassVectorA >();
+ arrayView1d< real32 const > const damping = nodeManager.getField< fields::DampingVector >();
+
+ arrayView1d< real32 > const p_nm1 = nodeManager.getField< fields::Pressure_nm1 >();
+ arrayView1d< real32 > const p_n = nodeManager.getField< fields::Pressure_n >();
+ arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
- auto kernelFactory = acousticWaveEquationSEMKernels::ExplicitAcousticSEMFactory( dt );
+ arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicatorA >();
+ arrayView1d< real32 > const stiffnessVector = nodeManager.getField< fields::StiffnessVector >();
+ arrayView1d< real32 > const rhs = nodeManager.getField< fields::ForcingRHS >();
- finiteElement::
- regionBasedKernelApplication< EXEC_POLICY,
- constitutive::NullModel,
- CellElementSubRegion >( mesh,
- regionNames,
- getDiscretizationName(),
- "",
- kernelFactory );
+ auto kernelFactory = acousticWaveEquationSEMKernels::ExplicitAcousticSEMFactory( dt );
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
- real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() );
- integer const cycleForSource = int(round( -minTime/dt + cycleNumber ));
- //std::cout<<"cycle GEOSX = "<( mesh,
+ regionNames,
+ getDiscretizationName(),
+ "",
+ kernelFactory );
- /// calculate your time integrators
- real64 const dt2 = dt*dt;
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
+ real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() );
+ integer const cycleForSource = int(round( -minTime/dt + cycleNumber ));
+ addSourceToRightHandSide( cycleForSource, rhs );
- if( !usePML )
+ /// calculate your time integrators
+ real64 const dt2 = pow( dt, 2 );
+
+ SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst();
+ if( !m_usePML )
+ {
+ GEOS_MARK_SCOPE ( updateP );
+ forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
{
- GEOS_MARK_SCOPE ( updateP );
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
+ localIndex const a = solverTargetNodesSet[n];
+ if( freeSurfaceNodeIndicator[a] != 1 )
{
- if( freeSurfaceNodeIndicator[a] != 1 )
+ p_np1[a] = p_n[a];
+ p_np1[a] *= 2.0*mass[a];
+ p_np1[a] -= (mass[a]-0.5*dt*damping[a])*p_nm1[a];
+ p_np1[a] += dt2*(rhs[a]-stiffnessVector[a]);
+ p_np1[a] /= mass[a]+0.5*dt*damping[a];
+ }
+ } );
+ }
+ else
+ {
+ parametersPML const & param = getReference< parametersPML >( viewKeyStruct::parametersPMLString() );
+ arrayView2d< real32 > const v_n = nodeManager.getField< fields::AuxiliaryVar1PML >();
+ arrayView2d< real32 > const grad_n = nodeManager.getField< fields::AuxiliaryVar2PML >();
+ arrayView1d< real32 > const divV_n = nodeManager.getField< fields::AuxiliaryVar3PML >();
+ arrayView1d< real32 > const u_n = nodeManager.getField< fields::AuxiliaryVar4PML >();
+ arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const
+ X32 = nodeManager.getField< fields::referencePosition32 >().toViewConst();
+
+ real32 const xMin[ 3 ] = {param.xMinPML[0], param.xMinPML[1], param.xMinPML[2]};
+ real32 const xMax[ 3 ] = {param.xMaxPML[0], param.xMaxPML[1], param.xMaxPML[2]};
+ real32 const dMin[ 3 ] = {param.thicknessMinXYZPML[0], param.thicknessMinXYZPML[1], param.thicknessMinXYZPML[2]};
+ real32 const dMax[ 3 ] = {param.thicknessMaxXYZPML[0], param.thicknessMaxXYZPML[1], param.thicknessMaxXYZPML[2]};
+ real32 const cMin[ 3 ] = {param.waveSpeedMinXYZPML[0], param.waveSpeedMinXYZPML[1], param.waveSpeedMinXYZPML[2]};
+ real32 const cMax[ 3 ] = {param.waveSpeedMaxXYZPML[0], param.waveSpeedMaxXYZPML[1], param.waveSpeedMaxXYZPML[2]};
+ real32 const r = param.reflectivityPML;
+
+ /// apply the main function to update some of the PML auxiliary variables
+ /// Compute (divV) and (B.pressureGrad - C.auxUGrad) vectors for the PML region
+ applyPML( time_n, domain );
+
+ GEOS_MARK_SCOPE ( updatePWithPML );
+ forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
+ {
+ localIndex const a = solverTargetNodesSet[n];
+ if( freeSurfaceNodeIndicator[a] != 1 )
+ {
+ real32 sigma[3];
+ real32 xLocal[ 3 ];
+
+ for( integer i=0; i<3; ++i )
{
- p_np1[a] = p_n[a];
- p_np1[a] *= 2.0*mass[a];
- p_np1[a] -= (mass[a]-0.5*dt*damping[a])*p_nm1[a];
- p_np1[a] += dt2*(rhs[a]-stiffnessVector[a]);
- p_np1[a] /= mass[a]+0.5*dt*damping[a];
+ xLocal[i] = X32[a][i];
}
- } );
- }
- else
- {
- parametersPML const & param = getReference< parametersPML >( viewKeyStruct::parametersPMLString() );
- arrayView2d< real32 > const v_n = nodeManager.getField< fields::AuxiliaryVar1PML >();
- arrayView2d< real32 > const grad_n = nodeManager.getField< fields::AuxiliaryVar2PML >();
- arrayView1d< real32 > const divV_n = nodeManager.getField< fields::AuxiliaryVar3PML >();
- arrayView1d< real32 > const u_n = nodeManager.getField< fields::AuxiliaryVar4PML >();
- arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X32 = nodeManager.getField< fields::referencePosition32 >().toViewConst();
-
- real32 const xMin[ 3 ] = {param.xMinPML[0], param.xMinPML[1], param.xMinPML[2]};
- real32 const xMax[ 3 ] = {param.xMaxPML[0], param.xMaxPML[1], param.xMaxPML[2]};
- real32 const dMin[ 3 ] = {param.thicknessMinXYZPML[0], param.thicknessMinXYZPML[1], param.thicknessMinXYZPML[2]};
- real32 const dMax[ 3 ] = {param.thicknessMaxXYZPML[0], param.thicknessMaxXYZPML[1], param.thicknessMaxXYZPML[2]};
- real32 const cMin[ 3 ] = {param.waveSpeedMinXYZPML[0], param.waveSpeedMinXYZPML[1], param.waveSpeedMinXYZPML[2]};
- real32 const cMax[ 3 ] = {param.waveSpeedMaxXYZPML[0], param.waveSpeedMaxXYZPML[1], param.waveSpeedMaxXYZPML[2]};
- real32 const r = param.reflectivityPML;
- /// apply the main function to update some of the PML auxiliary variables
- /// Compute (divV) and (B.pressureGrad - C.auxUGrad) vectors for the PML region
- applyPML( time_n, domain );
+ acousticWaveEquationSEMKernels::PMLKernelHelper::computeDampingProfilePML(
+ xLocal,
+ xMin,
+ xMax,
+ dMin,
+ dMax,
+ cMin,
+ cMax,
+ r,
+ sigma );
- GEOS_MARK_SCOPE ( updatePWithPML );
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
- {
- if( freeSurfaceNodeIndicator[a] != 1 )
- {
- real32 sigma[3];
- real32 xLocal[ 3 ];
+ real32 const alpha = sigma[0] + sigma[1] + sigma[2];
- for( integer i=0; i<3; ++i )
- {
- xLocal[i] = X32[a][i];
- }
+ p_np1[a] = dt2*( (rhs[a] - stiffnessVector[a])/mass[a] - divV_n[a])
+ - (1 - 0.5*alpha*dt)*p_nm1[a]
+ + 2*p_n[a];
- acousticWaveEquationSEMKernels::PMLKernelHelper::computeDampingProfilePML(
- xLocal,
- xMin,
- xMax,
- dMin,
- dMax,
- cMin,
- cMax,
- r,
- sigma );
+ p_np1[a] = p_np1[a] / (1 + 0.5*alpha*dt);
- real32 const alpha = sigma[0] + sigma[1] + sigma[2];
+ for( integer i=0; i<3; ++i )
+ {
+ v_n[a][i] = (1 - dt*sigma[i])*v_n[a][i] - dt*grad_n[a][i];
+ }
+ u_n[a] += dt*p_n[a];
+ }
+ } );
+ }
+}
- p_np1[a] = dt2*( (rhs[a] - stiffnessVector[a])/mass[a] - divV_n[a])
- - (1 - 0.5*alpha*dt)*p_nm1[a]
- + 2*p_n[a];
+void AcousticWaveEquationSEM::synchronizeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer const,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & )
+{
+ NodeManager & nodeManager = mesh.getNodeManager();
- p_np1[a] = p_np1[a] / (1 + 0.5*alpha*dt);
+ arrayView1d< real32 > const p_n = nodeManager.getField< fields::Pressure_n >();
+ arrayView1d< real32 > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
- for( integer i=0; i<3; ++i )
- {
- v_n[a][i] = (1 - dt*sigma[i])*v_n[a][i] - dt*grad_n[a][i];
- }
- u_n[a] += dt*p_n[a];
- }
- } );
- }
+ arrayView1d< real32 > const stiffnessVector = nodeManager.getField< fields::StiffnessVector >();
+ arrayView1d< real32 > const rhs = nodeManager.getField< fields::ForcingRHS >();
- /// synchronize pressure fields
- FieldIdentifiers fieldsToBeSync;
- fieldsToBeSync.addFields( FieldLocation::Node, { fields::Pressure_np1::key() } );
+ /// synchronize pressure fields
+ FieldIdentifiers fieldsToBeSync;
+ fieldsToBeSync.addFields( FieldLocation::Node, { fields::Pressure_np1::key() } );
- if( usePML )
- {
- fieldsToBeSync.addFields( FieldLocation::Node, {
- fields::AuxiliaryVar1PML::key(),
- fields::AuxiliaryVar4PML::key() } );
- }
+ if( m_usePML )
+ {
+ fieldsToBeSync.addFields( FieldLocation::Node, {
+ fields::AuxiliaryVar1PML::key(),
+ fields::AuxiliaryVar4PML::key() } );
+ }
- CommunicationTools & syncFields = CommunicationTools::getInstance();
- syncFields.synchronizeFields( fieldsToBeSync,
- mesh,
- domain.getNeighbors(),
- true );
+ CommunicationTools & syncFields = CommunicationTools::getInstance();
+ syncFields.synchronizeFields( fieldsToBeSync,
+ mesh,
+ domain.getNeighbors(),
+ true );
+ /// compute the seismic traces since last step.
+ arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
- /// compute the seismic traces since last step.
- arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
- if( time_n >= 0 )
- {
- computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers );
- }
- /// prepare next step
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
- {
- stiffnessVector[a] = 0.0;
- rhs[a] = 0.0;
- } );
+ computeAllSeismoTraces( time_n, dt, p_np1, p_n, pReceivers );
+ incrementIndexSeismoTrace( time_n );
- if( usePML )
- {
- arrayView2d< real32 > const grad_n = nodeManager.getField< fields::AuxiliaryVar2PML >();
- arrayView1d< real32 > const divV_n = nodeManager.getField< fields::AuxiliaryVar3PML >();
- grad_n.zero();
- divV_n.zero();
- }
+ if( m_usePML )
+ {
+ arrayView2d< real32 > const grad_n = nodeManager.getField< fields::AuxiliaryVar2PML >();
+ arrayView1d< real32 > const divV_n = nodeManager.getField< fields::AuxiliaryVar3PML >();
+ grad_n.zero();
+ divV_n.zero();
+ }
+}
+real64 AcousticWaveEquationSEM::explicitStepInternal( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain )
+{
+ GEOS_MARK_FUNCTION;
+
+ GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc );
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+ {
+ computeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames );
+ synchronizeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames );
} );
return dt;
@@ -1116,40 +1140,12 @@ void AcousticWaveEquationSEM::cleanup( real64 const time_n,
NodeManager & nodeManager = mesh.getNodeManager();
arrayView1d< real32 const > const p_n = nodeManager.getField< fields::Pressure_n >();
arrayView1d< real32 const > const p_np1 = nodeManager.getField< fields::Pressure_np1 >();
- arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
- computeAllSeismoTraces( time_n, 0, p_np1, p_n, pReceivers );
- } );
-}
-
-void AcousticWaveEquationSEM::computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers )
-{
+ arrayView2d< real32 > const pReceivers = m_pressureNp1AtReceivers.toView();
+ computeAllSeismoTraces( time_n, 0.0, p_np1, p_n, pReceivers );
- /*
- * In forward case we compute seismo if time_n + dt is the first time
- * step after the timeSeismo to write.
- *
- * time_n timeSeismo time_n + dt
- * ---|--------------|-------------|
- *
- * In backward (time_n goes decreasing) case we compute seismo if
- * time_n is the last time step before the timeSeismo to write.
- *
- * time_n - dt timeSeismo time_n
- * ---|--------------|-------------|
- */
- for( real64 timeSeismo;
- (m_forward)?((timeSeismo = m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + dt + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace):
- ((timeSeismo = m_dtSeismoTrace*(m_nsamplesSeismoTrace-m_indexSeismoTrace-1)) >= (time_n - dt - epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace);
- m_indexSeismoTrace++ )
- {
- WaveSolverUtils::computeSeismoTrace( time_n, (m_forward)?dt:-dt, timeSeismo, (m_forward)?m_indexSeismoTrace:(m_nsamplesSeismoTrace-m_indexSeismoTrace-1), m_receiverNodeIds, m_receiverConstants,
- m_receiverIsLocal,
- m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers );
- }
+ WaveSolverUtils::writeSeismoTrace( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, pReceivers );
+ } );
}
REGISTER_CATALOG_ENTRY( SolverBase, AcousticWaveEquationSEM, string const &, dataRepository::Group * const )
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp
index 657bac9a1b5..94e06c01660 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEM.hpp
@@ -31,6 +31,9 @@ class AcousticWaveEquationSEM : public WaveSolverBase
{
public:
+ /// String used to form the solverName used to register solvers in CoupledSolver
+ static string coupledSolverAttributePrefix() { return "acoustic"; }
+
using EXEC_POLICY = parallelDevicePolicy< >;
using ATOMIC_POLICY = AtomicPolicy< EXEC_POLICY >;
@@ -81,21 +84,6 @@ class AcousticWaveEquationSEM : public WaveSolverBase
*/
virtual void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs );
- /**
- * TODO: move implementation into WaveSolverBase
- * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
- * @param time_n the time corresponding to the field values pressure_n
- * @param dt the simulation timestep
- * @param var_np1 the field values at time_n + dt
- * @param var_n the field values at time_n
- * @param varAtReceivers the array holding the trace values, where the output is written
- */
- virtual void computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers );
-
/**
* @brief Initialize Perfectly Matched Layer (PML) information
@@ -110,14 +98,6 @@ class AcousticWaveEquationSEM : public WaveSolverBase
struct viewKeyStruct : WaveSolverBase::viewKeyStruct
{
- static constexpr char const * sourceNodeIdsString() { return "sourceNodeIds"; }
- static constexpr char const * sourceConstantsString() { return "sourceConstants"; }
- static constexpr char const * sourceIsAccessibleString() { return "sourceIsAccessible"; }
-
- static constexpr char const * receiverNodeIdsString() { return "receiverNodeIds"; }
- static constexpr char const * receiverConstantsString() {return "receiverConstants"; }
- static constexpr char const * receiverIsLocalString() { return "receiverIsLocal"; }
-
static constexpr char const * pressureNp1AtReceiversString() { return "pressureNp1AtReceivers"; }
} waveEquationViewKeys;
@@ -136,6 +116,22 @@ class AcousticWaveEquationSEM : public WaveSolverBase
integer const cycleNumber,
DomainPartition & domain );
+ void computeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames );
+
+ void synchronizeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames );
+
+ void prepareNextTimestep( MeshLevel & mesh );
+
protected:
virtual void postProcessInput() override final;
@@ -206,8 +202,8 @@ DECLARE_FIELD( ForcingRHS,
WRITE_AND_READ,
"RHS" );
-DECLARE_FIELD( MassVector,
- "massVector",
+DECLARE_FIELD( MassVectorA,
+ "massVectorA",
array1d< real32 >,
0,
NOPLOT,
@@ -230,6 +226,14 @@ DECLARE_FIELD( MediumVelocity,
WRITE_AND_READ,
"Medium velocity of the cell" );
+DECLARE_FIELD( MediumDensityA,
+ "mediumDensityA",
+ array1d< real32 >,
+ 0,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Medium density of the cell" );
+
DECLARE_FIELD( StiffnessVector,
"stiffnessVector",
array1d< real32 >,
@@ -238,16 +242,16 @@ DECLARE_FIELD( StiffnessVector,
WRITE_AND_READ,
"Stiffness vector contains R_h*Pressure_n." );
-DECLARE_FIELD( FreeSurfaceFaceIndicator,
- "freeSurfaceFaceIndicator",
+DECLARE_FIELD( FreeSurfaceFaceIndicatorA,
+ "freeSurfaceFaceIndicatorA",
array1d< localIndex >,
0,
NOPLOT,
WRITE_AND_READ,
"Free surface indicator, 1 if a face is on free surface 0 otherwise." );
-DECLARE_FIELD( FreeSurfaceNodeIndicator,
- "freeSurfaceNodeIndicator",
+DECLARE_FIELD( FreeSurfaceNodeIndicatorA,
+ "freeSurfaceNodeIndicatorA",
array1d< localIndex >,
0,
NOPLOT,
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp
index 7293a2a79b6..4de6fb92cb3 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/AcousticWaveEquationSEMKernel.hpp
@@ -40,7 +40,6 @@ struct PrecomputeSourceAndReceiverKernel
* @tparam EXEC_POLICY execution policy
* @tparam FE_TYPE finite element type
* @param[in] size the number of cells in the subRegion
- * @param[in] numNodesPerElem number of nodes per element
* @param[in] X coordinates of the nodes
* @param[in] elemsToNodes map from element to nodes
* @param[in] elemsToFaces map from element to faces
@@ -58,7 +57,6 @@ struct PrecomputeSourceAndReceiverKernel
template< typename EXEC_POLICY, typename FE_TYPE >
static void
launch( localIndex const size,
- localIndex const numNodesPerElem,
localIndex const numFacesPerElem,
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
arrayView1d< integer const > const elemGhostRank,
@@ -78,11 +76,12 @@ struct PrecomputeSourceAndReceiverKernel
arrayView2d< real32 > const sourceValue,
real64 const dt,
real32 const timeSourceFrequency,
+ real32 const timeSourceDelay,
localIndex const rickerOrder )
{
-
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
{
+ constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
real64 const center[3] = { elemCenter[k][0],
elemCenter[k][1],
elemCenter[k][2] };
@@ -121,14 +120,13 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- sourceNodeIds[isrc][a] = elemsToNodes[k][a];
+ sourceNodeIds[isrc][a] = elemsToNodes( k, a );
sourceConstants[isrc][a] = Ntest[a];
}
for( localIndex cycle = 0; cycle < sourceValue.size( 0 ); ++cycle )
{
- real64 const time = cycle*dt;
- sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( time, timeSourceFrequency, rickerOrder );
+ sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( cycle * dt, timeSourceFrequency, timeSourceDelay, rickerOrder );
}
}
}
@@ -164,12 +162,12 @@ struct PrecomputeSourceAndReceiverKernel
receiverIsLocal[ircv] = 1;
- real64 Ntest[FE_TYPE::numNodes];
+ real64 Ntest[numNodesPerElem];
FE_TYPE::calcN( coordsOnRefElem, Ntest );
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- receiverNodeIds[ircv][a] = elemsToNodes[k][a];
+ receiverNodeIds[ircv][a] = elemsToNodes( k, a );
receiverConstants[ircv][a] = Ntest[a];
}
}
@@ -198,6 +196,7 @@ struct MassMatrixKernel
* @param[in] X coordinates of the nodes
* @param[in] elemsToNodes map from element to nodes
* @param[in] velocity cell-wise velocity
+ * @param[in] density cell-wise density
* @param[out] mass diagonal of the mass matrix
*/
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
@@ -206,28 +205,29 @@ struct MassMatrixKernel
arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes,
arrayView1d< real32 const > const velocity,
+ arrayView1d< real32 const > const density,
arrayView1d< real32 > const mass )
{
- forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
constexpr localIndex numQuadraturePointsPerElem = FE_TYPE::numQuadraturePoints;
- real32 const invC2 = 1.0 / ( velocity[k] * velocity[k] );
+ real32 const invC2 = 1.0 / ( density[e] * velocity[e] * velocity[e] );
real64 xLocal[ numNodesPerElem ][ 3 ];
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
for( localIndex i = 0; i < 3; ++i )
{
- xLocal[a][i] = X( elemsToNodes( k, a ), i );
+ xLocal[a][i] = X( elemsToNodes( e, a ), i );
}
}
for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
{
real32 const localIncrement = invC2 * m_finiteElement.computeMassTerm( q, xLocal );
- RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement );
}
} ); // end loop over element
}
@@ -237,6 +237,13 @@ struct MassMatrixKernel
};
+using KERNEL_SEQ_POLICY =
+ RAJA::KernelPolicy<
+ RAJA::statement::For< 0, RAJA::seq_exec,
+ RAJA::statement::Lambda< 0 >
+ >
+ >;
+
template< typename FE_TYPE >
struct DampingMatrixKernel
{
@@ -250,56 +257,54 @@ struct DampingMatrixKernel
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
- * @param[in] X coordinates of the nodes
+ * @param[in] nodeCoords coordinates of the nodes
* @param[in] facesToElems map from faces to elements
* @param[in] facesToNodes map from face to nodes
* @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise
* @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise
* @param[in] velocity cell-wise velocity
+ * @param[in] density cell-wise density
* @param[out] damping diagonal of the damping matrix
*/
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
void
launch( localIndex const size,
- arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
- arrayView2d< localIndex const > const facesToElems,
+ arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
+ arrayView2d< localIndex const > const elemsToFaces,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView1d< integer const > const facesDomainBoundaryIndicator,
arrayView1d< localIndex const > const freeSurfaceFaceIndicator,
arrayView1d< real32 const > const velocity,
+ arrayView1d< real32 const > const density,
arrayView1d< real32 > const damping )
{
- forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
- // face on the domain boundary and not on free surface
- if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
+ for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i )
{
- localIndex k = facesToElems( f, 0 );
- if( k < 0 )
- {
- k = facesToElems( f, 1 );
- }
-
- real32 const alpha = 1.0 / velocity[k];
-
- constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
-
- real64 xLocal[ numNodesPerFace ][ 3 ];
- for( localIndex a = 0; a < numNodesPerFace; ++a )
+ localIndex const f = elemsToFaces( e, i );
+ // face on the domain boundary and not on free surface
+ if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
{
- for( localIndex i = 0; i < 3; ++i )
+ constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
+ real64 xLocal[ numNodesPerFace ][ 3 ];
+ for( localIndex a = 0; a < numNodesPerFace; ++a )
{
- xLocal[a][i] = X( facesToNodes( f, a ), i );
+ for( localIndex d = 0; d < 3; ++d )
+ {
+ xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d );
+ }
}
- }
+ real32 const alpha = 1.0 / (density[e] * velocity[e]);
- for( localIndex q = 0; q < numNodesPerFace; ++q )
- {
- real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
- RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes[f][q]], localIncrement );
+ for( localIndex q = 0; q < numNodesPerFace; ++q )
+ {
+ real32 const localIncrement = alpha * m_finiteElement.computeDampingTerm( q, xLocal );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &damping[facesToNodes( f, q )], localIncrement );
+ }
}
}
- } ); // end loop over element
+ } );
}
/// The finite element space/discretization object for the element type in the subRegion
@@ -323,7 +328,7 @@ struct PMLKernelHelper
*/
GEOS_HOST_DEVICE
inline
- static void computeDampingProfilePML( real32 const (&xLocal)[3],
+ static void computeDampingProfilePML( real32 const (&xLocal )[3],
real32 const (&xMin)[3],
real32 const (&xMax)[3],
real32 const (&dMin)[3],
@@ -738,6 +743,7 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE,
m_X( nodeManager.getField< fields::referencePosition32 >() ),
m_p_n( nodeManager.getField< fields::Pressure_n >() ),
m_stiffnessVector( nodeManager.getField< fields::StiffnessVector >() ),
+ m_density( elementSubRegion.template getField< fields::MediumDensityA >() ),
m_dt( dt )
{
GEOS_UNUSED_VAR( edgeManager );
@@ -802,8 +808,9 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE,
{
m_finiteElementSpace.template computeStiffnessTerm( q, stack.xLocal, [&] ( int i, int j, real64 val )
{
- real32 const localIncrement = val*m_p_n[m_elemsToNodes[k][j]];
- RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector[m_elemsToNodes[k][i]], localIncrement );
+ real32 invDensity = 1./m_density[k];
+ real32 const localIncrement = invDensity*val*m_p_n[m_elemsToNodes( k, j )];
+ RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVector[m_elemsToNodes( k, i )], localIncrement );
} );
}
@@ -817,6 +824,9 @@ class ExplicitAcousticSEM : public finiteElement::KernelBase< SUBREGION_TYPE,
/// The array containing the product of the stiffness matrix and the nodal pressure.
arrayView1d< real32 > const m_stiffnessVector;
+ /// The array containing the cell-wise density
+ arrayView1d< real32 const > const m_density;
+
/// The time increment for this time integration step.
real64 const m_dt;
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp
index 1eb94a47555..0b7035dd437 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.cpp
@@ -89,21 +89,11 @@ ElasticFirstOrderWaveEquationSEM::ElasticFirstOrderWaveEquationSEM( const std::s
setSizedFromParent( 0 ).
setDescription( "Element containing the sources" );
- registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ).
- setInputFlag( InputFlags::FALSE ).
- setSizedFromParent( 0 ).
- setDescription( "Element containing the receivers" );
-
registerWrapper( viewKeyStruct::sourceRegionString(), &m_sourceRegion ).
setInputFlag( InputFlags::FALSE ).
setSizedFromParent( 0 ).
setDescription( "Region containing the sources" );
- registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ).
- setInputFlag( InputFlags::FALSE ).
- setSizedFromParent( 0 ).
- setDescription( "Region containing the receivers" );
-
}
ElasticFirstOrderWaveEquationSEM::~ElasticFirstOrderWaveEquationSEM()
@@ -136,27 +126,27 @@ void ElasticFirstOrderWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
wavesolverfields::DampingVectorx,
wavesolverfields::DampingVectory,
wavesolverfields::DampingVectorz,
- wavesolverfields::FreeSurfaceNodeIndicator >( this->getName() );
+ wavesolverfields::FreeSurfaceNodeIndicator >( getName() );
FaceManager & faceManager = mesh.getFaceManager();
- faceManager.registerField< wavesolverfields::FreeSurfaceFaceIndicator >( this->getName() );
+ faceManager.registerField< wavesolverfields::FreeSurfaceFaceIndicator >( getName() );
ElementRegionManager & elemManager = mesh.getElemManager();
elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion )
{
- subRegion.registerField< wavesolverfields::MediumVelocityVp >( this->getName() );
- subRegion.registerField< wavesolverfields::MediumVelocityVs >( this->getName() );
- subRegion.registerField< wavesolverfields::MediumDensity >( this->getName() );
- subRegion.registerField< wavesolverfields::Lambda >( this->getName() );
- subRegion.registerField< wavesolverfields::Mu >( this->getName() );
-
- subRegion.registerField< wavesolverfields::Stresstensorxx >( this->getName());
- subRegion.registerField< wavesolverfields::Stresstensoryy >( this->getName());
- subRegion.registerField< wavesolverfields::Stresstensorzz >( this->getName());
- subRegion.registerField< wavesolverfields::Stresstensorxy >( this->getName());
- subRegion.registerField< wavesolverfields::Stresstensorxz >( this->getName());
- subRegion.registerField< wavesolverfields::Stresstensoryz >( this->getName());
+ subRegion.registerField< wavesolverfields::MediumVelocityVp >( getName() );
+ subRegion.registerField< wavesolverfields::MediumVelocityVs >( getName() );
+ subRegion.registerField< wavesolverfields::MediumDensity >( getName() );
+ subRegion.registerField< wavesolverfields::Lambda >( getName() );
+ subRegion.registerField< wavesolverfields::Mu >( getName() );
+
+ subRegion.registerField< wavesolverfields::Stresstensorxx >( getName());
+ subRegion.registerField< wavesolverfields::Stresstensoryy >( getName());
+ subRegion.registerField< wavesolverfields::Stresstensorzz >( getName());
+ subRegion.registerField< wavesolverfields::Stresstensorxy >( getName());
+ subRegion.registerField< wavesolverfields::Stresstensorxz >( getName());
+ subRegion.registerField< wavesolverfields::Stresstensoryz >( getName());
finiteElement::FiniteElementBase const & fe = subRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
@@ -181,8 +171,6 @@ void ElasticFirstOrderWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
} );
}
-
-
void ElasticFirstOrderWaveEquationSEM::postProcessInput()
{
@@ -196,18 +184,16 @@ void ElasticFirstOrderWaveEquationSEM::postProcessInput()
m_rcvElem.resize( numReceiversGlobal );
m_receiverRegion.resize( numReceiversGlobal );
- m_displacementxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_displacementyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_displacementzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
-
- m_sigmaxxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_sigmayyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_sigmazzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_sigmaxyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_sigmaxzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_sigmayzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
-
+ m_displacementxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_displacementyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_displacementzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_sigmaxxNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_sigmayyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_sigmazzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_sigmaxyNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_sigmaxzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_sigmayzNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
}
@@ -241,16 +227,14 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve
receiverConstants.setValues< serialPolicy >( -1 );
receiverIsLocal.zero();
- real32 const timeSourceFrequency = this->m_timeSourceFrequency;
- localIndex const rickerOrder = this->m_rickerOrder;
arrayView2d< real32 > const sourceValue = m_sourceValue.toView();
real64 dt = 0;
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
{
EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
- if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+ if( subEvent->getEventName() == "/Solvers/" + getName() )
{
dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
}
@@ -306,8 +290,9 @@ void ElasticFirstOrderWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLeve
receiverRegion,
sourceValue,
dt,
- timeSourceFrequency,
- rickerOrder );
+ m_timeSourceFrequency,
+ m_timeSourceDelay,
+ m_rickerOrder );
} );
} );
}
@@ -340,10 +325,9 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou
WaveSolverBase::initializePostInitialConditionsPreSubGroups();
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
- real64 const time = 0.0;
- applyFreeSurfaceBC( time, domain );
+ applyFreeSurfaceBC( 0.0, domain );
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
@@ -381,7 +365,7 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou
{
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
- arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
+ arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();
arrayView1d< real32 > const density = elementSubRegion.getField< wavesolverfields::MediumDensity >();
arrayView1d< real32 > const velocityVp = elementSubRegion.getField< wavesolverfields::MediumVelocityVp >();
arrayView1d< real32 > const velocityVs = elementSubRegion.getField< wavesolverfields::MediumVelocityVs >();
@@ -405,9 +389,9 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou
elasticFirstOrderWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );
- kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
+ kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
X,
- facesToElements,
+ elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
@@ -421,6 +405,8 @@ void ElasticFirstOrderWaveEquationSEM::initializePostInitialConditionsPreSubGrou
} );
} );
} );
+
+ WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
}
@@ -625,15 +611,14 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti
uy_np1,
uz_np1 );
-
} );
- arrayView2d< real32 > const sigmaxxReceivers = m_sigmaxxNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmayyReceivers = m_sigmayyNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmazzReceivers = m_sigmazzNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmaxyReceivers = m_sigmaxyNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmaxzReceivers = m_sigmaxzNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmayzReceivers = m_sigmayzNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmaxxReceivers = m_sigmaxxNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmayyReceivers = m_sigmayyNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmazzReceivers = m_sigmazzNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmaxyReceivers = m_sigmaxyNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmaxzReceivers = m_sigmaxzNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmayzReceivers = m_sigmayzNp1AtReceivers.toView();
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressxx, stressxx, sigmaxxReceivers );
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressyy, stressyy, sigmayyReceivers );
@@ -641,8 +626,6 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressxy, stressxy, sigmaxyReceivers );
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressxz, stressxz, sigmaxzReceivers );
compute2dVariableAllSeismoTraces( regionIndex, time_n, dt, stressyz, stressyz, sigmayzReceivers );
-
-
} );
FieldIdentifiers fieldsToBeSync;
@@ -659,22 +642,18 @@ real64 ElasticFirstOrderWaveEquationSEM::explicitStepInternal( real64 const & ti
true );
// compute the seismic traces since last step.
- arrayView2d< real32 > const uxReceivers = m_displacementxNp1AtReceivers.toView();
- arrayView2d< real32 > const uyReceivers = m_displacementyNp1AtReceivers.toView();
- arrayView2d< real32 > const uzReceivers = m_displacementzNp1AtReceivers.toView();
+ arrayView2d< real32 > const uxReceivers = m_displacementxNp1AtReceivers.toView();
+ arrayView2d< real32 > const uyReceivers = m_displacementyNp1AtReceivers.toView();
+ arrayView2d< real32 > const uzReceivers = m_displacementzNp1AtReceivers.toView();
computeAllSeismoTraces( time_n, dt, ux_np1, ux_np1, uxReceivers );
computeAllSeismoTraces( time_n, dt, uy_np1, uy_np1, uyReceivers );
computeAllSeismoTraces( time_n, dt, uz_np1, uz_np1, uzReceivers );
- // increment m_indexSeismoTrace
- while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
- {
- m_indexSeismoTrace++;
- }
-
} );
+ incrementIndexSeismoTrace( time_n );
+
return dt;
}
@@ -704,75 +683,41 @@ void ElasticFirstOrderWaveEquationSEM::cleanup( real64 const time_n,
arrayView2d< real32 const > const stressxz = elementSubRegion.getField< wavesolverfields::Stresstensorxz >();
arrayView2d< real32 const > const stressyz = elementSubRegion.getField< wavesolverfields::Stresstensoryz >();
- arrayView2d< real32 > const sigmaxxReceivers = m_sigmaxxNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmayyReceivers = m_sigmayyNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmazzReceivers = m_sigmazzNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmaxyReceivers = m_sigmaxyNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmaxzReceivers = m_sigmaxzNp1AtReceivers.toView();
- arrayView2d< real32 > const sigmayzReceivers = m_sigmayzNp1AtReceivers.toView();
-
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressxx, stressxx, sigmaxxReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressyy, stressyy, sigmayyReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stresszz, stresszz, sigmazzReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressxy, stressxy, sigmaxyReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressxz, stressxz, sigmaxzReceivers );
- compute2dVariableAllSeismoTraces( regionIndex, time_n, 0, stressyz, stressyz, sigmayzReceivers );
+ arrayView2d< real32 > const sigmaxxReceivers = m_sigmaxxNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmayyReceivers = m_sigmayyNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmazzReceivers = m_sigmazzNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmaxyReceivers = m_sigmaxyNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmaxzReceivers = m_sigmaxzNp1AtReceivers.toView();
+ arrayView2d< real32 > const sigmayzReceivers = m_sigmayzNp1AtReceivers.toView();
+
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressxx, stressxx, sigmaxxReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressyy, stressyy, sigmayyReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stresszz, stresszz, sigmazzReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressxy, stressxy, sigmaxyReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressxz, stressxz, sigmaxzReceivers );
+ compute2dVariableAllSeismoTraces( regionIndex, time_n, 0.0, stressyz, stressyz, sigmayzReceivers );
+
+ WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, sigmaxxReceivers, sigmayyReceivers, sigmazzReceivers );
+ WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, sigmaxyReceivers, sigmaxzReceivers, sigmayzReceivers );
} );
arrayView1d< real32 > const ux_np1 = nodeManager.getField< wavesolverfields::Displacementx_np1 >();
arrayView1d< real32 > const uy_np1 = nodeManager.getField< wavesolverfields::Displacementy_np1 >();
arrayView1d< real32 > const uz_np1 = nodeManager.getField< wavesolverfields::Displacementz_np1 >();
// compute the seismic traces since last step.
- arrayView2d< real32 > const uxReceivers = m_displacementxNp1AtReceivers.toView();
- arrayView2d< real32 > const uyReceivers = m_displacementyNp1AtReceivers.toView();
- arrayView2d< real32 > const uzReceivers = m_displacementzNp1AtReceivers.toView();
+ arrayView2d< real32 > const uxReceivers = m_displacementxNp1AtReceivers.toView();
+ arrayView2d< real32 > const uyReceivers = m_displacementyNp1AtReceivers.toView();
+ arrayView2d< real32 > const uzReceivers = m_displacementzNp1AtReceivers.toView();
- computeAllSeismoTraces( time_n, 0, ux_np1, ux_np1, uxReceivers );
- computeAllSeismoTraces( time_n, 0, uy_np1, uy_np1, uyReceivers );
- computeAllSeismoTraces( time_n, 0, uz_np1, uz_np1, uzReceivers );
+ computeAllSeismoTraces( time_n, 0.0, ux_np1, ux_np1, uxReceivers );
+ computeAllSeismoTraces( time_n, 0.0, uy_np1, uy_np1, uyReceivers );
+ computeAllSeismoTraces( time_n, 0.0, uz_np1, uz_np1, uzReceivers );
+ WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, uxReceivers, uyReceivers, uzReceivers );
} );
-
-// increment m_indexSeismoTrace
- while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
- {
- m_indexSeismoTrace++;
- }
-
-}
-
-void ElasticFirstOrderWaveEquationSEM::computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers )
-{
-
- localIndex indexSeismoTrace = m_indexSeismoTrace;
- for( real64 timeSeismo;
- (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
- indexSeismoTrace++ )
- {
- WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal,
- m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers );
- }
-}
-
-void ElasticFirstOrderWaveEquationSEM::compute2dVariableAllSeismoTraces( localIndex const regionIndex,
- real64 const time_n,
- real64 const dt,
- arrayView2d< real32 const > const var_np1,
- arrayView2d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers )
-{
- localIndex indexSeismoTrace = m_indexSeismoTrace;
- for( real64 timeSeismo;
- (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
- indexSeismoTrace++ )
- {
- WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dt, regionIndex, m_receiverRegion, timeSeismo, indexSeismoTrace, m_rcvElem, m_receiverConstants, m_receiverIsLocal,
- m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers );
- }
}
void ElasticFirstOrderWaveEquationSEM::initializePML()
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp
index a4c3c6a3266..a45cb94d0a3 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEM.hpp
@@ -21,10 +21,8 @@
#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_WAVEPROPAGATION_ELASTICFIRSTORDERWAVEEQUATIONSEM_HPP_
#include "mesh/MeshFields.hpp"
-#include "WaveSolverUtils.hpp"
#include "WaveSolverBaseFields.hpp"
-
-
+#include "WaveSolverBase.hpp"
namespace geos
{
@@ -77,36 +75,6 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase
*/
void addSourceToRightHandSide( integer const & cycleNumber, arrayView1d< real32 > const rhs );
- /**
- * TODO: move implementation into WaveSolverBase
- * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
- * @param time_n the time corresponding to the field values pressure_n
- * @param dt the simulation timestep
- * @param var_np1 the field values at time_n + dt
- * @param var_n the field values at time_n
- * @param varAtreceivers the array holding the trace values, where the output is written
- */
- virtual void computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers );
-
- /**
- * TODO: move implementation into WaveSolverBase
- * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
- * @param time_n the time corresponding to the field values pressure_n
- * @param dt the simulation timestep
- * @param var_np1 the field values at time_n + dt
- * @param var_n the field values at time_n
- * @param varAtreceivers the array holding the trace values, where the output is written
- */
- virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex,
- real64 const time_n,
- real64 const dt,
- arrayView2d< real32 const > const var_np1,
- arrayView2d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers );
/**
* @brief Initialize Perfectly Matched Layer (PML) information
@@ -135,8 +103,6 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase
static constexpr char const * sourceElemString() { return "sourceElem"; }
static constexpr char const * sourceRegionString() { return "sourceRegion"; }
- static constexpr char const * receiverElemString() { return "rcvElem"; }
- static constexpr char const * receiverRegionString() { return "receiverRegion"; }
} waveEquationViewKeys;
@@ -216,12 +182,6 @@ class ElasticFirstOrderWaveEquationSEM : public WaveSolverBase
/// Array containing the elements which contain the region which the source belongs
array1d< localIndex > m_sourceRegion;
- /// Array containing the elements which contain a receiver
- array1d< localIndex > m_rcvElem;
-
- /// Array containing the elements which contain the region which the receiver belongs
- array1d< localIndex > m_receiverRegion;
-
};
} /* namespace geos */
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp
index c4564bbf791..dbb2ba45057 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticFirstOrderWaveEquationSEMKernel.hpp
@@ -50,12 +50,11 @@ struct PrecomputeSourceAndReceiverKernel
* @param[in] receiverCoordinates coordinates of the receiver terms
* @param[in] dt time-step
* @param[in] timeSourceFrequency Peak frequency of the source
+ * @param[in] timeSourceDelay the time delay of the source
* @param[in] rickerOrder Order of the Ricker wavelet
* @param[out] sourceIsLocal flag indicating whether the source is local or not
* @param[out] sourceNodeIds indices of the nodes of the element where the source is located
- * @param[out] sourceConstantsx constant part of the source terms in x-direction
- * @param[out] sourceConstantsy constant part of the source terms in y-direction
- * @param[out] sourceConstantsz constant part of the source terms in z-direction
+ * @param[out] sourceConstants constant part of the source terms in x-direction
* @param[out] receiverIsLocal flag indicating whether the receiver is local or not
* @param[out] receiverNodeIds indices of the nodes of the element where the receiver is located
* @param[out] receiverNodeConstants constant part of the receiver term
@@ -89,6 +88,7 @@ struct PrecomputeSourceAndReceiverKernel
arrayView2d< real32 > const sourceValue,
real64 const dt,
real32 const timeSourceFrequency,
+ real32 const timeSourceDelay,
localIndex const rickerOrder )
{
@@ -133,14 +133,13 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- sourceNodeIds[isrc][a] = elemsToNodes[k][a];
+ sourceNodeIds[isrc][a] = elemsToNodes( k, a );
sourceConstants[isrc][a] = Ntest[a];
}
for( localIndex cycle = 0; cycle < sourceValue.size( 0 ); ++cycle )
{
- real64 const time = cycle*dt;
- sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( time, timeSourceFrequency, rickerOrder );
+ sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( cycle * dt, timeSourceFrequency, timeSourceDelay, rickerOrder );
}
}
@@ -182,7 +181,7 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- receiverNodeIds[ircv][a] = elemsToNodes[k][a];
+ receiverNodeIds[ircv][a] = elemsToNodes( k, a );
receiverConstants[ircv][a] = Ntest[a];
}
}
@@ -239,7 +238,7 @@ struct MassMatrixKernel
for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
{
real32 const localIncrement = density[k] * m_finiteElement.computeMassTerm( q, xLocal );
- RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( k, q )], localIncrement );
}
} ); // end loop over element
}
@@ -262,8 +261,8 @@ struct DampingMatrixKernel
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
- * @param[in] X coordinates of the nodes
- * @param[in] facesToElems map from faces to elements
+ * @param[in] nodeCoords coordinates of the nodes
+ * @param[in] elemsToFaces map from elements to faces
* @param[in] facesToNodes map from face to nodes
* @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise
* @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise
@@ -273,8 +272,8 @@ struct DampingMatrixKernel
template< typename EXEC_POLICY, typename ATOMIC_POLICY >
void
launch( localIndex const size,
- arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
- arrayView2d< localIndex const > const facesToElems,
+ arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
+ arrayView2d< localIndex const > const elemsToFaces,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView1d< integer const > const facesDomainBoundaryIndicator,
arrayView1d< localIndex const > const freeSurfaceFaceIndicator,
@@ -286,46 +285,37 @@ struct DampingMatrixKernel
arrayView1d< real32 > const dampingy,
arrayView1d< real32 > const dampingz )
{
- forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
- // face on the domain boundary and not on free surface
- if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
+ for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i )
{
- localIndex k = facesToElems( f, 0 );
- if( k < 0 )
+ localIndex const f = elemsToFaces( e, i );
+ // face on the domain boundary and not on free surface
+ if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
{
- k = facesToElems( f, 1 );
- }
-
- constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
+ constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
- real64 xLocal[ numNodesPerFace ][ 3 ];
- for( localIndex a = 0; a < numNodesPerFace; ++a )
- {
- for( localIndex i = 0; i < 3; ++i )
+ real64 xLocal[ numNodesPerFace ][ 3 ];
+ for( localIndex a = 0; a < numNodesPerFace; ++a )
{
- xLocal[a][i] = X( facesToNodes( f, a ), i );
+ for( localIndex d = 0; d < 3; ++d )
+ {
+ xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d );
+ }
}
- }
- for( localIndex q = 0; q < numNodesPerFace; ++q )
- {
- real32 const localIncrementx = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][0] ) + velocityVs[k]*sqrt( faceNormal[f][1]*faceNormal[f][1] +
- faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm(
- q,
- xLocal );
- real32 const localIncrementy = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][1] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] +
- faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm(
- q,
- xLocal );
- real32 const localIncrementz = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][2] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] +
- faceNormal[f][1]*faceNormal[f][1] ) )* m_finiteElement.computeDampingTerm(
- q,
- xLocal );
-
- RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes[f][q]], localIncrementx );
- RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes[f][q]], localIncrementy );
- RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes[f][q]], localIncrementz );
+ real32 const nx = faceNormal( f, 0 ), ny = faceNormal( f, 1 ), nz = faceNormal( f, 2 );
+ for( localIndex q = 0; q < numNodesPerFace; ++q )
+ {
+ real32 const aux = density[e] * m_finiteElement.computeDampingTerm( q, xLocal );
+ real32 const localIncrementx = density[e] * (velocityVp[e] * abs( nx ) + velocityVs[e] * sqrt( pow( ny, 2 ) + pow( nz, 2 ) ) ) * aux;
+ real32 const localIncrementy = density[e] * (velocityVp[e] * abs( ny ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( nz, 2 ) ) ) * aux;
+ real32 const localIncrementz = density[e] * (velocityVp[e] * abs( nz ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( ny, 2 ) ) ) * aux;
+
+ RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes( f, q )], localIncrementx );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes( f, q )], localIncrementy );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes( f, q )], localIncrementz );
+ }
}
}
} ); // end loop over element
@@ -425,34 +415,34 @@ struct StressComputation
//Volume integral
m_finiteElement.template computeFirstOrderStiffnessTermX( q, xLocal, [&] ( int i, int j, real32 dfx1, real32 dfx2, real32 dfx3 )
{
- auxx[j]+= dfx1*ux_np1[elemsToNodes[k][i]];
- auyy[j]+= dfx2*uy_np1[elemsToNodes[k][i]];
- auzz[j]+= dfx3*uz_np1[elemsToNodes[k][i]];
- auxy[j]+= dfx1*uy_np1[elemsToNodes[k][i]]+dfx2*ux_np1[elemsToNodes[k][i]];
- auxz[j]+= dfx1*uz_np1[elemsToNodes[k][i]]+dfx3*ux_np1[elemsToNodes[k][i]];
- auyz[j]+= dfx2*uz_np1[elemsToNodes[k][i]]+dfx3*uy_np1[elemsToNodes[k][i]];
+ auxx[j]+= dfx1*ux_np1[elemsToNodes( k, i )];
+ auyy[j]+= dfx2*uy_np1[elemsToNodes( k, i )];
+ auzz[j]+= dfx3*uz_np1[elemsToNodes( k, i )];
+ auxy[j]+= dfx1*uy_np1[elemsToNodes( k, i )]+dfx2*ux_np1[elemsToNodes( k, i )];
+ auxz[j]+= dfx1*uz_np1[elemsToNodes( k, i )]+dfx3*ux_np1[elemsToNodes( k, i )];
+ auyz[j]+= dfx2*uz_np1[elemsToNodes( k, i )]+dfx3*uy_np1[elemsToNodes( k, i )];
} );
m_finiteElement.template computeFirstOrderStiffnessTermY( q, xLocal, [&] ( int i, int j, real32 dfy1, real32 dfy2, real32 dfy3 )
{
- auxx[j]+= dfy1*ux_np1[elemsToNodes[k][i]];
- auyy[j]+= dfy2*uy_np1[elemsToNodes[k][i]];
- auzz[j]+= dfy3*uz_np1[elemsToNodes[k][i]];
- auxy[j]+= dfy1*uy_np1[elemsToNodes[k][i]]+dfy2*ux_np1[elemsToNodes[k][i]];
- auxz[j]+= dfy1*uz_np1[elemsToNodes[k][i]]+dfy3*ux_np1[elemsToNodes[k][i]];
- auyz[j]+= dfy2*uz_np1[elemsToNodes[k][i]]+dfy3*uy_np1[elemsToNodes[k][i]];
+ auxx[j]+= dfy1*ux_np1[elemsToNodes( k, i )];
+ auyy[j]+= dfy2*uy_np1[elemsToNodes( k, i )];
+ auzz[j]+= dfy3*uz_np1[elemsToNodes( k, i )];
+ auxy[j]+= dfy1*uy_np1[elemsToNodes( k, i )]+dfy2*ux_np1[elemsToNodes( k, i )];
+ auxz[j]+= dfy1*uz_np1[elemsToNodes( k, i )]+dfy3*ux_np1[elemsToNodes( k, i )];
+ auyz[j]+= dfy2*uz_np1[elemsToNodes( k, i )]+dfy3*uy_np1[elemsToNodes( k, i )];
} );
m_finiteElement.template computeFirstOrderStiffnessTermZ( q, xLocal, [&] ( int i, int j, real32 dfz1, real32 dfz2, real32 dfz3 )
{
- auxx[j]+= dfz1*ux_np1[elemsToNodes[k][i]];
- auyy[j]+= dfz2*uy_np1[elemsToNodes[k][i]];
- auzz[j]+= dfz3*uz_np1[elemsToNodes[k][i]];
- auxy[j]+= dfz1*uy_np1[elemsToNodes[k][i]]+dfz2*ux_np1[elemsToNodes[k][i]];
- auxz[j]+= dfz1*uz_np1[elemsToNodes[k][i]]+dfz3*ux_np1[elemsToNodes[k][i]];
- auyz[j]+= dfz2*uz_np1[elemsToNodes[k][i]]+dfz3*uy_np1[elemsToNodes[k][i]];
+ auxx[j]+= dfz1*ux_np1[elemsToNodes( k, i )];
+ auyy[j]+= dfz2*uy_np1[elemsToNodes( k, i )];
+ auzz[j]+= dfz3*uz_np1[elemsToNodes( k, i )];
+ auxy[j]+= dfz1*uy_np1[elemsToNodes( k, i )]+dfz2*ux_np1[elemsToNodes( k, i )];
+ auxz[j]+= dfz1*uz_np1[elemsToNodes( k, i )]+dfz3*ux_np1[elemsToNodes( k, i )];
+ auyz[j]+= dfz2*uz_np1[elemsToNodes( k, i )]+dfz3*uy_np1[elemsToNodes( k, i )];
} );
@@ -613,12 +603,12 @@ struct VelocityComputation
// Mult by inverse mass matrix + damping matrix
for( localIndex i = 0; i < numNodesPerElem; ++i )
{
- real32 localIncrement1 = uelemx[i]/mass[elemsToNodes[k][i]];
- real32 localIncrement2 = uelemy[i]/mass[elemsToNodes[k][i]];
- real32 localIncrement3 = uelemz[i]/mass[elemsToNodes[k][i]];
- RAJA::atomicAdd< ATOMIC_POLICY >( &ux_np1[elemsToNodes[k][i]], localIncrement1 );
- RAJA::atomicAdd< ATOMIC_POLICY >( &uy_np1[elemsToNodes[k][i]], localIncrement2 );
- RAJA::atomicAdd< ATOMIC_POLICY >( &uz_np1[elemsToNodes[k][i]], localIncrement3 );
+ real32 localIncrement1 = uelemx[i]/mass[elemsToNodes( k, i )];
+ real32 localIncrement2 = uelemy[i]/mass[elemsToNodes( k, i )];
+ real32 localIncrement3 = uelemz[i]/mass[elemsToNodes( k, i )];
+ RAJA::atomicAdd< ATOMIC_POLICY >( &ux_np1[elemsToNodes( k, i )], localIncrement1 );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &uy_np1[elemsToNodes( k, i )], localIncrement2 );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &uz_np1[elemsToNodes( k, i )], localIncrement3 );
}
} );
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp
index 3416e82f7ef..fa5ca666925 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.cpp
@@ -64,16 +64,6 @@ ElasticWaveEquationSEM::ElasticWaveEquationSEM( const std::string & name,
setSizedFromParent( 0 ).
setDescription( "Flag that indicates whether the source is accessible to this MPI rank" );
- registerWrapper( viewKeyStruct::receiverNodeIdsString(), &m_receiverNodeIds ).
- setInputFlag( InputFlags::FALSE ).
- setSizedFromParent( 0 ).
- setDescription( "Indices of the nodes (in the right order) for each receiver point" );
-
- registerWrapper( viewKeyStruct::sourceConstantsString(), &m_receiverConstants ).
- setInputFlag( InputFlags::FALSE ).
- setSizedFromParent( 0 ).
- setDescription( "Constant part of the receiver for the nodes listed in m_receiverNodeIds" );
-
registerWrapper( viewKeyStruct::receiverIsLocalString(), &m_receiverIsLocal ).
setInputFlag( InputFlags::FALSE ).
setSizedFromParent( 0 ).
@@ -116,7 +106,7 @@ ElasticWaveEquationSEM::~ElasticWaveEquationSEM()
localIndex ElasticWaveEquationSEM::getNumNodesPerElem()
{
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
@@ -199,25 +189,25 @@ void ElasticWaveEquationSEM::registerDataOnMesh( Group & meshBodies )
fields::ForcingRHSx,
fields::ForcingRHSy,
fields::ForcingRHSz,
- fields::MassVector,
+ fields::MassVectorE,
fields::DampingVectorx,
fields::DampingVectory,
fields::DampingVectorz,
fields::StiffnessVectorx,
fields::StiffnessVectory,
fields::StiffnessVectorz,
- fields::FreeSurfaceNodeIndicator >( this->getName() );
+ fields::FreeSurfaceNodeIndicatorE >( getName() );
FaceManager & faceManager = mesh.getFaceManager();
- faceManager.registerField< fields::FreeSurfaceFaceIndicator >( this->getName() );
+ faceManager.registerField< fields::FreeSurfaceFaceIndicatorE >( getName() );
ElementRegionManager & elemManager = mesh.getElemManager();
elemManager.forElementSubRegions< CellElementSubRegion >( [&]( CellElementSubRegion & subRegion )
{
- subRegion.registerField< fields::MediumVelocityVp >( this->getName() );
- subRegion.registerField< fields::MediumVelocityVs >( this->getName() );
- subRegion.registerField< fields::MediumDensity >( this->getName() );
+ subRegion.registerField< fields::MediumVelocityVp >( getName() );
+ subRegion.registerField< fields::MediumVelocityVs >( getName() );
+ subRegion.registerField< fields::MediumDensityE >( getName() );
} );
} );
@@ -230,19 +220,19 @@ void ElasticWaveEquationSEM::postProcessInput()
WaveSolverBase::postProcessInput();
- GEOS_ERROR_IF( m_sourceCoordinates.size( 1 ) != 3,
+ GEOS_ERROR_IF( m_sourceCoordinates.size( 0 ) > 0 && m_sourceCoordinates.size( 1 ) != 3,
"Invalid number of physical coordinates for the sources" );
- GEOS_ERROR_IF( m_receiverCoordinates.size( 1 ) != 3,
+ GEOS_ERROR_IF( m_receiverCoordinates.size( 0 ) > 0 && m_receiverCoordinates.size( 1 ) != 3,
"Invalid number of physical coordinates for the receivers" );
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() );
real64 dt = 0;
for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
{
EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
- if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+ if( subEvent->getEventName() == "/Solvers/" + getName() )
{
dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
}
@@ -266,9 +256,9 @@ void ElasticWaveEquationSEM::postProcessInput()
localIndex const numReceiversGlobal = m_receiverCoordinates.size( 0 );
m_receiverIsLocal.resize( numReceiversGlobal );
- m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
- m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
+ m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
+ m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
m_sourceValue.resize( nsamples, numSourcesGlobal );
/// The receivers are initialized to zero.
@@ -310,16 +300,14 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
receiverConstants.setValues< EXEC_POLICY >( 0 );
receiverIsLocal.zero();
- real32 const timeSourceFrequency = this->m_timeSourceFrequency;
- localIndex const rickerOrder = this->m_rickerOrder;
arrayView2d< real32 > const sourceValue = m_sourceValue.toView();
real64 dt = 0;
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
{
EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
- if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+ if( subEvent->getEventName() == "/Solvers/" + getName() )
{
dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
}
@@ -370,17 +358,18 @@ void ElasticWaveEquationSEM::precomputeSourceAndReceiverTerm( MeshLevel & mesh,
receiverConstants,
sourceValue,
dt,
- timeSourceFrequency,
- rickerOrder,
+ m_timeSourceFrequency,
+ m_timeSourceDelay,
+ m_rickerOrder,
m_sourceForce,
m_sourceMoment );
} );
} );
}
-void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv,
- arrayView2d< real32 > const yCompRcv,
- arrayView2d< real32 > const zCompRcv )
+void ElasticWaveEquationSEM::computeDAS( arrayView2d< real32 > const xCompRcv,
+ arrayView2d< real32 > const yCompRcv,
+ arrayView2d< real32 > const zCompRcv )
{
arrayView2d< real64 const > const linearDASGeometry = m_linearDASGeometry.toViewConst();
@@ -434,23 +423,6 @@ void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv,
} );
}
- /// temporary output to txt
- if( this->m_outputSeismoTrace == 1 )
- {
- forAll< serialPolicy >( numReceiversGlobal, [=] ( localIndex const ircv )
- {
- if( receiverIsLocal[ircv] == 1 )
- {
- std::ofstream f( GEOS_FMT( "dasTraceReceiver{:03}.txt", ircv ), std::ios::app );
- for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample )
- {
- f<< iSample << " " << zCompRcv[iSample][ircv] << std::endl;
- }
- f.close();
- }
- } );
- }
-
/// resize the receiver arrays by dropping the extra pair to avoid confusion
/// the remaining x-component contains DAS data, the other components are set to zero
m_displacementXNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
@@ -468,9 +440,9 @@ void ElasticWaveEquationSEM::computeDAS ( arrayView2d< real32 > const xCompRcv,
}
} );
/// set the y and z components to zero to avoid any confusion
- m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
+ m_displacementYNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
m_displacementYNp1AtReceivers.zero();
- m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal );
+ m_displacementZNp1AtReceivers.resize( m_nsamplesSeismoTrace, numReceiversGlobal + 1 );
m_displacementZNp1AtReceivers.zero();
}
@@ -507,13 +479,11 @@ void ElasticWaveEquationSEM::addSourceToRightHandSide( integer const & cycleNumb
void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
{
-
WaveSolverBase::initializePostInitialConditionsPreSubGroups();
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
- real64 const time = 0.0;
- applyFreeSurfaceBC( time, domain );
+ applyFreeSurfaceBC( 0.0, domain );
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
MeshLevel & mesh,
@@ -523,17 +493,12 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
NodeManager & nodeManager = mesh.getNodeManager();
FaceManager & faceManager = mesh.getFaceManager();
+ ElementRegionManager & elemManager = mesh.getElemManager();
- /// get the array of indicators: 1 if the face is on the boundary; 0 otherwise
- arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator();
- arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.getField< fields::referencePosition32 >().toViewConst();
- arrayView2d< real64 const > const faceNormal = faceManager.faceNormal();
-
- /// get face to node map
- ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
+ arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords = nodeManager.getField< fields::referencePosition32 >().toViewConst();
// mass matrix to be computed in this function
- arrayView1d< real32 > const mass = nodeManager.getField< fields::MassVector >();
+ arrayView1d< real32 > const mass = nodeManager.getField< fields::MassVectorE >();
mass.zero();
/// damping matrix to be computed for each dof in the boundary of the mesh
arrayView1d< real32 > const dampingx = nodeManager.getField< fields::DampingVectorx >();
@@ -544,40 +509,42 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
dampingz.zero();
/// get array of indicators: 1 if face is on the free surface; 0 otherwise
- arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >();
+ arrayView1d< localIndex const > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicatorE >();
+ arrayView1d< integer const > const & facesDomainBoundaryIndicator = faceManager.getDomainBoundaryIndicator();
+ ArrayOfArraysView< localIndex const > const facesToNodes = faceManager.nodeList().toViewConst();
+ arrayView2d< real64 const > const faceNormal = faceManager.faceNormal();
+
- mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
- CellElementSubRegion & elementSubRegion )
+ elemManager.forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
+ CellElementSubRegion & elementSubRegion )
{
+ finiteElement::FiniteElementBase const &
+ fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
arrayView2d< localIndex const, cells::NODE_MAP_USD > const elemsToNodes = elementSubRegion.nodeList();
- arrayView2d< localIndex const > const facesToElements = faceManager.elementList();
- arrayView1d< real32 > const density = elementSubRegion.getField< fields::MediumDensity >();
- arrayView1d< real32 > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >();
- arrayView1d< real32 > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >();
+ arrayView2d< localIndex const > const elemsToFaces = elementSubRegion.faceList();
- finiteElement::FiniteElementBase const &
- fe = elementSubRegion.getReference< finiteElement::FiniteElementBase >( getDiscretizationName() );
+ computeTargetNodeSet( elemsToNodes, elementSubRegion.size(), fe.getNumQuadraturePoints() );
+
+ arrayView1d< real32 const > const density = elementSubRegion.getField< fields::MediumDensityE >();
+ arrayView1d< real32 const > const velocityVp = elementSubRegion.getField< fields::MediumVelocityVp >();
+ arrayView1d< real32 const > const velocityVs = elementSubRegion.getField< fields::MediumVelocityVs >();
- finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe,
- [&]
- ( auto const finiteElement )
+ finiteElement::FiniteElementDispatchHandler< SEM_FE_TYPES >::dispatch3D( fe, [&] ( auto const finiteElement )
{
using FE_TYPE = TYPEOFREF( finiteElement );
elasticWaveEquationSEMKernels::MassMatrixKernel< FE_TYPE > kernelM( finiteElement );
-
kernelM.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
- X,
+ nodeCoords,
elemsToNodes,
density,
mass );
elasticWaveEquationSEMKernels::DampingMatrixKernel< FE_TYPE > kernelD( finiteElement );
-
- kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( faceManager.size(),
- X,
- facesToElements,
+ kernelD.template launch< EXEC_POLICY, ATOMIC_POLICY >( elementSubRegion.size(),
+ nodeCoords,
+ elemsToFaces,
facesToNodes,
facesDomainBoundaryIndicator,
freeSurfaceFaceIndicator,
@@ -592,9 +559,10 @@ void ElasticWaveEquationSEM::initializePostInitialConditionsPreSubGroups()
} );
} );
+ WaveSolverUtils::initTrace( "seismoTraceReceiver", getName(), m_receiverConstants.size( 0 ), m_receiverIsLocal );
+ WaveSolverUtils::initTrace( "dasTraceReceiver", getName(), m_linearDASGeometry.size( 0 ), m_receiverIsLocal );
}
-
void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartition & domain )
{
FieldSpecificationManager & fsManager = FieldSpecificationManager::getInstance();
@@ -606,9 +574,9 @@ void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartit
arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
- arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >();
- arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >();
- arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >();
+ arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >();
+ arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >();
+ arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >();
arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >();
arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >();
arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >();
@@ -616,10 +584,10 @@ void ElasticWaveEquationSEM::applyFreeSurfaceBC( real64 const time, DomainPartit
ArrayOfArraysView< localIndex const > const faceToNodeMap = faceManager.nodeList().toViewConst();
/// set array of indicators: 1 if a face is on on free surface; 0 otherwise
- arrayView1d< localIndex > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicator >();
+ arrayView1d< localIndex > const freeSurfaceFaceIndicator = faceManager.getField< fields::FreeSurfaceFaceIndicatorE >();
/// set array of indicators: 1 if a node is on on free surface; 0 otherwise
- arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicator >();
+ arrayView1d< localIndex > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicatorE >();
fsManager.apply( time,
@@ -675,8 +643,7 @@ real64 ElasticWaveEquationSEM::explicitStepForward( real64 const & time_n,
DomainPartition & domain,
bool GEOS_UNUSED_PARAM( computeGradient ) )
{
- real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain );
- return dtOut;
+ return explicitStepInternal( time_n, dt, cycleNumber, domain );
}
@@ -688,137 +655,192 @@ real64 ElasticWaveEquationSEM::explicitStepBackward( real64 const & time_n,
bool GEOS_UNUSED_PARAM( computeGradient ) )
{
GEOS_ERROR( "Backward propagation for the elastic wave propagator not yet implemented" );
- real64 dtOut = explicitStepInternal( time_n, dt, cycleNumber, domain );
- return dtOut;
+ return explicitStepInternal( time_n, dt, cycleNumber, domain );
}
-real64 ElasticWaveEquationSEM::explicitStepInternal( real64 const & time_n,
- real64 const & dt,
- integer const cycleNumber,
- DomainPartition & domain )
+void ElasticWaveEquationSEM::computeUnknowns( real64 const &,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
{
- GEOS_MARK_FUNCTION;
+ NodeManager & nodeManager = mesh.getNodeManager();
- GEOS_UNUSED_VAR( time_n, dt, cycleNumber );
+ arrayView1d< real32 const > const mass = nodeManager.getField< fields::MassVectorE >();
+ arrayView1d< real32 const > const dampingx = nodeManager.getField< fields::DampingVectorx >();
+ arrayView1d< real32 const > const dampingy = nodeManager.getField< fields::DampingVectory >();
+ arrayView1d< real32 const > const dampingz = nodeManager.getField< fields::DampingVectorz >();
+ arrayView1d< real32 > const stiffnessVectorx = nodeManager.getField< fields::StiffnessVectorx >();
+ arrayView1d< real32 > const stiffnessVectory = nodeManager.getField< fields::StiffnessVectory >();
+ arrayView1d< real32 > const stiffnessVectorz = nodeManager.getField< fields::StiffnessVectorz >();
- GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc );
- forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
- MeshLevel & mesh,
- arrayView1d< string const > const & regionNames )
+ arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >();
+ arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >();
+ arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >();
+ arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >();
+ arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >();
+ arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >();
+ arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
+ arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
+ arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
+
+ /// get array of indicators: 1 if node on free surface; 0 otherwise
+ arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicatorE >();
+
+ arrayView1d< real32 > const rhsx = nodeManager.getField< fields::ForcingRHSx >();
+ arrayView1d< real32 > const rhsy = nodeManager.getField< fields::ForcingRHSy >();
+ arrayView1d< real32 > const rhsz = nodeManager.getField< fields::ForcingRHSz >();
+
+ auto kernelFactory = elasticWaveEquationSEMKernels::ExplicitElasticSEMFactory( dt );
+
+ finiteElement::
+ regionBasedKernelApplication< EXEC_POLICY,
+ constitutive::NullModel,
+ CellElementSubRegion >( mesh,
+ regionNames,
+ getDiscretizationName(),
+ "",
+ kernelFactory );
+
+
+ addSourceToRightHandSide( cycleNumber, rhsx, rhsy, rhsz );
+
+ real64 const dt2 = pow( dt, 2 );
+ SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst();
+ forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
{
- NodeManager & nodeManager = mesh.getNodeManager();
+ localIndex const a = solverTargetNodesSet[n];
+ if( freeSurfaceNodeIndicator[a] != 1 )
+ {
+ ux_np1[a] = ux_n[a];
+ ux_np1[a] *= 2.0*mass[a];
+ ux_np1[a] -= (mass[a]-0.5*dt*dampingx[a])*ux_nm1[a];
+ ux_np1[a] += dt2*(rhsx[a]-stiffnessVectorx[a]);
+ ux_np1[a] /= mass[a]+0.5*dt*dampingx[a];
+ uy_np1[a] = uy_n[a];
+ uy_np1[a] *= 2.0*mass[a];
+ uy_np1[a] -= (mass[a]-0.5*dt*dampingy[a])*uy_nm1[a];
+ uy_np1[a] += dt2*(rhsy[a]-stiffnessVectory[a]);
+ uy_np1[a] /= mass[a]+0.5*dt*dampingy[a];
+ uz_np1[a] = uz_n[a];
+ uz_np1[a] *= 2.0*mass[a];
+ uz_np1[a] -= (mass[a]-0.5*dt*dampingz[a])*uz_nm1[a];
+ uz_np1[a] += dt2*(rhsz[a]-stiffnessVectorz[a]);
+ uz_np1[a] /= mass[a]+0.5*dt*dampingz[a];
+ }
+ } );
+}
- arrayView1d< real32 const > const mass = nodeManager.getField< fields::MassVector >();
- arrayView1d< real32 const > const dampingx = nodeManager.getField< fields::DampingVectorx >();
- arrayView1d< real32 const > const dampingy = nodeManager.getField< fields::DampingVectory >();
- arrayView1d< real32 const > const dampingz = nodeManager.getField< fields::DampingVectorz >();
- arrayView1d< real32 > const stiffnessVectorx = nodeManager.getField< fields::StiffnessVectorx >();
- arrayView1d< real32 > const stiffnessVectory = nodeManager.getField< fields::StiffnessVectory >();
- arrayView1d< real32 > const stiffnessVectorz = nodeManager.getField< fields::StiffnessVectorz >();
+void ElasticWaveEquationSEM::synchronizeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer const,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & )
+{
+ NodeManager & nodeManager = mesh.getNodeManager();
+ arrayView1d< real32 > const stiffnessVectorx = nodeManager.getField< fields::StiffnessVectorx >();
+ arrayView1d< real32 > const stiffnessVectory = nodeManager.getField< fields::StiffnessVectory >();
+ arrayView1d< real32 > const stiffnessVectorz = nodeManager.getField< fields::StiffnessVectorz >();
- arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >();
- arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >();
- arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >();
- arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >();
- arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >();
- arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >();
- arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
- arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
- arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
+ arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >();
+ arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >();
+ arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >();
+ arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >();
+ arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >();
+ arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >();
+ arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
+ arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
+ arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
- /// get array of indicators: 1 if node on free surface; 0 otherwise
- arrayView1d< localIndex const > const freeSurfaceNodeIndicator = nodeManager.getField< fields::FreeSurfaceNodeIndicator >();
+ arrayView1d< real32 > const rhsx = nodeManager.getField< fields::ForcingRHSx >();
+ arrayView1d< real32 > const rhsy = nodeManager.getField< fields::ForcingRHSy >();
+ arrayView1d< real32 > const rhsz = nodeManager.getField< fields::ForcingRHSz >();
- arrayView1d< real32 > const rhsx = nodeManager.getField< fields::ForcingRHSx >();
- arrayView1d< real32 > const rhsy = nodeManager.getField< fields::ForcingRHSy >();
- arrayView1d< real32 > const rhsz = nodeManager.getField< fields::ForcingRHSz >();
+ /// synchronize displacement fields
+ FieldIdentifiers fieldsToBeSync;
+ fieldsToBeSync.addFields( FieldLocation::Node, { fields::Displacementx_np1::key(), fields::Displacementy_np1::key(), fields::Displacementz_np1::key() } );
- auto kernelFactory = elasticWaveEquationSEMKernels::ExplicitElasticSEMFactory( dt );
+ CommunicationTools & syncFields = CommunicationTools::getInstance();
+ syncFields.synchronizeFields( fieldsToBeSync,
+ domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ),
+ domain.getNeighbors(),
+ true );
- finiteElement::
- regionBasedKernelApplication< EXEC_POLICY,
- constitutive::NullModel,
- CellElementSubRegion >( mesh,
- regionNames,
- getDiscretizationName(),
- "",
- kernelFactory );
+ // compute the seismic traces since last step.
+ arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView();
+ arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView();
+ arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView();
+ computeAllSeismoTraces( time_n, dt, ux_np1, ux_n, uXReceivers );
+ computeAllSeismoTraces( time_n, dt, uy_np1, uy_n, uYReceivers );
+ computeAllSeismoTraces( time_n, dt, uz_np1, uz_n, uZReceivers );
- addSourceToRightHandSide( cycleNumber, rhsx, rhsy, rhsz );
+ incrementIndexSeismoTrace( time_n );
+}
+void ElasticWaveEquationSEM::prepareNextTimestep( MeshLevel & mesh )
+{
+ NodeManager & nodeManager = mesh.getNodeManager();
+ arrayView1d< real32 > const ux_nm1 = nodeManager.getField< fields::Displacementx_nm1 >();
+ arrayView1d< real32 > const uy_nm1 = nodeManager.getField< fields::Displacementy_nm1 >();
+ arrayView1d< real32 > const uz_nm1 = nodeManager.getField< fields::Displacementz_nm1 >();
+ arrayView1d< real32 > const ux_n = nodeManager.getField< fields::Displacementx_n >();
+ arrayView1d< real32 > const uy_n = nodeManager.getField< fields::Displacementy_n >();
+ arrayView1d< real32 > const uz_n = nodeManager.getField< fields::Displacementz_n >();
+ arrayView1d< real32 > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
+ arrayView1d< real32 > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
+ arrayView1d< real32 > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
- real64 const dt2 = dt*dt;
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
- {
- if( freeSurfaceNodeIndicator[a] != 1 )
- {
- ux_np1[a] = ux_n[a];
- ux_np1[a] *= 2.0*mass[a];
- ux_np1[a] -= (mass[a]-0.5*dt*dampingx[a])*ux_nm1[a];
- ux_np1[a] += dt2*(rhsx[a]-stiffnessVectorx[a]);
- ux_np1[a] /= mass[a]+0.5*dt*dampingx[a];
- uy_np1[a] = uy_n[a];
- uy_np1[a] *= 2.0*mass[a];
- uy_np1[a] -= (mass[a]-0.5*dt*dampingy[a])*uy_nm1[a];
- uy_np1[a] += dt2*(rhsy[a]-stiffnessVectory[a]);
- uy_np1[a] /= mass[a]+0.5*dt*dampingy[a];
- uz_np1[a] = uz_n[a];
- uz_np1[a] *= 2.0*mass[a];
- uz_np1[a] -= (mass[a]-0.5*dt*dampingz[a])*uz_nm1[a];
- uz_np1[a] += dt2*(rhsz[a]-stiffnessVectorz[a]);
- uz_np1[a] /= mass[a]+0.5*dt*dampingz[a];
- }
- } );
+ arrayView1d< real32 > const stiffnessVectorx = nodeManager.getField< fields::StiffnessVectorx >();
+ arrayView1d< real32 > const stiffnessVectory = nodeManager.getField< fields::StiffnessVectory >();
+ arrayView1d< real32 > const stiffnessVectorz = nodeManager.getField< fields::StiffnessVectorz >();
- /// synchronize pressure fields
- FieldIdentifiers fieldsToBeSync;
- fieldsToBeSync.addFields( FieldLocation::Node, { fields::Displacementx_np1::key(), fields::Displacementy_np1::key(), fields::Displacementz_np1::key() } );
+ arrayView1d< real32 > const rhsx = nodeManager.getField< fields::ForcingRHSx >();
+ arrayView1d< real32 > const rhsy = nodeManager.getField< fields::ForcingRHSy >();
+ arrayView1d< real32 > const rhsz = nodeManager.getField< fields::ForcingRHSz >();
- CommunicationTools & syncFields = CommunicationTools::getInstance();
- syncFields.synchronizeFields( fieldsToBeSync,
- domain.getMeshBody( 0 ).getMeshLevel( m_discretizationName ),
- domain.getNeighbors(),
- true );
+ SortedArrayView< localIndex const > const solverTargetNodesSet = m_solverTargetNodesSet.toViewConst();
- // compute the seismic traces since last step.
- arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView();
- arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView();
- arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView();
+ forAll< EXEC_POLICY >( solverTargetNodesSet.size(), [=] GEOS_HOST_DEVICE ( localIndex const n )
+ {
+ localIndex const a = solverTargetNodesSet[n];
+ ux_nm1[a] = ux_n[a];
+ uy_nm1[a] = uy_n[a];
+ uz_nm1[a] = uz_n[a];
+ ux_n[a] = ux_np1[a];
+ uy_n[a] = uy_np1[a];
+ uz_n[a] = uz_np1[a];
+
+ stiffnessVectorx[a] = stiffnessVectory[a] = stiffnessVectorz[a] = 0.0;
+ rhsx[a] = rhsy[a] = rhsz[a] = 0.0;
+ } );
- computeAllSeismoTraces( time_n, dt, ux_np1, ux_n, uXReceivers );
- computeAllSeismoTraces( time_n, dt, uy_np1, uy_n, uYReceivers );
- computeAllSeismoTraces( time_n, dt, uz_np1, uz_n, uZReceivers );
+}
- forAll< EXEC_POLICY >( nodeManager.size(), [=] GEOS_HOST_DEVICE ( localIndex const a )
- {
- ux_nm1[a] = ux_n[a];
- uy_nm1[a] = uy_n[a];
- uz_nm1[a] = uz_n[a];
- ux_n[a] = ux_np1[a];
- uy_n[a] = uy_np1[a];
- uz_n[a] = uz_np1[a];
-
- stiffnessVectorx[a] = 0.0;
- stiffnessVectory[a] = 0.0;
- stiffnessVectorz[a] = 0.0;
- rhsx[a] = 0.0;
- rhsy[a] = 0.0;
- rhsz[a] = 0.0;
- } );
+real64 ElasticWaveEquationSEM::explicitStepInternal( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain )
+{
+ GEOS_MARK_FUNCTION;
- // increment m_indexSeismoTrace
- while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
- {
- m_indexSeismoTrace++;
- }
+ GEOS_LOG_RANK_0_IF( dt < epsilonLoc, "Warning! Value for dt: " << dt << "s is smaller than local threshold: " << epsilonLoc );
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+ {
+ computeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames );
+ synchronizeUnknowns( time_n, dt, cycleNumber, domain, mesh, regionNames );
+ prepareNextTimestep( mesh );
} );
- return dt;
+ return dt;
}
void ElasticWaveEquationSEM::cleanup( real64 const time_n,
@@ -836,50 +858,32 @@ void ElasticWaveEquationSEM::cleanup( real64 const time_n,
arrayView1d< string const > const & )
{
NodeManager & nodeManager = mesh.getNodeManager();
- arrayView1d< real32 const > const ux_n = nodeManager.getField< fields::Displacementx_n >();
+ arrayView1d< real32 const > const ux_n = nodeManager.getField< fields::Displacementx_n >();
arrayView1d< real32 const > const ux_np1 = nodeManager.getField< fields::Displacementx_np1 >();
- arrayView1d< real32 const > const uy_n = nodeManager.getField< fields::Displacementy_n >();
+ arrayView1d< real32 const > const uy_n = nodeManager.getField< fields::Displacementy_n >();
arrayView1d< real32 const > const uy_np1 = nodeManager.getField< fields::Displacementy_np1 >();
- arrayView1d< real32 const > const uz_n = nodeManager.getField< fields::Displacementz_n >();
+ arrayView1d< real32 const > const uz_n = nodeManager.getField< fields::Displacementz_n >();
arrayView1d< real32 const > const uz_np1 = nodeManager.getField< fields::Displacementz_np1 >();
- arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView();
- arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView();
- arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView();
+ arrayView2d< real32 > const uXReceivers = m_displacementXNp1AtReceivers.toView();
+ arrayView2d< real32 > const uYReceivers = m_displacementYNp1AtReceivers.toView();
+ arrayView2d< real32 > const uZReceivers = m_displacementZNp1AtReceivers.toView();
- computeAllSeismoTraces( time_n, 0, ux_np1, ux_n, uXReceivers );
- computeAllSeismoTraces( time_n, 0, uy_np1, uy_n, uYReceivers );
- computeAllSeismoTraces( time_n, 0, uz_np1, uz_n, uZReceivers );
+ computeAllSeismoTraces( time_n, 0.0, ux_np1, ux_n, uXReceivers );
+ computeAllSeismoTraces( time_n, 0.0, uy_np1, uy_n, uYReceivers );
+ computeAllSeismoTraces( time_n, 0.0, uz_np1, uz_n, uZReceivers );
+
+ WaveSolverUtils::writeSeismoTraceVector( "seismoTraceReceiver", getName(), m_outputSeismoTrace, m_receiverConstants.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, uXReceivers, uYReceivers, uZReceivers );
/// Compute DAS data if requested
/// Pairs of receivers are assumed to be modeled ( see WaveSolverBase::initializeDAS() )
if( m_useDAS )
{
computeDAS( uXReceivers, uYReceivers, uZReceivers );
+ WaveSolverUtils::writeSeismoTrace( "dasTraceReceiver", getName(), m_outputSeismoTrace, m_linearDASGeometry.size( 0 ),
+ m_receiverIsLocal, m_nsamplesSeismoTrace, uZReceivers );
}
} );
-
- // increment m_indexSeismoTrace
- while( (m_dtSeismoTrace*m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
- {
- m_indexSeismoTrace++;
- }
-}
-
-void ElasticWaveEquationSEM::computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers )
-{
- localIndex indexSeismoTrace = m_indexSeismoTrace;
-
- for( real64 timeSeismo;
- (timeSeismo = m_dtSeismoTrace*indexSeismoTrace) <= (time_n + epsilonLoc) && indexSeismoTrace < m_nsamplesSeismoTrace;
- indexSeismoTrace++ )
- {
- WaveSolverUtils::computeSeismoTrace( time_n, dt, timeSeismo, indexSeismoTrace, m_receiverNodeIds, m_receiverConstants, m_receiverIsLocal,
- m_nsamplesSeismoTrace, m_outputSeismoTrace, var_np1, var_n, varAtReceivers );
- }
}
void ElasticWaveEquationSEM::initializePML()
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp
index 15f12013b93..1088d2b3490 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEM.hpp
@@ -31,6 +31,9 @@ class ElasticWaveEquationSEM : public WaveSolverBase
{
public:
+ /// String used to form the solverName used to register solvers in CoupledSolver
+ static string coupledSolverAttributePrefix() { return "elastic"; }
+
using EXEC_POLICY = parallelDevicePolicy< >;
using ATOMIC_POLICY = parallelDeviceAtomic;
@@ -89,21 +92,6 @@ class ElasticWaveEquationSEM : public WaveSolverBase
/**
* TODO: move implementation into WaveSolverBase
- * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
- * @param time_n the time corresponding to the field values at iteration n
- * @param dt the simulation timestep
- * @param var_np1 the field values at time_n + dt
- * @param var_n the field values at time_n
- * @param varAtreceivers the array holding the trace values, where the output is written
- */
- virtual void computeAllSeismoTraces( real64 const time_n,
- real64 const dt,
- arrayView1d< real32 const > const var_np1,
- arrayView1d< real32 const > const var_n,
- arrayView2d< real32 > varAtReceivers );
-
- /**
- * TODO: move implementation into WaveSolverBase once 'm_receiverIsLocal' is also moved
* @brief Compute DAS data from the appropriate three-component receiver pairs
* @param xCompRcv the array holding the x-component of pairs of receivers
* @param yCompRcv the array holding the y-component of pairs of receivers
@@ -123,16 +111,8 @@ class ElasticWaveEquationSEM : public WaveSolverBase
real64 const eventProgress,
DomainPartition & domain ) override;
- struct viewKeyStruct : SolverBase::viewKeyStruct
+ struct viewKeyStruct : WaveSolverBase::viewKeyStruct
{
- static constexpr char const * sourceNodeIdsString() { return "sourceNodeIds"; }
- static constexpr char const * sourceConstantsString() { return "sourceConstants"; }
- static constexpr char const * sourceIsAccessibleString() { return "sourceIsAccessible"; }
-
- static constexpr char const * receiverNodeIdsString() { return "receiverNodeIds"; }
- static constexpr char const * receiverConstantsString() {return "receiverConstants"; }
- static constexpr char const * receiverIsLocalString() { return "receiverIsLocal"; }
-
static constexpr char const * displacementXNp1AtReceiversString() { return "displacementXNp1AtReceivers"; }
static constexpr char const * displacementYNp1AtReceiversString() { return "displacementYNp1AtReceivers"; }
static constexpr char const * displacementZNp1AtReceiversString() { return "displacementZNp1AtReceivers"; }
@@ -155,6 +135,23 @@ class ElasticWaveEquationSEM : public WaveSolverBase
real64 const & dt,
integer const cycleNumber,
DomainPartition & domain );
+
+ void computeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames );
+
+ void synchronizeUnknowns( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames );
+
+ void prepareNextTimestep( MeshLevel & mesh );
+
protected:
virtual void postProcessInput() override final;
@@ -198,24 +195,15 @@ class ElasticWaveEquationSEM : public WaveSolverBase
/// Constant part of the source for the nodes listed in m_sourceNodeIds in x-direction
array2d< real64 > m_sourceConstantsx;
- /// Constant part of the source for the nodes listed in m_sourceNodeIds in x-direction
+ /// Constant part of the source for the nodes listed in m_sourceNodeIds in y-direction
array2d< real64 > m_sourceConstantsy;
- /// Constant part of the source for the nodes listed in m_sourceNodeIds in x-direction
+ /// Constant part of the source for the nodes listed in m_sourceNodeIds in z-direction
array2d< real64 > m_sourceConstantsz;
/// Flag that indicates whether the source is accessible or not to the MPI rank
array1d< localIndex > m_sourceIsAccessible;
- /// Indices of the element nodes (in the right order) for each receiver point
- array2d< localIndex > m_receiverNodeIds;
-
- /// Basis function evaluated at the receiver for the nodes listed in m_receiverNodeIds
- array2d< real64 > m_receiverConstants;
-
- /// Flag that indicates whether the receiver is local or not to the MPI rank
- array1d< localIndex > m_receiverIsLocal;
-
/// Displacement_np1 at the receiver location for each time step for each receiver (x-component)
array2d< real32 > m_displacementXNp1AtReceivers;
@@ -333,8 +321,8 @@ DECLARE_FIELD( ForcingRHSz,
WRITE_AND_READ,
"RHS for z-direction" );
-DECLARE_FIELD( MassVector,
- "massVector",
+DECLARE_FIELD( MassVectorE,
+ "massVectorE",
array1d< real32 >,
0,
NOPLOT,
@@ -405,24 +393,24 @@ DECLARE_FIELD( MediumVelocityVs,
WRITE_AND_READ,
"S-waves speed in the cell" );
-DECLARE_FIELD( MediumDensity,
- "mediumDensity",
+DECLARE_FIELD( MediumDensityE,
+ "mediumDensityE",
array1d< real32 >,
0,
NOPLOT,
WRITE_AND_READ,
"Medium density of the cell" );
-DECLARE_FIELD( FreeSurfaceFaceIndicator,
- "freeSurfaceFaceIndicator",
+DECLARE_FIELD( FreeSurfaceFaceIndicatorE,
+ "freeSurfaceFaceIndicatorE",
array1d< localIndex >,
0,
NOPLOT,
WRITE_AND_READ,
"Free surface indicator, 1 if a face is on free surface 0 otherwise." );
-DECLARE_FIELD( FreeSurfaceNodeIndicator,
- "freeSurfaceNodeIndicator",
+DECLARE_FIELD( FreeSurfaceNodeIndicatorE,
+ "freeSurfaceNodeIndicatorE",
array1d< localIndex >,
0,
NOPLOT,
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp
index a776dec8a6c..fe9035d0dc5 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/ElasticWaveEquationSEMKernel.hpp
@@ -32,7 +32,6 @@ namespace elasticWaveEquationSEMKernels
struct PrecomputeSourceAndReceiverKernel
{
-
/**
* @brief Launches the precomputation of the source and receiver terms
* @tparam EXEC_POLICY execution policy
@@ -85,14 +84,13 @@ struct PrecomputeSourceAndReceiverKernel
arrayView2d< real32 > const sourceValue,
real64 const dt,
real32 const timeSourceFrequency,
+ real32 const timeSourceDelay,
localIndex const rickerOrder,
R1Tensor const sourceForce,
R2SymTensor const sourceMoment )
{
-
forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
{
-
constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
real64 const center[3] = { elemCenter[k][0],
@@ -148,7 +146,7 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex q=0; q< numNodesPerElem; ++q )
{
real64 inc[3] = { 0, 0, 0 };
- sourceNodeIds[isrc][q] = elemsToNodes[k][q];
+ sourceNodeIds[isrc][q] = elemsToNodes( k, q );
inc[0] += sourceForce[0] * N[q];
inc[1] += sourceForce[1] * N[q];
inc[2] += sourceForce[2] * N[q];
@@ -161,8 +159,7 @@ struct PrecomputeSourceAndReceiverKernel
for( localIndex cycle = 0; cycle < sourceValue.size( 0 ); ++cycle )
{
- real64 const time = cycle*dt;
- sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( time, timeSourceFrequency, rickerOrder );
+ sourceValue[cycle][isrc] = WaveSolverUtils::evaluateRicker( cycle * dt, timeSourceFrequency, timeSourceDelay, rickerOrder );
}
}
@@ -200,12 +197,11 @@ struct PrecomputeSourceAndReceiverKernel
receiverIsLocal[ircv] = 1;
real64 Ntest[numNodesPerElem];
-
FE_TYPE::calcN( coordsOnRefElem, Ntest );
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
- receiverNodeIds[ircv][a] = elemsToNodes[k][a];
+ receiverNodeIds[ircv][a] = elemsToNodes( k, a );
receiverConstants[ircv][a] = Ntest[a];
}
}
@@ -245,7 +241,7 @@ struct MassMatrixKernel
arrayView1d< real32 > const mass )
{
- forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const k )
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
constexpr localIndex numNodesPerElem = FE_TYPE::numNodes;
@@ -256,14 +252,14 @@ struct MassMatrixKernel
{
for( localIndex i = 0; i < 3; ++i )
{
- xLocal[a][i] = X( elemsToNodes( k, a ), i );
+ xLocal[a][i] = X( elemsToNodes( e, a ), i );
}
}
for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
{
- real32 const localIncrement = density[k] * m_finiteElement.computeMassTerm( q, xLocal );
- RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes[k][q]], localIncrement );
+ real32 const localIncrement = density[e] * m_finiteElement.computeMassTerm( q, xLocal );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &mass[elemsToNodes( e, q )], localIncrement );
}
} ); // end loop over element
}
@@ -276,7 +272,6 @@ struct MassMatrixKernel
template< typename FE_TYPE >
struct DampingMatrixKernel
{
-
DampingMatrixKernel( FE_TYPE const & finiteElement )
: m_finiteElement( finiteElement )
{}
@@ -286,8 +281,8 @@ struct DampingMatrixKernel
* @tparam EXEC_POLICY the execution policy
* @tparam ATOMIC_POLICY the atomic policy
* @param[in] size the number of cells in the subRegion
- * @param[in] X coordinates of the nodes
- * @param[in] facesToElems map from faces to elements
+ * @param[in] nodeCoords coordinates of the nodes
+ * @param[in] elemsToFaces map from elements to faces
* @param[in] facesToNodes map from face to nodes
* @param[in] facesDomainBoundaryIndicator flag equal to 1 if the face is on the boundary, and to 0 otherwise
* @param[in] freeSurfaceFaceIndicator flag equal to 1 if the face is on the free surface, and to 0 otherwise
@@ -298,8 +293,8 @@ struct DampingMatrixKernel
//std::enable_if_t< geos::is_sem_formulation< std::remove_cv_t< FE_TYPE_ > >::value, void >
void
launch( localIndex const size,
- arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
- arrayView2d< localIndex const > const facesToElems,
+ arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const nodeCoords,
+ arrayView2d< localIndex const > const elemsToFaces,
ArrayOfArraysView< localIndex const > const facesToNodes,
arrayView1d< integer const > const facesDomainBoundaryIndicator,
arrayView1d< localIndex const > const freeSurfaceFaceIndicator,
@@ -311,49 +306,40 @@ struct DampingMatrixKernel
arrayView1d< real32 > const dampingy,
arrayView1d< real32 > const dampingz )
{
- forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const f )
+
+ forAll< EXEC_POLICY >( size, [=] GEOS_HOST_DEVICE ( localIndex const e )
{
- // face on the domain boundary and not on free surface
- if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
+ for( localIndex i = 0; i < elemsToFaces.size( 1 ); ++i )
{
- localIndex k = facesToElems( f, 0 );
- if( k < 0 )
- {
- k = facesToElems( f, 1 );
- }
-
- constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
-
- real64 xLocal[ numNodesPerFace ][ 3 ];
- for( localIndex a = 0; a < numNodesPerFace; ++a )
+ localIndex const f = elemsToFaces( e, i );
+ // face on the domain boundary and not on free surface
+ if( facesDomainBoundaryIndicator[f] == 1 && freeSurfaceFaceIndicator[f] != 1 )
{
- for( localIndex i = 0; i < 3; ++i )
+ constexpr localIndex numNodesPerFace = FE_TYPE::numNodesPerFace;
+ real64 xLocal[ numNodesPerFace ][ 3 ];
+ for( localIndex a = 0; a < numNodesPerFace; ++a )
{
- xLocal[a][i] = X( facesToNodes( f, a ), i );
+ for( localIndex d = 0; d < 3; ++d )
+ {
+ xLocal[a][d] = nodeCoords( facesToNodes( f, a ), d );
+ }
}
- }
- for( localIndex q = 0; q < numNodesPerFace; ++q )
- {
- real32 const localIncrementx = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][0] ) + velocityVs[k]*sqrt( faceNormal[f][1]*faceNormal[f][1] +
- faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm(
- q,
- xLocal );
- real32 const localIncrementy = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][1] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] +
- faceNormal[f][2]*faceNormal[f][2] ) )* m_finiteElement.computeDampingTerm(
- q,
- xLocal );
- real32 const localIncrementz = density[k] * (velocityVp[k]*LvArray::math::abs( faceNormal[f][2] ) + velocityVs[k]*sqrt( faceNormal[f][0]*faceNormal[f][0] +
- faceNormal[f][1]*faceNormal[f][1] ) )* m_finiteElement.computeDampingTerm(
- q,
- xLocal );
-
- RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes[f][q]], localIncrementx );
- RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes[f][q]], localIncrementy );
- RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes[f][q]], localIncrementz );
+ real32 const nx = faceNormal( f, 0 ), ny = faceNormal( f, 1 ), nz = faceNormal( f, 2 );
+ for( localIndex q = 0; q < numNodesPerFace; ++q )
+ {
+ real32 const aux = density[e] * m_finiteElement.computeDampingTerm( q, xLocal );
+ real32 const localIncrementx = aux * ( velocityVp[e] * abs( nx ) + velocityVs[e] * sqrt( pow( ny, 2 ) + pow( nz, 2 ) ) );
+ real32 const localIncrementy = aux * ( velocityVp[e] * abs( ny ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( nz, 2 ) ) );
+ real32 const localIncrementz = aux * ( velocityVp[e] * abs( nz ) + velocityVs[e] * sqrt( pow( nx, 2 ) + pow( ny, 2 ) ) );
+
+ RAJA::atomicAdd< ATOMIC_POLICY >( &dampingx[facesToNodes( f, q )], localIncrementx );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &dampingy[facesToNodes( f, q )], localIncrementy );
+ RAJA::atomicAdd< ATOMIC_POLICY >( &dampingz[facesToNodes( f, q )], localIncrementz );
+ }
}
}
- } ); // end loop over element
+ } );
}
/// The finite element space/discretization object for the element type in the subRegion
@@ -435,7 +421,7 @@ class ExplicitElasticSEM : public finiteElement::KernelBase< SUBREGION_TYPE,
m_stiffnessVectorx( nodeManager.getField< fields::StiffnessVectorx >() ),
m_stiffnessVectory( nodeManager.getField< fields::StiffnessVectory >() ),
m_stiffnessVectorz( nodeManager.getField< fields::StiffnessVectorz >() ),
- m_density( elementSubRegion.template getField< fields::MediumDensity >() ),
+ m_density( elementSubRegion.template getField< fields::MediumDensityE >() ),
m_velocityVp( elementSubRegion.template getField< fields::MediumVelocityVp >() ),
m_velocityVs( elementSubRegion.template getField< fields::MediumVelocityVs >() ),
m_dt( dt )
@@ -518,13 +504,13 @@ class ExplicitElasticSEM : public finiteElement::KernelBase< SUBREGION_TYPE,
real32 const Ryz_ij = val*(stack.lambda*J[p][1]*J[r][2]+stack.mu*J[p][2]*J[r][1]);
real32 const Rzy_ij = val*(stack.mu*J[p][1]*J[r][2]+stack.lambda*J[p][2]*J[r][1]);
- real32 const localIncrementx = (Rxx_ij * m_ux_n[m_elemsToNodes[k][j]] + Rxy_ij*m_uy_n[m_elemsToNodes[k][j]] + Rxz_ij*m_uz_n[m_elemsToNodes[k][j]]);
- real32 const localIncrementy = (Ryx_ij * m_ux_n[m_elemsToNodes[k][j]] + Ryy_ij*m_uy_n[m_elemsToNodes[k][j]] + Ryz_ij*m_uz_n[m_elemsToNodes[k][j]]);
- real32 const localIncrementz = (Rzx_ij * m_ux_n[m_elemsToNodes[k][j]] + Rzy_ij*m_uy_n[m_elemsToNodes[k][j]] + Rzz_ij*m_uz_n[m_elemsToNodes[k][j]]);
+ real32 const localIncrementx = (Rxx_ij * m_ux_n[m_elemsToNodes( k, j )] + Rxy_ij*m_uy_n[m_elemsToNodes( k, j )] + Rxz_ij*m_uz_n[m_elemsToNodes( k, j )]);
+ real32 const localIncrementy = (Ryx_ij * m_ux_n[m_elemsToNodes( k, j )] + Ryy_ij*m_uy_n[m_elemsToNodes( k, j )] + Ryz_ij*m_uz_n[m_elemsToNodes( k, j )]);
+ real32 const localIncrementz = (Rzx_ij * m_ux_n[m_elemsToNodes( k, j )] + Rzy_ij*m_uy_n[m_elemsToNodes( k, j )] + Rzz_ij*m_uz_n[m_elemsToNodes( k, j )]);
- RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorx[m_elemsToNodes[k][i]], localIncrementx );
- RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectory[m_elemsToNodes[k][i]], localIncrementy );
- RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorz[m_elemsToNodes[k][i]], localIncrementz );
+ RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorx[m_elemsToNodes( k, i )], localIncrementx );
+ RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectory[m_elemsToNodes( k, i )], localIncrementy );
+ RAJA::atomicAdd< parallelDeviceAtomic >( &m_stiffnessVectorz[m_elemsToNodes( k, i )], localIncrementz );
} );
}
@@ -542,13 +528,13 @@ class ExplicitElasticSEM : public finiteElement::KernelBase< SUBREGION_TYPE,
/// The array containing the nodal displacement array in z direction.
arrayView1d< real32 > const m_uz_n;
- /// The array containing the product of the stiffness matrix and the nodal pressure.
+ /// The array containing the product of the stiffness matrix and the nodal displacement.
arrayView1d< real32 > const m_stiffnessVectorx;
- /// The array containing the product of the stiffness matrix and the nodal pressure.
+ /// The array containing the product of the stiffness matrix and the nodal displacement.
arrayView1d< real32 > const m_stiffnessVectory;
- /// The array containing the product of the stiffness matrix and the nodal pressure.
+ /// The array containing the product of the stiffness matrix and the nodal displacement.
arrayView1d< real32 > const m_stiffnessVectorz;
/// The array containing the density of the medium
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp
index 64b69f3c103..52876ebf6b9 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.cpp
@@ -51,15 +51,21 @@ WaveSolverBase::WaveSolverBase( const std::string & name,
setSizedFromParent( 0 ).
setDescription( "Source Value of the sources" );
- registerWrapper( viewKeyStruct::timeSourceFrequencyString(), &m_timeSourceFrequency ).
- setInputFlag( InputFlags::REQUIRED ).
- setDescription( "Central frequency for the time source" );
-
registerWrapper( viewKeyStruct::receiverCoordinatesString(), &m_receiverCoordinates ).
setInputFlag( InputFlags::REQUIRED ).
setSizedFromParent( 0 ).
setDescription( "Coordinates (x,y,z) of the receivers" );
+ registerWrapper( viewKeyStruct::timeSourceFrequencyString(), &m_timeSourceFrequency ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setApplyDefaultValue( 0 ).
+ setDescription( "Central frequency for the time source" );
+
+ registerWrapper( viewKeyStruct::timeSourceDelayString(), &m_timeSourceDelay ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setApplyDefaultValue( -1 ).
+ setDescription( "Source time delay (1 / f0 by default)" );
+
registerWrapper( viewKeyStruct::rickerOrderString(), &m_rickerOrder ).
setInputFlag( InputFlags::OPTIONAL ).
setApplyDefaultValue( 2 ).
@@ -115,7 +121,6 @@ WaveSolverBase::WaveSolverBase( const std::string & name,
setApplyDefaultValue( -80 ).
setDescription( "Set the capacity of the lifo host storage (if negative, opposite of percentage of remaining memory)" );
-
registerWrapper( viewKeyStruct::usePMLString(), &m_usePML ).
setInputFlag( InputFlags::FALSE ).
setApplyDefaultValue( 0 ).
@@ -151,7 +156,7 @@ WaveSolverBase::WaveSolverBase( const std::string & name,
setSizedFromParent( 0 ).
setDescription( "Indices of the nodes (in the right order) for each receiver point" );
- registerWrapper( viewKeyStruct::sourceConstantsString(), &m_sourceConstants ).
+ registerWrapper( viewKeyStruct::receiverConstantsString(), &m_receiverConstants ).
setInputFlag( InputFlags::FALSE ).
setSizedFromParent( 0 ).
setDescription( "Constant part of the receiver for the nodes listed in m_receiverNodeIds" );
@@ -161,7 +166,15 @@ WaveSolverBase::WaveSolverBase( const std::string & name,
setSizedFromParent( 0 ).
setDescription( "Flag that indicates whether the receiver is local to this MPI rank" );
+ registerWrapper( viewKeyStruct::receiverRegionString(), &m_receiverRegion ).
+ setInputFlag( InputFlags::FALSE ).
+ setSizedFromParent( 0 ).
+ setDescription( "Region containing the receivers" );
+ registerWrapper( viewKeyStruct::receiverElemString(), &m_rcvElem ).
+ setInputFlag( InputFlags::FALSE ).
+ setSizedFromParent( 0 ).
+ setDescription( "Element containing the receivers" );
}
WaveSolverBase::~WaveSolverBase()
@@ -184,7 +197,7 @@ void WaveSolverBase::registerDataOnMesh( Group & meshBodies )
{
NodeManager & nodeManager = mesh.getNodeManager();
- nodeManager.registerField< fields::referencePosition32 >( this->getName() );
+ nodeManager.registerField< fields::referencePosition32 >( getName() );
arrayView2d< real64 const, nodes::REFERENCE_POSITION_USD > const X = nodeManager.referencePosition().toViewConst();
nodeManager.getField< fields::referencePosition32 >().resizeDimension< 1 >( X.size( 1 ) );
@@ -236,9 +249,7 @@ void WaveSolverBase::postProcessInput()
m_usePML = counter;
if( m_linearDASGeometry.size( 1 ) > 0 )
- {
m_useDAS = 1;
- }
if( m_useDAS )
{
@@ -255,22 +266,22 @@ void WaveSolverBase::postProcessInput()
}
- GEOS_THROW_IF( m_sourceCoordinates.size( 1 ) != 3,
+ GEOS_THROW_IF( m_sourceCoordinates.size( 0 ) > 0 && m_sourceCoordinates.size( 1 ) != 3,
"Invalid number of physical coordinates for the sources",
InputError );
- GEOS_THROW_IF( m_receiverCoordinates.size( 1 ) != 3,
+ GEOS_THROW_IF( m_receiverCoordinates.size( 0 ) > 0 && m_receiverCoordinates.size( 1 ) != 3,
"Invalid number of physical coordinates for the receivers",
InputError );
- EventManager const & event = this->getGroupByPath< EventManager >( "/Problem/Events" );
+ EventManager const & event = getGroupByPath< EventManager >( "/Problem/Events" );
real64 const & maxTime = event.getReference< real64 >( EventManager::viewKeyStruct::maxTimeString() );
real64 const & minTime = event.getReference< real64 >( EventManager::viewKeyStruct::minTimeString() );
real64 dt = 0;
for( localIndex numSubEvent = 0; numSubEvent < event.numSubGroups(); ++numSubEvent )
{
EventBase const * subEvent = static_cast< EventBase const * >( event.getSubGroups()[numSubEvent] );
- if( subEvent->getEventName() == "/Solvers/" + this->getName() )
+ if( subEvent->getEventName() == "/Solvers/" + getName() )
{
dt = subEvent->getReference< real64 >( EventBase::viewKeyStruct::forceDtString() );
}
@@ -286,7 +297,7 @@ void WaveSolverBase::postProcessInput()
{
m_nsamplesSeismoTrace = 0;
}
- localIndex const nsamples = int( (maxTime-minTime) /dt) + 1;
+ localIndex const nsamples = int( (maxTime-minTime) / dt) + 1;
localIndex const numSourcesGlobal = m_sourceCoordinates.size( 0 );
m_sourceValue.resize( nsamples, numSourcesGlobal );
@@ -348,7 +359,7 @@ real64 WaveSolverBase::explicitStep( real64 const & time_n,
localIndex WaveSolverBase::getNumNodesPerElem()
{
- DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+ DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" );
NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
@@ -387,7 +398,86 @@ localIndex WaveSolverBase::getNumNodesPerElem()
} );
return numNodesPerElem;
+}
+
+void WaveSolverBase::incrementIndexSeismoTrace( real64 const time_n )
+{
+ while( (m_dtSeismoTrace * m_indexSeismoTrace) <= (time_n + epsilonLoc) && m_indexSeismoTrace < m_nsamplesSeismoTrace )
+ {
+ m_indexSeismoTrace++;
+ }
+}
+
+void WaveSolverBase::computeAllSeismoTraces( real64 const time_n,
+ real64 const dt,
+ arrayView1d< real32 const > const var_np1,
+ arrayView1d< real32 const > const var_n,
+ arrayView2d< real32 > varAtReceivers )
+{
+ /*
+ * In forward case we compute seismo if time_n + dt is the first time
+ * step after the timeSeismo to write.
+ *
+ * time_n timeSeismo time_n + dt
+ * ---|--------------|-------------|
+ *
+ * In backward (time_n goes decreasing) case we compute seismo if
+ * time_n is the last time step before the timeSeismo to write.
+ *
+ * time_n - dt timeSeismo time_n
+ * ---|--------------|-------------|
+ */
+
+ if( m_nsamplesSeismoTrace == 0 )
+ return;
+ integer const dir = m_forward ? +1 : -1;
+ for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ )
+ {
+ real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo);
+ if( dir * timeSeismo > dir * (time_n + epsilonLoc))
+ break;
+ WaveSolverUtils::computeSeismoTrace( time_n, dir * dt, timeSeismo, iSeismo, m_receiverNodeIds,
+ m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers );
+ }
+}
+
+void WaveSolverBase::compute2dVariableAllSeismoTraces( localIndex const regionIndex,
+ real64 const time_n,
+ real64 const dt,
+ arrayView2d< real32 const > const var_np1,
+ arrayView2d< real32 const > const var_n,
+ arrayView2d< real32 > varAtReceivers )
+{
+ if( m_nsamplesSeismoTrace == 0 )
+ return;
+ integer const dir = m_forward ? +1 : -1;
+ for( localIndex iSeismo = m_indexSeismoTrace; iSeismo < m_nsamplesSeismoTrace; iSeismo++ )
+ {
+ real64 const timeSeismo = m_dtSeismoTrace * (m_forward ? iSeismo : (m_nsamplesSeismoTrace - 1) - iSeismo);
+ if( dir * timeSeismo > dir * (time_n + epsilonLoc))
+ break;
+ WaveSolverUtils::compute2dVariableSeismoTrace( time_n, dir * dt, regionIndex, m_receiverRegion, timeSeismo, iSeismo, m_rcvElem,
+ m_receiverConstants, m_receiverIsLocal, var_np1, var_n, varAtReceivers );
+ }
+}
+
+
+void WaveSolverBase::computeTargetNodeSet( arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
+ localIndex const subRegionSize,
+ localIndex const numQuadraturePointsPerElem )
+{
+ array1d< localIndex > scratch( subRegionSize * numQuadraturePointsPerElem );
+ localIndex i = 0;
+ for( localIndex e = 0; e < subRegionSize; ++e )
+ {
+ for( localIndex q = 0; q < numQuadraturePointsPerElem; ++q )
+ {
+ scratch[i++] = elemsToNodes( e, q );
+ }
+ }
+ std::ptrdiff_t const numUniqueValues = LvArray::sortedArrayManipulation::makeSortedUnique( scratch.begin(), scratch.end() );
+ m_solverTargetNodesSet.insert( scratch.begin(), scratch.begin() + numUniqueValues );
}
} /* namespace geos */
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp
index 14beb4b005f..4f58f18f304 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverBase.hpp
@@ -27,7 +27,10 @@
#if !defined( GEOS_USE_HIP )
#include "finiteElement/elementFormulations/Qk_Hexahedron_Lagrange_GaussLobatto.hpp"
#endif
-
+#include "WaveSolverUtils.hpp"
+#include
+#include
+#include
#if !defined( GEOS_USE_HIP )
#define SEM_FE_TYPES \
@@ -49,8 +52,9 @@ class WaveSolverBase : public SolverBase
{
public:
- using EXEC_POLICY = parallelDevicePolicy< >;
- using wsCoordType = real32;
+ static constexpr real64 epsilonLoc = WaveSolverUtils::epsilonLoc;
+ using EXEC_POLICY = WaveSolverUtils::EXEC_POLICY;
+ using wsCoordType = WaveSolverUtils::wsCoordType;
WaveSolverBase( const std::string & name,
Group * const parent );
@@ -83,6 +87,7 @@ class WaveSolverBase : public SolverBase
static constexpr char const * sourceValueString() { return "sourceValue"; }
static constexpr char const * timeSourceFrequencyString() { return "timeSourceFrequency"; }
+ static constexpr char const * timeSourceDelayString() { return "timeSourceDelay"; }
static constexpr char const * receiverCoordinatesString() { return "receiverCoordinates"; }
@@ -112,18 +117,21 @@ class WaveSolverBase : public SolverBase
static constexpr char const * usePMLString() { return "usePML"; }
static constexpr char const * parametersPMLString() { return "parametersPML"; }
+ static constexpr char const * receiverElemString() { return "rcvElem"; }
+ static constexpr char const * receiverRegionString() { return "receiverRegion"; }
};
- /**
- * @brief Safeguard for timeStep. Used to avoid memory issue due to too small value.
- */
- static constexpr real64 epsilonLoc = 1e-8;
-
/**
* @brief Re-initialize source and receivers positions in the mesh, and resize the pressureNp1_at_receivers array
*/
void reinit() override final;
+ SortedArray< localIndex > const & getSolverNodesSet() { return m_solverTargetNodesSet; }
+
+ void computeTargetNodeSet( arrayView2d< localIndex const, cells::NODE_MAP_USD > const & elemsToNodes,
+ localIndex const subRegionSize,
+ localIndex const numQuadraturePointsPerElem );
+
protected:
virtual void postProcessInput() override;
@@ -145,6 +153,36 @@ class WaveSolverBase : public SolverBase
*/
virtual void initializePML() = 0;
+ virtual void incrementIndexSeismoTrace( real64 const time_n );
+
+ /**
+ * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
+ * @param time_n the time corresponding to the field values pressure_n
+ * @param dt the simulation timestep
+ * @param var_np1 the field values at time_n + dt
+ * @param var_n the field values at time_n
+ * @param varAtreceivers the array holding the trace values, where the output is written
+ */
+ virtual void computeAllSeismoTraces( real64 const time_n,
+ real64 const dt,
+ arrayView1d< real32 const > const var_np1,
+ arrayView1d< real32 const > const var_n,
+ arrayView2d< real32 > varAtReceivers );
+ /**
+ * @brief Computes the traces on all receivers (see @computeSeismoTraces) up to time_n+dt
+ * @param time_n the time corresponding to the field values pressure_n
+ * @param dt the simulation timestep
+ * @param var_np1 the field values at time_n + dt
+ * @param var_n the field values at time_n
+ * @param varAtreceivers the array holding the trace values, where the output is written
+ */
+ virtual void compute2dVariableAllSeismoTraces( localIndex const regionIndex,
+ real64 const time_n,
+ real64 const dt,
+ arrayView2d< real32 const > const var_np1,
+ arrayView2d< real32 const > const var_n,
+ arrayView2d< real32 > varAtReceivers );
+
/**
* @brief Apply Perfectly Matched Layer (PML) to the regions defined in the geometry box from the xml
* @param time the time to apply the BC
@@ -152,8 +190,6 @@ class WaveSolverBase : public SolverBase
*/
virtual void applyPML( real64 const time, DomainPartition & domain ) = 0;
-
-
/**
* @brief Locate sources and receivers positions in the mesh elements, evaluate the basis functions at each point and save them to the
* corresponding elements nodes.
@@ -193,7 +229,6 @@ class WaveSolverBase : public SolverBase
virtual void registerDataOnMesh( Group & meshBodies ) override;
-
localIndex getNumNodesPerElem();
/// Coordinates of the sources in the mesh
@@ -205,6 +240,9 @@ class WaveSolverBase : public SolverBase
/// Central frequency for the Ricker time source
real32 m_timeSourceFrequency;
+ /// Source time delay (1 / f0 by default)
+ real32 m_timeSourceDelay;
+
/// Coordinates of the receivers in the mesh
array2d< real64 > m_receiverCoordinates;
@@ -212,7 +250,7 @@ class WaveSolverBase : public SolverBase
localIndex m_rickerOrder;
/// Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise
- localIndex m_outputSeismoTrace;
+ integer m_outputSeismoTrace;
/// Time step for seismoTrace output
real64 m_dtSeismoTrace;
@@ -259,6 +297,12 @@ class WaveSolverBase : public SolverBase
/// Flag that indicates whether the receiver is local or not to the MPI rank
array1d< localIndex > m_receiverIsLocal;
+ /// Array containing the elements which contain a receiver
+ array1d< localIndex > m_rcvElem;
+
+ /// Array containing the elements which contain the region which the receiver belongs
+ array1d< localIndex > m_receiverRegion;
+
/// Flag to enable LIFO
localIndex m_enableLifo;
@@ -294,11 +338,13 @@ class WaveSolverBase : public SolverBase
R1Tensor32 waveSpeedMaxXYZPML;
};
+ /// A set of target nodes IDs that will be handled by the current solver
+ SortedArray< localIndex > m_solverTargetNodesSet;
};
namespace fields
{
-using reference32Type = array2d< WaveSolverBase::wsCoordType, nodes::REFERENCE_POSITION_PERM >;
+using reference32Type = array2d< WaveSolverUtils::wsCoordType, nodes::REFERENCE_POSITION_PERM >;
DECLARE_FIELD( referencePosition32,
"referencePosition32",
reference32Type,
diff --git a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp
index fddb4470bf8..7c52e702482 100644
--- a/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp
+++ b/src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp
@@ -20,112 +20,171 @@
#ifndef GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_
#define GEOS_PHYSICSSOLVERS_WAVEPROPAGATION_WAVESOLVERUTILS_HPP_
-#include "WaveSolverBase.hpp"
+#include "mesh/utilities/ComputationalGeometry.hpp"
+#include "fileIO/Outputs/OutputBase.hpp"
+#include "LvArray/src/tensorOps.hpp"
namespace geos
{
struct WaveSolverUtils
{
+ static constexpr real64 epsilonLoc = 1e-8;
+ static constexpr real64 eps64 = std::numeric_limits< real64 >::epsilon();
+ static constexpr real32 eps32 = std::numeric_limits< real32 >::epsilon();
+ using EXEC_POLICY = parallelDevicePolicy< >;
+ using wsCoordType = real32;
GEOS_HOST_DEVICE
- static real32 evaluateRicker( real64 const & time_n, real32 const & f0, localIndex order )
+ static real32 evaluateRicker( real64 const time_n, real32 const f0, real32 const t0, localIndex order )
{
- real32 const o_tpeak = 1.0/f0;
+ real32 const delay = t0 > 0 ? t0 : 1 / f0;
real32 pulse = 0.0;
- if((time_n <= -0.9*o_tpeak) || (time_n >= 2.9*o_tpeak))
+ if( time_n <= -0.9 * delay || time_n >= 2.9 * delay )
{
return pulse;
}
- constexpr real32 pi = M_PI;
- real32 const lam = (f0*pi)*(f0*pi);
+ real32 const alpha = -pow( f0 * M_PI, 2 );
+ real32 const time_d = time_n - delay;
+ real32 const gaussian = exp( alpha * pow( time_d, 2 ));
+ localIndex const sgn = pow( -1, order + 1 );
switch( order )
{
- case 2:
- {
- pulse = 2.0*lam*(2.0*lam*(time_n-o_tpeak)*(time_n-o_tpeak)-1.0)*exp( -lam*(time_n-o_tpeak)*(time_n-o_tpeak));
- }
- break;
- case 1:
- {
- pulse = -2.0*lam*(time_n-o_tpeak)*exp( -lam*(time_n-o_tpeak)*(time_n-o_tpeak));
- }
- break;
case 0:
- {
- pulse = -(time_n-o_tpeak)*exp( -2*lam*(time_n-o_tpeak)*(time_n-o_tpeak) );
- }
- break;
+ pulse = sgn * gaussian;
+ break;
+ case 1:
+ pulse = sgn * (2 * alpha * time_d) * gaussian;
+ break;
+ case 2:
+ pulse = sgn * (2 * alpha + 4 * pow( alpha, 2 ) * pow( time_d, 2 )) * gaussian;
+ break;
+ case 3:
+ pulse = sgn * (12 * pow( alpha, 2 ) * time_d + 8 * pow( alpha, 3 ) * pow( time_d, 3 )) * gaussian;
+ break;
+ case 4:
+ pulse = sgn * (12 * pow( alpha, 2 ) + 48 * pow( alpha, 3 ) * pow( time_d, 2 ) + 16 * pow( alpha, 4 ) * pow( time_d, 4 )) * gaussian;
+ break;
default:
- GEOS_ERROR( "This option is not supported yet, rickerOrder must be 0, 1 or 2" );
+ GEOS_ERROR( "This option is not supported yet, rickerOrder must be in range {0:4}" );
}
return pulse;
}
+ /**
+ * @brief Initialize the trace file (touch)
+ */
+ static void initTrace( char const * prefix,
+ string const & name,
+ localIndex const nReceivers,
+ arrayView1d< localIndex const > const receiverIsLocal )
+ {
+ string const outputDir = OutputBase::getOutputDirectory();
+ RAJA::ReduceSum< ReducePolicy< serialPolicy >, localIndex > count( 0 );
+
+ forAll< serialPolicy >( nReceivers, [=] ( localIndex const ircv )
+ {
+ if( receiverIsLocal[ircv] == 1 )
+ {
+ count += 1;
+ string const fn = joinPath( outputDir, GEOS_FMT( "{}_{}_{:03}.txt", prefix, name, ircv ) );
+ // std::cout << "touch " << fn << std::endl;
+ std::ofstream f( fn, std::ios::out | std::ios::trunc );
+ }
+ } );
+
+ localIndex const total = MpiWrapper::sum( count.get() );
+ GEOS_ERROR_IF( nReceivers != total, GEOS_FMT( "Invalid distribution of receivers: nReceivers={} != MPI::sum={}.", nReceivers, total ) );
+ }
+
+ static void writeSeismoTraceVector( char const * prefix,
+ string const & name,
+ bool const outputSeismoTrace,
+ localIndex const nReceivers,
+ arrayView1d< localIndex const > const receiverIsLocal,
+ localIndex const nsamplesSeismoTrace,
+ arrayView2d< real32 const > const varAtReceiversx,
+ arrayView2d< real32 const > const varAtReceiversy,
+ arrayView2d< real32 const > const varAtReceiversz )
+ {
+ writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversx );
+ writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversy );
+ writeSeismoTrace( prefix, name, outputSeismoTrace, nReceivers, receiverIsLocal, nsamplesSeismoTrace, varAtReceiversz );
+ }
+
+ static void writeSeismoTrace( char const * prefix,
+ string const & name,
+ bool const outputSeismoTrace,
+ localIndex const nReceivers,
+ arrayView1d< localIndex const > const receiverIsLocal,
+ localIndex const nsamplesSeismoTrace,
+ arrayView2d< real32 const > const varAtReceivers )
+ {
+ if( !outputSeismoTrace ) return;
+
+ string const outputDir = OutputBase::getOutputDirectory();
+ forAll< serialPolicy >( nReceivers, [=] ( localIndex const ircv )
+ {
+ if( receiverIsLocal[ircv] == 1 )
+ {
+ string const fn = joinPath( outputDir, GEOS_FMT( "{}_{}_{:03}.txt", prefix, name, ircv ) );
+ // std::cout << "writing " << fn << std::endl;
+ std::ofstream f( fn, std::ios::app );
+ if( f )
+ {
+ for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample )
+ {
+ // index - time - value
+ f << iSample << " " << varAtReceivers[iSample][nReceivers] << " " << varAtReceivers[iSample][ircv] << std::endl;
+ }
+ f.close();
+ }
+ else
+ {
+ GEOS_WARNING( GEOS_FMT( "Failed to open output file {}", fn ) );
+ }
+ }
+ } );
+ }
+
static void computeSeismoTrace( real64 const time_n,
real64 const dt,
real64 const timeSeismo,
- localIndex iSeismo,
+ localIndex const iSeismo,
arrayView2d< localIndex const > const receiverNodeIds,
arrayView2d< real64 const > const receiverConstants,
arrayView1d< localIndex const > const receiverIsLocal,
- localIndex const nsamplesSeismoTrace,
- localIndex const outputSeismoTrace,
arrayView1d< real32 const > const var_np1,
arrayView1d< real32 const > const var_n,
arrayView2d< real32 > varAtReceivers )
{
real64 const time_np1 = time_n + dt;
- real32 const a1 = (LvArray::math::abs( dt ) < WaveSolverBase::epsilonLoc ) ? 1.0 : (time_np1 - timeSeismo)/dt;
+ real32 const a1 = abs( dt ) < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt;
real32 const a2 = 1.0 - a1;
- if( nsamplesSeismoTrace > 0 )
- {
- forAll< WaveSolverBase::EXEC_POLICY >( receiverConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const ircv )
- {
- if( receiverIsLocal[ircv] == 1 )
- {
- varAtReceivers[iSeismo][ircv] = 0.0;
- real32 vtmp_np1 = 0.0;
- real32 vtmp_n = 0.0;
- for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode )
- {
- vtmp_np1 += var_np1[receiverNodeIds[ircv][inode]] * receiverConstants[ircv][inode];
- vtmp_n += var_n[receiverNodeIds[ircv][inode]] * receiverConstants[ircv][inode];
- }
- // linear interpolation between the pressure value at time_n and time_(n+1)
- varAtReceivers[iSeismo][ircv] = a1*vtmp_n + a2*vtmp_np1;
- }
- } );
- }
+ localIndex const nReceivers = receiverConstants.size( 0 );
- // TODO DEBUG: the following output is only temporary until our wave propagation kernels are finalized.
- // Output will then only be done via the previous code.
- if( iSeismo == nsamplesSeismoTrace - 1 )
+ forAll< EXEC_POLICY >( nReceivers, [=] GEOS_HOST_DEVICE ( localIndex const ircv )
{
- forAll< serialPolicy >( receiverConstants.size( 0 ), [=] ( localIndex const ircv )
+ if( receiverIsLocal[ircv] == 1 )
{
- if( outputSeismoTrace == 1 )
+ real32 vtmp_np1 = 0.0, vtmp_n = 0.0;
+ for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode )
{
- if( receiverIsLocal[ircv] == 1 )
- {
- // Note: this "manual" output to file is temporary
- // It should be removed as soon as we can use TimeHistory to output data not registered on the mesh
- // TODO: remove saveSeismo and replace with TimeHistory
- std::ofstream f( GEOS_FMT( "seismoTraceReceiver{:03}.txt", ircv ), std::ios::app );
- for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample )
- {
- f << iSample << " " << varAtReceivers[iSample][ircv] << std::endl;
- }
- f.close();
- }
+ vtmp_np1 += var_np1[receiverNodeIds( ircv, inode )] * receiverConstants( ircv, inode );
+ vtmp_n += var_n[receiverNodeIds( ircv, inode )] * receiverConstants( ircv, inode );
}
- } );
- }
+ // linear interpolation between the pressure value at time_n and time_{n+1}
+ varAtReceivers( iSeismo, ircv ) = a1 * vtmp_n + a2 * vtmp_np1;
+ // NOTE: varAtReceivers has size(1) = numReceiversGlobal + 1, this does not OOB
+ // left in the forAll loop for sync issues
+ varAtReceivers( iSeismo, nReceivers ) = a1 * time_n + a2 * time_np1;
+ }
+ } );
}
static void compute2dVariableSeismoTrace( real64 const time_n,
@@ -133,70 +192,41 @@ struct WaveSolverUtils
localIndex const regionIndex,
arrayView1d< localIndex const > const receiverRegion,
real64 const timeSeismo,
- localIndex iSeismo,
+ localIndex const iSeismo,
arrayView1d< localIndex const > const rcvElem,
arrayView2d< real64 const > const receiverConstants,
arrayView1d< localIndex const > const receiverIsLocal,
- localIndex const nsamplesSeismoTrace,
- localIndex const outputSeismoTrace,
arrayView2d< real32 const > const var_np1,
arrayView2d< real32 const > const var_n,
arrayView2d< real32 > varAtReceivers )
{
real64 const time_np1 = time_n+dt;
- real32 const a1 = (dt < WaveSolverBase::epsilonLoc) ? 1.0 : (time_np1 - timeSeismo)/dt;
+ real32 const a1 = dt < epsilonLoc ? 1.0 : (time_np1 - timeSeismo) / dt;
real32 const a2 = 1.0 - a1;
- if( nsamplesSeismoTrace > 0 )
- {
- forAll< WaveSolverBase::EXEC_POLICY >( receiverConstants.size( 0 ), [=] GEOS_HOST_DEVICE ( localIndex const ircv )
- {
- if( receiverIsLocal[ircv] == 1 )
- {
- if( receiverRegion[ircv] == regionIndex )
- {
- varAtReceivers[iSeismo][ircv] = 0.0;
- real32 vtmp_np1 = 0.0;
- real32 vtmp_n = 0.0;
- for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode )
- {
- vtmp_np1 += var_np1[rcvElem[ircv]][inode] * receiverConstants[ircv][inode];
- vtmp_n += var_n[rcvElem[ircv]][inode] * receiverConstants[ircv][inode];
- }
- // linear interpolation between the pressure value at time_n and time_(n+1)
- varAtReceivers[iSeismo][ircv] = a1*vtmp_n + a2*vtmp_np1;
- }
- }
- } );
- }
+ localIndex const nReceivers = receiverConstants.size( 0 );
- // TODO DEBUG: the following output is only temporary until our wave propagation kernels are finalized.
- // Output will then only be done via the previous code.
- if( iSeismo == nsamplesSeismoTrace - 1 )
+ forAll< EXEC_POLICY >( nReceivers, [=] GEOS_HOST_DEVICE ( localIndex const ircv )
{
- if( outputSeismoTrace == 1 )
+ if( receiverIsLocal[ircv] == 1 )
{
- forAll< serialPolicy >( receiverConstants.size( 0 ), [=] ( localIndex const ircv )
+ if( receiverRegion[ircv] == regionIndex )
{
- if( receiverIsLocal[ircv] == 1 )
+ real32 vtmp_np1 = 0.0, vtmp_n = 0.0;
+ for( localIndex inode = 0; inode < receiverConstants.size( 1 ); ++inode )
{
- // Note: this "manual" output to file is temporary
- // It should be removed as soon as we can use TimeHistory to output data not registered on the mesh
- // TODO: remove saveSeismo and replace with TimeHistory
- if( receiverRegion[ircv] == regionIndex )
- {
- std::ofstream f( GEOS_FMT( "seismoTraceReceiver{:03}.txt", ircv ), std::ios::app );
- for( localIndex iSample = 0; iSample < nsamplesSeismoTrace; ++iSample )
- {
- f << iSample << " " << varAtReceivers[iSample][ircv] << std::endl;
- }
- f.close();
- }
+ vtmp_np1 += var_np1( rcvElem[ircv], inode ) * receiverConstants( ircv, inode );
+ vtmp_n += var_n( rcvElem[ircv], inode ) * receiverConstants( ircv, inode );
}
- } );
+ // linear interpolation between the pressure value at time_n and time_{n+1}
+ varAtReceivers( iSeismo, ircv ) = a1 * vtmp_n + a2 * vtmp_np1;
+ // NOTE: varAtReceivers has size(1) = numReceiversGlobal + 1, this does not OOB
+ // left in the forAll loop for sync issues
+ varAtReceivers( iSeismo, nReceivers ) = a1 * time_n + a2 * time_np1;
+ }
}
- }
+ } );
}
/**
@@ -250,9 +280,7 @@ struct WaveSolverUtils
// all dot products should be non-negative (we enforce outward normals)
if( s < 0 )
- {
return false;
- }
}
return true;
@@ -271,7 +299,7 @@ struct WaveSolverUtils
static void
computeCoordinatesOnReferenceElement( real64 const (&coords)[3],
arraySlice1d< localIndex const, cells::NODE_MAP_USD - 1 > const elemsToNodes,
- arrayView2d< WaveSolverBase::wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
+ arrayView2d< wsCoordType const, nodes::REFERENCE_POSITION_USD > const X,
real64 (& coordsOnRefElem)[3] )
{
real64 xLocal[FE_TYPE::numNodes][3]{};
@@ -293,8 +321,6 @@ struct WaveSolverUtils
}
}
-
-
};
} /* namespace geos */
diff --git a/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst b/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst
index 000b0ace852..2ad303a4306 100644
--- a/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst
+++ b/src/coreComponents/schema/docs/AcousticFirstOrderSEM.rst
@@ -17,7 +17,7 @@ logLevel integer 0 Log level
name string required A name is required for any non-unique nodes
outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise
receiverCoordinates real64_array2d required Coordinates (x,y,z) of the receivers
-rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default
+rickerOrder integer 2 Flag that indicates the order of the Ricker to be used: {0:4}. Order 2 by default
saveFields integer 0 Set to 1 to save fields during forward and restore them during backward
shotIndex integer 0 Set the current shot for temporary files
sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources
diff --git a/src/coreComponents/schema/docs/AcousticSEM.rst b/src/coreComponents/schema/docs/AcousticSEM.rst
index 000b0ace852..2ad303a4306 100644
--- a/src/coreComponents/schema/docs/AcousticSEM.rst
+++ b/src/coreComponents/schema/docs/AcousticSEM.rst
@@ -17,7 +17,7 @@ logLevel integer 0 Log level
name string required A name is required for any non-unique nodes
outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise
receiverCoordinates real64_array2d required Coordinates (x,y,z) of the receivers
-rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default
+rickerOrder integer 2 Flag that indicates the order of the Ricker to be used: {0:4}. Order 2 by default
saveFields integer 0 Set to 1 to save fields during forward and restore them during backward
shotIndex integer 0 Set the current shot for temporary files
sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources
diff --git a/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst b/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst
index 000b0ace852..2ad303a4306 100644
--- a/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst
+++ b/src/coreComponents/schema/docs/ElasticFirstOrderSEM.rst
@@ -17,7 +17,7 @@ logLevel integer 0 Log level
name string required A name is required for any non-unique nodes
outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise
receiverCoordinates real64_array2d required Coordinates (x,y,z) of the receivers
-rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default
+rickerOrder integer 2 Flag that indicates the order of the Ricker to be used: {0:4}. Order 2 by default
saveFields integer 0 Set to 1 to save fields during forward and restore them during backward
shotIndex integer 0 Set the current shot for temporary files
sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources
diff --git a/src/coreComponents/schema/docs/ElasticSEM.rst b/src/coreComponents/schema/docs/ElasticSEM.rst
index acc22024dcc..f6fa8df2106 100644
--- a/src/coreComponents/schema/docs/ElasticSEM.rst
+++ b/src/coreComponents/schema/docs/ElasticSEM.rst
@@ -17,7 +17,7 @@ logLevel integer 0 Log level
name string required A name is required for any non-unique nodes
outputSeismoTrace integer 0 Flag that indicates if we write the seismo trace in a file .txt, 0 no output, 1 otherwise
receiverCoordinates real64_array2d required Coordinates (x,y,z) of the receivers
-rickerOrder integer 2 Flag that indicates the order of the Ricker to be used o, 1 or 2. Order 2 by default
+rickerOrder integer 2 Flag that indicates the order of the Ricker to be used: {0:4}. Order 2 by default
saveFields integer 0 Set to 1 to save fields during forward and restore them during backward
shotIndex integer 0 Set the current shot for temporary files
sourceCoordinates real64_array2d required Coordinates (x,y,z) of the sources
diff --git a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp
index 1f2a6c6d57b..2f93fa26a2b 100644
--- a/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp
+++ b/src/coreComponents/unitTests/wavePropagationTests/testWavePropagation.cpp
@@ -127,6 +127,13 @@ char const * xmlInput =
fieldName="mediumVelocity"
scale="1500"
setNames="{ all }"/>
+