Unknown AxiSEM3D running issue and Error locating source in mesh.

[Yishharu]

Hi Kuangdai,

I'm trying to simulate some marine seismology shots for the ocean boundary layer. We define two fluid layers with undulation in the cartesian model. AxiSEM3D reads the model well but seems to have some unknown map issues and errors locating the source in the mesh.

Could you help have a look?

Best wishes, Zhi

Returing error:

$ mpirun -np 4 ./axisem3d 

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      A            |i .|'''.|'||''''E'||    ||'  ____'||''|.    
     |||   ... ...... ||..  ' ||  .   |||  |||   ` // ||   ||   
    |  ||   '|..'  ||  ''|||. ||''|   |'|..'||    //  ||    ||  
   .''''|.   .x.   ||      '||||      | 'M' ||    \\  ||    ||  
  .|.  .||..|  ||..||.|'...|S.||....|.|. | .||.    3' D|...|'   
  .............................................   //            
                                                 /'     v 1.0   

  Copyright (c) 2019 Kuangdai Leng & friends, MIT License       
  Source, docs and issues: github.com/kuangdai/AxiSEM-3D        
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

********************* WARNING FROM AxiSEM3D ********************
FROM: io::verifyDirectories
WHAT: Output directory exists; old output renamed to
      ./output__backup@2021-07-23T23:43:05
TIME: 2021-07-23T23:43:05
****************************************************************

========================== Exodus Mesh =========================
Overview________________________________________________________
  Exodus file           =  local_mesh__ocean_topo__20Hz.e
  mesh CS               =  cartesian
  model name            =  ocean_topo
  isotropic             =  true
  attenuation           =  false
  storage type          =  element_nodes
  # discontinuities     =  3
Mesh geometry___________________________________________________
  # elements            =  420
  # nodes               =  464
  range of s-axis       =  [0, 1000]
  range of z-axis       =  [6.3705e+06, 6.371e+06]
Global variables________________________________________________
  dist_tolerance        =  0.0333333
  dt (nPol = 1)         =  0.0130719
  min_edge_length       =  33.3333
  minimum_period        =  0.0485909
  crdsys                =  cartesian
  model                 =  ocean_topo
Element geometry________________________________________________
  linear                :  420
  semi-spherical        :  0
  spherical             :  0
Element medium__________________________________________________
  fluid                 :  420
  solid                 :  0
Material properties_____________________________________________
  RHO                   ∈  [1029, 1050]
  VP                    ∈  [1470, 1530]
  VS                    ∈  [0, 0]
Side sets_______________________________________________________
  x0                    :  15 sides
  x1                    :  15 sides
  y0                    :  28 sides
  y1                    :  28 sides
Miscellaneous___________________________________________________
  ellipticity curve     :  2 knots
----------------------------------------------------------------
Mesh generation command line:
>> python -m salvus_mesh_lite.interface AxiSEMCartesian --basic.model
   ocean_topo.bm --basic.period 0.05 --cartesian2Daxisem.x 1.
   --cartesian2Daxisem.min_z 6370.5 --output_filename
   local_mesh__ocean_topo__20Hz.e
================================================================

============================ Geodesy ===========================
model in Cartesian     =  true
surface radius         =  6.371e+06
solid surface radius   =  6.37047e+06
bottom radius          =  6.3705e+06
a reference geographic location on the positive z-axis
  latitude             =  90
  longitude            =  0
  radius               =  6.37099e+06
* No ellipticity correction for Cartesian models.
================================================================

====================== Absorbing Boundary ======================
user-specified boundaries  =  [RIGHT, BOTTOM]
those contained in mesh    =  [RIGHT, BOTTOM]
Clayton-Enquist enabled    =  true
Kosloff-Kosloff enabled    =  true
Parameters for Kosloff-Kosloff__________________________________
  * RIGHT:
    boundary location      =  1000
    span of sponge layer   =  50
    range of sponge layer  =  [950, 1000]
  * BOTTOM:
    boundary location      =  6.3705e+06
    span of sponge layer   =  25
    range of sponge layer  =  [6.37052e+06, 6.3705e+06]
  * Expression for γ-factor:
    in solid: 1.1 / T0 * (VS / VP)^2 * exp(-0.04 * SPAN / (VP * T0))
    in fluid: 0.88 / T0 * exp(-0.04 * SPAN / (VP * T0))
================================================================

============================ Nr(s,z) ===========================
type   =  CONSTANT
value  =  50
================================================================

====================== Computed Nr on Mesh =====================
min Nr  =  5
max Nr  =  50
ave Nr  =  48
sum Nr  =  22096
* Nr has been limited by inplane resolution.
* Nr has been rounded up to FFTW lucky numbers.
================================================================

=========================== 3D Models ==========================
SEG_C3_volumetric ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  class name   =  StructuredGridV3D
  NetCDF file  =  ocean_topo_data/ocean_topo_fluid.nc
  Coordinates___________________________________________________
      COORD | NC-VAR | SCOPE
    * x     | x      | [-1050, 1050]
    * y     | y      | [-1050, 1050]
    * depth | depth  | [0, 510]
    depth below solid surface      =  false
    force element center in scope  =  true
  Properties____________________________________________________
      KEY | NC-VAR | REF KIND | RANGE
    * VP  | VP     | ABS      | [1.5e+06, 1.53e+06]
    * RHO | RHO    | ABS      | [1029, 1050]
    leader-only storage  =  true
SEG_C3_geometric ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  class name   =  StructuredGridG3D
  NetCDF file  =  ocean_topo_data/ocean_topo_undulation.nc
  Coordinates___________________________________________________
      COORD | NC-VAR | SCOPE
    * x     | x      | [-1050, 1050]
    * y     | y      | [-1050, 1050]
    * depth | N/A    | [140, 260]
    interface depth            =  200
    depth below solid surface  =  false
  Undulation data_______________________________________________
    NetCDF variable      =  depth
    data range           =  [-49.9497, 49.9497]
    leader-only storage  =  true
================================================================

=========================== Time Step ==========================
Δt determined by mesh  =  1.575e-07
   Courant number      =  0.5
   location (s,z)      =  [982.143, 6.37078e+06]
Δt enforced by user    =  NONE
Δt to be used          =  1.575e-07
================================================================

========================== Attenuation =========================
* Attenuation is turned off.
================================================================

===================== Computational Domain =====================
GLL points______________________________________________________
  FluidPoint$Mass1D              =  4989
  FluidPoint$Mass3D              =  1904
  Σ                              =  6893
Spectral elements_______________________________________________
  FluidElement$Acoustic1D        =  308
  FluidElement$PRT3D$Acoustic1D  =  4
  FluidElement$PRT3D$Acoustic3D  =  108
  Σ                              =  420
Axial boundary__________________________________________________
  Solid                          =  0
  Fluid                          =  61
  Σ                              =  61
Absorbing boundary______________________________________________
  ClaytonFluid1D                 =  194
  ClaytonFluid3D                 =  19
  SpongeFluid                    =  687
  Σ                              =  900
Fluid surface boundary__________________________________________
  Σ                              =  113
Solid-fluid boundary____________________________________________
  Σ                              =  0
Domain decomposition____________________________________________
  # sub-domains (nproc)          =  4
  # rank-to-rank communications  =  5
================================================================

!!!!!!!!!!!! AxiSEM3D ABORTED UPON RUNTIME EXCEPTION !!!!!!!!!!!
FROM: Unknown
WHAT: map::at:  key not found
TIME: 2021-07-23T23:43:06
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!! AxiSEM3D ABORTED UPON RUNTIME EXCEPTION !!!!!!!!!!!
FROM: Source::release
WHAT: Error locating source in mesh.
      Source index = 0
      Source name  = the_only_source
TIME: 2021-07-23T23:43:06
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!! AxiSEM3D ABORTED UPON RUNTIME EXCEPTION !!!!!!!!!!!
FROM: Unknown
WHAT: map::at:  key not found
TIME: 2021-07-23T23:43:06
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 3 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
!!!!!!!!!!!! AxiSEM3D ABORTED UPON RUNTIME EXCEPTION !!!!!!!!!!!
FROM: Unknown
WHAT: map::at:  key not found
TIME: 2021-07-23T23:43:06
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

[MBP:23069] PMIX ERROR: UNREACHABLE in file server/pmix_server.c at line 2198
[MBP:23069] 3 more processes have sent help message help-mpi-api.txt / mpi-abort
[MBP:23069] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

$$

This is our input directory: 04_Cartesian_ocean_topo_20Hz.zip

[kuangdai]

You cannot add a moment tensor in fluid. Use FLUID_PRESSURE instead. Unit can be considered as Pascal.

mechanism:
            # what: type of source mechanism
            # type: string
            # only: MOMENT_TENSOR, FORCE_VECTOR, FLUID_PRESSURE
            type: FLUID_PRESSURE
            # what: data for the source mechanism
            # type: array of double
            # note: 1) use [M11, M22, M33, M12, M13, M23] for MOMENT_TENSOR;
            #              [F1, F2, F3] for FORCE_VECTOR;
            #              [P] for FLUID_PRESSURE,
            #          where 123 stands for ZRT (vertical, radial, transpose)
            #       2) if horizontal location is given by "latitude_longitude",
            #          the RT-axes are determined w.r.t. the north pole;
            #          the moment tensor of an earthquake then follows the same
            #          order as in the CMTSOLUTION format (globalcmt.org)
            #       3) if horizontal location is given by "distance_azimuth",
            #          the RT-axes are determined w.r.t. the source (mesh axis)
            data: [1e20]
            # what: unit of data
            # type: double
            # note: use 1e-7 to convert dyn*cm (in CMTSOLUTION) to N*m
            unit: 1e-7

[kuangdai]

Mind you. Are you sure you are adding an undulating interface between two fluid layers? I cannot see the physical background but maybe you have your own reasons.

[tnissen]

Many applications can exist for such cases, most prominently the entire seismic industry: almost all of their production runs are still based on acoustics. Ideally, there would be no assumption from our developmental side as to what scenario might be plausible or not, as in most likelihood our own imagination is the limitation.

Tarje

On 26/07/2021 15:29, Kuangdai Leng wrote:

Mind you. Are you sure you are adding an undulating interface between two fluid layers? I cannot see the physical background but maybe you have your own reasons.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/kuangdai/AxiSEM-3D/issues/20#issuecomment-886704579, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABUMH5FQOC37KNHSS6JJRKDTZVPKDANCNFSM5A42A52Q.

-- <>--<>--<>--<>--<>--<>--<>--<>--<>--<>--<>--<> Associate Professor of Geophysics Dept. of Earth Sciences, Oxford University South Parks Road, Oxford OX1 3AN; UK www.earth.ox.ac.uk/people/tarje-nissen-meyer <>--<>--<>--<>--<>--<>--<>--<>--<>--<>--<>--<>

pronouns: he/him/his

[Yishharu]

Thank you Kuangdai and Tarje for the answer. And yes maybe the undulating boundary between fluid is not necessary. The stimulation now runs successfully. However all the recievers returns 0 in each traces. Wondering is there anything I'm setting wrong?

Many thanks!

slurm-output.txt input.zip

[Yishharu]

Hi Kuangdai,

We've managed to simulate a marine seismic shot in the cartesian mesh. This is the animation of 2s wave propagation in the fluid domain at one azimuthal slice. However, the left boundary looks a bit strange. The left boundary should be an unphysical boundary but seems a lot of energy is accumulated there. What do you think?

https://user-images.githubusercontent.com/13760732/131660458-9850fcec-92f3-4799-8f9a-cdd7ec485d8f.mp4

input.zip

Best, Zhi

[kuangdai]

Can you deactivate the 3D model and use constant Nr=1, and redo the simulation and show the animation here?

[Yishharu]

This is the wavefield after we deactivate the 3D model and use constant Nr=1. The edge oscillation still exists.

The 1D model is also a fluid domain. ocean_topo.bm.txt

https://user-images.githubusercontent.com/13760732/132880160-13de472c-51d4-4bc0-9dc2-c7a8e7fd7d82.mp4

[kuangdai]

Reopen this

[Yishharu]

animation setting;

inparam.output copy.yaml.txt

[kuangdai]

I think there could be an issue with computing displacement near the axis in fluid. In fluid, we actually compute pressure, so this issue does not affect far-field displacement results.

Can you try

channels: [P]

to animate pressure instead of [U]? Let see if heal the problem.

[Yishharu]

Thanks Kuangdai. I thought I'm plotting [P] in these animation already? you could check in the inparam.output file.

[kuangdai]

I see. So let's check the next possibility. The wavefield may vary quickly near the axis, so using a downsampled mesh grid is probably not enough. Can you use FULL for GLL_points_one_edge instead of [0, 2, 4]?

Also, in the post-processing, you need to connect all 25 GLL points to make 16 quads for animation. The following example connect 9 GLL points to make 4 quads for animation. You need to change the bit in the screenshot.

https://github.com/kuangdai/AxiSEM-3D/blob/master/examples/03_Cartesian_SEG_EAGE_salt_5Hz/post_processing.ipynb

[Yishharu]

This is another color normalisation of pressure from 0 to 1e+11. But indeed if this does not the far-field displacement results, the result should be fine. Thanks

https://user-images.githubusercontent.com/13760732/134338213-1ad631df-8f54-465c-a513-a69da8802945.mp4

!

[kuangdai]

The wavefield looks quite fuzzy, reminding me of another possible issue.

What is the period of the mesh? For animation, we cannot use a delta function as source-time function and convolve and filter later, as we do for seismograms. For animation, we use mesh period as half_duration. Here you are using half_duration: 0.2, is this comparable to mesh period?

[Yishharu]

Yes this might be the issue. Our mesh period is 0.05s. indeed the 0.2 half duration is a bit longer. I will adjust this first. And in terms of the full GLL points, could you remind me of the mesh connectivity in that case, shall I modify the code as followed:

# connectivity list, shared by all slices
# with GLL_points_one_edge = Full,
# the GLL point layout on each element is
# 0--1--2--3--4
# |  |  |  |  |
# 5--6--7--8--9
# |  |  |  |  |
# 10--11--12--13--14
# |   |   |   |   |
# 15--16--17--18--19
# |   |   |   |   |
# 20--21--22--23--24

connectivity = []
for ielem in np.arange(nelem):
    start = ielem * 25
    connectivity.append([start + 0, start + 1, start + 6, start + 5])
    connectivity.append([start + 1, start + 2, start + 7, start + 6])
    connectivity.append([start + 2, start + 3, start + 8, start + 7])
    connectivity.append([start + 3, start + 4, start + 9, start + 8])
    connectivity.append([start + 5, start + 6, start + 11, start + 10])
    connectivity.append([start + 6, start + 7, start + 12, start + 11])
    connectivity.append([start + 7, start + 8, start + 13, start + 12])
    connectivity.append([start + 8, start + 9, start + 14, start + 13])
    connectivity.append([start + 10, start + 11, start + 16, start + 15])
    connectivity.append([start + 11, start + 12, start + 17, start + 16])
    connectivity.append([start + 12, start + 13, start + 18, start + 17])
    connectivity.append([start + 13, start + 14, start + 19, start + 18])
    connectivity.append([start + 15, start + 16, start + 21, start + 20])
    connectivity.append([start + 16, start + 17, start + 22, start + 21])
    connectivity.append([start + 17, start + 18, start + 23, start + 22])
    connectivity.append([start + 18, start + 19, start + 24, start + 23])

[kuangdai]

Yes, that's correct.

[Yishharu]

It reminds me that in the previous animation, we've already set the half duration to 0. In this new 3D animation (still 9 points per element), we change the half duration to 0.05s to match the meshing period (20Hz). The wavefield of pressure indeed looks clearer. But the oscillation at the left boundary still exists. So you think this left boundary issue will not affect the far-field displacement results?

https://user-images.githubusercontent.com/13760732/134664002-1476f99c-28d8-4ace-8237-0627229516b0.mp4

[kuangdai]

That looks better. I think this is a visualization problem. Seismograms should be correct.