Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unknown AxiSEM3D running issue and Error locating source in mesh. #20

Closed
Yishharu opened this issue Jul 23, 2021 · 18 comments
Closed

Unknown AxiSEM3D running issue and Error locating source in mesh. #20

Yishharu opened this issue Jul 23, 2021 · 18 comments

Comments

@Yishharu
Copy link

Yishharu commented Jul 23, 2021

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
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Collaborator

tnissen commented Jul 26, 2021 via email

@Yishharu
Copy link
Author

Yishharu commented Aug 7, 2021

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
Copy link
Author

Yishharu commented Sep 1, 2021

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?

animation.mp4

input.zip

Best,
Zhi

@kuangdai
Copy link
Contributor

kuangdai commented Sep 1, 2021

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

@Yishharu
Copy link
Author

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

animation.mp4

@kuangdai
Copy link
Contributor

Reopen this

@kuangdai kuangdai reopened this Sep 22, 2021
@Yishharu
Copy link
Author

animation setting;

inparam.output copy.yaml.txt

@kuangdai
Copy link
Contributor

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
Copy link
Author

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

@kuangdai
Copy link
Contributor

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

image

@Yishharu
Copy link
Author

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

animation.mp4

!

@kuangdai
Copy link
Contributor

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
Copy link
Author

Yishharu commented Sep 22, 2021

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
Copy link
Contributor

Yes, that's correct.

@Yishharu
Copy link
Author

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?

animation.mp4

@kuangdai
Copy link
Contributor

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants