Locate seismic point sources and point receivers in the mesh and calculate constant source and receiver terms

[shenchengyi]

Point sources/receivers shared by several elements The current algorithm to find the sources and receivers will not stop looking for them once a point is located. For example, if a source or a receiver is on an internal edge shared by 4 elements, the point will be found on 4 neighbour elements and the constant terms will be calculated 4 times. With GPU this leads to incorrect source/receiver constant values. Note that these incorrect values are not reproducible: they vary from one run to another.

To Reproduce Steps to reproduce the behavior:

1. Add some printf in AcousticWaveEquationSEMKernel.hpp or ElasticWaveEquationSEMKernel.hpp like: printf("Elem %d source %d sourceIsAccessible %d\n", k, isrc, sourceIsAccessible[isrc]); around L98 and printf("Node=%d Elem=%d sourceCoords=%.3f %.3f %.3f sourceConstants=%.3e \n", a,k,coords[0],coords[1],coords[2],sourceConstants[isrc][a]); inside for( localIndex a = 0; a < numNodesPerElem; ++a ){}
2. Recompile GEOS (with or without pyGEOS)
3. use the xml input below (the mesh file path is available in the xml)
4. run geos or pygeos
5. in the log we will have the following lines (the sourceConstants may not be reproducible) Node=0 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=0 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=0 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=0 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=2.000e-01
   Node=1 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=1 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=1 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=2 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=2 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=2.000e-01
   Node=2 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=3 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=2.000e-01
   Node=3 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=3 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=4 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=4 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=5 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=5 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=4 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=8.000e-01
   Node=6 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=8.000e-01
   Node=5 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=6 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00
   Node=1 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=2.000e-01
   Node=6 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=7 Elem=8127 sourceCoords=6325.000 3162.500 408.000 sourceConstants=8.000e-01 Node=7 Elem=8128 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=2 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=7 Elem=8256 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=3 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=4 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=5 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=8.000e-01 Node=6 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00 Node=7 Elem=8255 sourceCoords=6325.000 3162.500 408.000 sourceConstants=0.000e+00

Expected behavior Preferably make the source/receiver search yield a unique and consistent result.

Platform (please complete the following information):

• Machine: pangea3
• GEOSX Version: any recent version, for example HEAD, sha1: 079ef8cec

xml input

<?xml version="1.0" ?>

<Problem>
  <Solvers>
    <!-- define the solver -->
    <!-- we'll have to discuss the cflFactor with Randy -->
    <!-- define the source coordinates: an explosion at 8m beneath the surface-->
    <!-- define the time source frequency -->
    <!-- define the receiver coordinates: they are geophones -->
       <AcousticSEM
      name="acousticSolver"
      cflFactor="0.8"
      discretization="FE1"
      targetRegions="{Region}"
      sourceCoordinates="{ { 6325.0000,3162.5000,408.00000 } }"
      timeSourceFrequency="10.0"
      receiverCoordinates="{ { 6325.0000,3162.5000,418.00000 } }"
      outputSeismoTrace="1"
      dtSeismoTrace="0.002" />
  </Solvers>

  <!-- hexahedral mesh generated with vtk from regular IJK grid -->
  <Mesh>
    <VTKMesh
        name="mesh_cut400"
        fieldsToImport="{mediumDensity,mediumVelocityVp}"
        fieldNamesInGEOSX="{mediumDensity,mediumVelocity}"
        partitionRefinement="0"
        useGlobalIds="-1"
    file="Y00_elastic_12x12x10_cut400.vtu"/>
  </Mesh>

  <Geometry>
    <!-- tag the free surface boundary using box -->
    <Box
      name="zpos"
      xMin="{ 5524.99, 2362.49, 399.99 }"
      xMax="{ 7125.01, 3962.51, 400.01 }"/>
  </Geometry>

  <!-- list of events-->
  <Events maxTime="3">
    <!-- trigger the application of the solver -->
    <!-- control the timestepping here with forceDt -->
    <!-- (there are other ways to control the time step size) -->
    <PeriodicEvent
      name="solverApplications"
      forceDt="0.001"
      targetExactStartStop="0"
      targetExactTimestep="0"
      target="/Solvers/acousticSolver"/>
    <PeriodicEvent
      name="vtk"
      timeFrequency="4.50"
      targetExactTimestep="0"
      target="/Outputs/vtkOutput"/>
  </Events>

  <!-- numerical methods for discretization-->
  <NumericalMethods>
    <FiniteElements>
      <FiniteElementSpace
        name="FE1"
        order="1"
        formulation="SEM"/>
    </FiniteElements>
  </NumericalMethods>

  <!-- Region -->
  <ElementRegions>
    <CellElementRegion
      name="Region"
      cellBlocks="{ hexahedra }"
      materialList="{ nullModel }"/>
  </ElementRegions>

  <Constitutive>
    <NullModel
      name="nullModel"/>
  </Constitutive>

  <!-- here use FieldSpecification to define: -->
  <FieldSpecifications>
    <FieldSpecification
      name="zposFreeSurface"
      objectPath="mesh_cut400/FE1/faceManager"
      fieldName="FreeSurface"
      scale="0.0"
    setNames="{zpos}"/>
  </FieldSpecifications>
  <Outputs>
    <VTK
      name="vtkOutput"
      fieldNames="{mediumDensity,mediumVelocity}"
      plotLevel="3"/>
  </Outputs>

</Problem>

[sframba]

Related to #3104

[sframba]

Solved by PR #3202