GEOS
GEOS copied to clipboard
GEOS-DEV
Locate seismic point sources and point receivers in the mesh and calculate constant source and receiver terms
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:
- 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 ){}
- Recompile GEOS (with or without pyGEOS)
- use the xml input below (the mesh file path is available in the xml)
- run geos or pygeos
- 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>
