Wave solvers: Source on a mesh node/edge results in wrong signal (same for receivers)

[Bertbk]

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!

[sframba]

Related to #2888