GEOS
GEOS copied to clipboard
GEOS-DEV
Wave solvers: Source on a mesh node/edge results in wrong signal (same for receivers)
Describe the bug When running a wave propagation solver (eg acoustic), if the source lies on the boundary of 2 or more elements, the amplitude of the source signal can be counted multiple times, or not a all, resulting in a wrong result. This has been reported on issue #2256 but closed for duplicate reason with #2227. However, the problem persists. Could it be due to parallel concurrency?
To Reproduce Launch numerous time an xml file with the following geometry and solver. The source is at the origin and the mesh is such that the origin is a mesh node. The complete file is provided at the end.
<Mesh>
<InternalMesh
name="mesh"
elementTypes="{ C3D8 }"
xCoords="{ -1500,1500}"
yCoords="{ -1500,1500 }"
zCoords="{ -1500,1500 }"
nx="{ 150 }"
ny="{ 150 }"
nz="{ 150 }"
cellBlockNames="{ cb }"/>
</Mesh>
<Solvers>
<AcousticSEM
name="acousticSolver"
cflFactor="0.25"
discretization="FE1"
targetRegions="{ Region }"
sourceCoordinates="{ { 0, 0, 0 } }"
outputSeismoTrace = "1"
dtSeismoTrace="0.0013"
timeSourceFrequency="5"
receiverCoordinates="{ { 1, 2, 3 }}"/>
</Solvers>
<NumericalMethods>
<FiniteElements>
<FiniteElementSpace
name="FE1"
order="1"
formulation="SEM"/>
</FiniteElements>
</NumericalMethods>
Screenshots
Here is a plot of different signal obtained through the same piece of code. The amplitude varies from 0 to the double. 99 simulations have been launched. The source is at the origin belonging to 8 elements, the receiver is close to at (1,2,3), inside an element.
Platform (please complete the following information):
- Machine: pangea3
- Compiler: tested with gcc
- GEOSX Version: develop
Additional information
If a source is on a mesh node, then it belongs to 8 elements. Currently, GEOS loops on each element and check if the source is inside the element (boundary included). If so, GEOS pre-computes the parameters and quantities for the source. As the loop on the element is achieved in parallel, concurrency can appear.
The loop on the element is in function Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage in the file
vti/src/coreComponents/physicsSolvers/wavePropagation/PrecomputeSourcesAndReceiversKermel.hpp
In case it helps, here is a piece of log of Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage obtained through (ugly) printf:
// file recomputeSourcesAndReceiversKernel.hpp
// Function Compute1DSourceAndReceiverConstants
for( localIndex a = 0; a < numNodesPerElem; ++a )
{
sourceNodeIds[isrc][a] = elemsToNodes[k][a];
sourceConstants[isrc][a] = Ntest[a];
// ugly but useful printf :-)
printf("sourceNodeIds[%d][%d] = %d \n", isrc, a, elemsToNodes[k][a]);
printf("sourceConstants[%d][%d] = %.2f \n", isrc, a, Ntest[a]);
}
Part of the resulting printf, for the case where the source went to zero:
sourceNodeIds[0][0] = 196374
sourceNodeIds[0][0] = 196375
sourceNodeIds[0][0] = 193773
sourceNodeIds[0][0] = 193774
sourceNodeIds[0][0] = 193722
sourceNodeIds[0][0] = 193723
sourceNodeIds[0][0] = 196323
sourceNodeIds[0][0] = 196324
sourceConstants[0][0] = 0.00
sourceConstants[0][0] = 1.00
sourceConstants[0][0] = 0.00
sourceConstants[0][0] = 0.00
sourceConstants[0][0] = 0.00
sourceConstants[0][0] = 0.00
sourceConstants[0][0] = 0.00
sourceConstants[0][0] = 0.00
...
The source information sourceConstants[0][0] and sourceNodeIds[0][0] have been changed 8 times... Same for 8 other nodes.
Same problem for Receivers
The same problem happens for receivers. In addition, if a receiver is at the boundary between 2 subdomains (of the domain decomposition method) then GEOS crashes unexpectedly. To observe that, simply change the receiver's position to, eg, the origin: receiverCoordinates="{ { 0, 0, 0 }}"
Launch then in parallel with a ddm 2 in the x direction, for example:
mpirun -n 8 geosx -i issue.xml -x 2 -y 2 -z 2
You should get something like
***** ERROR
***** LOCATION: src/coreComponents/physicsSolvers/wavePropagation/WaveSolverUtils.hpp:108
***** Controlling expression (should be false): nReceivers != total
***** Rank 0: : Invalid distribution of receivers: nReceivers=1 != MPI::sum=8.
Maybe due to the check on ghostElement in Compute1DSourceAndReceiverConstantsWithElementsAndRegionStorage ?
xml complete file
The complete xml is the following
<?xml version="1.0" ?>
<Problem>
<!-- hexahedral mesh generated internally by GEOSX -->
<Mesh>
<InternalMesh
name="mesh"
elementTypes="{ C3D8 }"
xCoords="{ -1500,1500}"
yCoords="{ -1500,1500 }"
zCoords="{ -1500,1500 }"
nx="{ 150 }"
ny="{ 150 }"
nz="{ 150 }"
cellBlockNames="{ cb }"/>
</Mesh>
<Geometry>
<Box
name="zpos"
xMin="{ -1500.1, -1500.1, 1499.9}"
xMax="{ 1500.1, 1500.1, 1500.1}"/>
<Box
name="zneg"
xMin="{ -1500.1, -1500.1, -1500.1}"
xMax="{ 1500.1, 1500.1, -1499.9}"/>
<Box
name="xpos"
xMin="{ 1499.9, -1500.1, -1500.1}"
xMax="{ 1500.1, 1500.1, 1500.1}"/>
<Box
name="xneg"
xMin="{ -1500.1, -1500.1, -1500.1}"
xMax="{ -1499.9, 1500.1, 1500.1}"/>
<Box
name="ypos"
xMin="{ -1500.1, 1499.9, -1500.1}"
xMax="{ 1500.1, 1500.1, 1500.1}"/>
<Box
name="yneg"
xMin="{ -1500.1, -1500.1, -1500.1}"
xMax="{ 1500.1, -1499.9, 1500.1}"/>
</Geometry>
<!-- The free surface condition the domain -->
<FieldSpecifications>
<FieldSpecification
name="zposFreeSurface"
objectPath="faceManager"
fieldName="acousticFreeSurface"
scale="0.0"
setNames="{ zpos }"/>
<FieldSpecification
name="znegBottomSurface"
objectPath="faceManager"
fieldName="acousticBottomSurface"
scale="0.0"
setNames="{ zneg }"/>
<FieldSpecification
name="LateralFreeSurface"
objectPath="faceManager"
fieldName="acousticLateralSurface"
scale="0.0"
setNames="{ xpos, xneg, ypos, yneg }"/>
</FieldSpecifications>
<Solvers>
<AcousticSEM
name="acousticSolver"
cflFactor="0.25"
discretization="FE1"
targetRegions="{ Region }"
sourceCoordinates="{ { 0, 0, 0 } }"
outputSeismoTrace = "1"
dtSeismoTrace="0.0013"
timeSourceFrequency="5"
receiverCoordinates="{ { 1, 2, 3 }}"/>
</Solvers>
<NumericalMethods>
<FiniteElements>
<FiniteElementSpace
name="FE1"
order="1"
formulation="SEM"/>
</FiniteElements>
</NumericalMethods>
<ElementRegions>
<CellElementRegion
name="Region"
cellBlocks="{ cb }"
materialList="{ nullModel }"/>
</ElementRegions>
<Constitutive>
<NullModel
name="nullModel"/>
</Constitutive>
<FieldSpecifications>
<!-- 1) The initial pressure field -->
<FieldSpecification
name="initialPressure_n"
initialCondition="1"
setNames="{ all }"
objectPath="nodeManager"
fieldName="pressure_n"
scale="0.0"/>
<FieldSpecification
name="initialPressure_nm1"
initialCondition="1"
setNames="{ all }"
objectPath="nodeManager"
fieldName="pressure_nm1"
scale="0.0"/>
<!-- 2) The velocity in the domain -->
<FieldSpecification
name="cellVelocity"
initialCondition="1"
objectPath="ElementRegions/Region/cb"
fieldName="acousticVelocity"
scale="3000"
setNames="{ all }"/>
<FieldSpecification
name="cellDensity"
initialCondition="1"
objectPath="ElementRegions/Region/cb"
fieldName="acousticDensity"
scale="1.0"
setNames="{ all }"/>
</FieldSpecifications>
<Events
maxTime="0.5">
<!-- trigger the application of the solver -->
<!-- control the timestepping here with forceDt -->
<PeriodicEvent
name="solverApplications"
forceDt="0.0013"
target="/Solvers/acousticSolver"/>
</Events>
<!-- collect the pressure values at the nodes -->
<Tasks>
<PackCollection
name="pressureCollection"
objectPath="/Problem/domain/MeshBodies/mesh/meshLevels/FE1/nodeManager"
fieldName="pressure_np1"/>
</Tasks>
</Problem>
Thank you and sorry for the long issue!
