Converting a regional 3D model to a specfem mesh (Addition to docs?)

[lsawade]

Hi,

Opening this issue because I believe the documentation is lacking here and I would like to do this. @clairedd has an external regional, 3D model that is based on a different SEM mesh and can probably write an output routine that converts it to something else.

Questions

I would like to have answered before diving into writing my own routines:

• What is the best format for specfem interpolation routines?
  • Should the mesh maybe be created first, and then the points' model parameters interpolated using Lagrange interpolation from the original SEM mesh?
• For something like this in particular, which is the best example to look at for loading an external model?

@danielpeter do you have suggestions here?

Related

• What flags in the constants.h.in should I look out for? For one model that I'm working with I'm using, e.g., EARTH_REGIONAL_MOHO_MESH

@rdno do you have input on this?

[lsawade]

Adding @pmoulik and @clairedd to the discussion.

[bch0w]

Hi @lsawade, just saw this issue and wanted to chime in that this feature would be really beneficial to our group and is something we've had on our to-do list for a while, and is likely the last key piece that would enable us to fully migrate from Cartesian to Globe (in most cases).

Our overall goal is to be able to take NetCDF models from IRIS EMC and use them directly in a one-chunk regional cut off mesh. But, even understanding the best way to implement external 3D models (in any format) into the code would be super useful.

Have you made any progress on this or is this something you're working on currently? @thurinj spent some time on this problem. Adding in @carltape to the conversation as well

[lsawade]

Hi @bch0w , so far I have only chatted with @danielpeter, and for my solution the best way was going to be just following the standard procedure of making crust and mantle separately, and reading the model from text files.

Do note that I'm still very interested in this because I think it would truly change the way people would use specfem, and reduce the hurdle of bringing your own model. The problem really is the meshing when it comes to land and Moho topography when velocities are first assigned. Maybe a focused meeting could be very helpful where we all chime in on what models need and what the easiest thing is read a model onto a mesh.

• topography divided meshes regions? crust/upper mantle, upper mantle lower mantle
• prescribe mesh topography?
• No topography
• local 1D, 2D (?), and 3D models.

Should these things be completely customizable how are flags set etc. I'm totally willing to spend substantial time on this, but I need to know whether it's worthwhile, and what people want.

[bch0w]

Hi @lsawade, this is great, it's nice to know we have overlapping interests here.

To provide context for our needs, we have an Alaska-based point velocity model with regular grid spacing in netCDF format that we are interested in using for regional-scale simulations.

@thurinj has written some FORTRAN code to read netCDF format, but it seems like the difficult decision is how to handle parts of the background model not defined by the regional model, which I think is related to the problem of assigning parameters to topography (?).

We had some discussion on the optimal way to drop the regional model into a background model, e.g., by providing perturbations to a 1D/3D background model (like PREM or s20rts), or by extrapolating the regional model into any gaps in the background, but this is as far as we got to my understanding.

Just to clarify, are you also interested in bringing in external mesh attributes like topography and moho topography, or just point model values? We are mainly interested in the model values and will likely be okay with using one of the in-built topography files like ETOPO1.

I think a focused discussion would be great to see if we have the same goals and could work together on this. I am also willing to dedicate some time to this as I can see it being a very useful feature for SPECFEM moving forward.

Perhaps we can make this a discussion topic at the next monthly Zoom, and go from there.

[lsawade]

Hi @bch0w, yes totally. I'm happy to work on this together!

@thurinj has written some FORTRAN code to read netCDF format.

Is the netcdf choice is based on simpler abstractions for gridding, and/or because IRIS models come in netcdf? How is the Fortran netcdf API in terms of reading Fortran/C ordered arrays. I recently started writing a fortran API for the GF database, and my main struggle was handling array order of multidimensional arrays written with C and read with Fortran. In my case the simple fix was fully transposing the numpy arrays I was writing with h5py, but other than that I couldn't really find a "good" fix.

but it seems like the difficult decision is how to handle parts of the background model not defined by the regional model, which I think is related to the problem of assigning parameters to topography (?).

Exactly this is the main issue, for which it is so hard to come up with some default behavior in terms of parameters of an Earth model. For _globe and other heterogenous global models, so many decisions have been made in terms of topographies of the layers, ellipticity, physics involved making it difficult to create a single fileformat.

I think what we are really missing is a way to handle these things. Especially in a way that can take care of future changes, e.g., converting the model is in simple anisotropy (vpv, vph, etc.) to full anisotropy (c_{ij}) without losing generality on meshing. So that the models can be added to later for higher complexity, without losing generality in the file format.

One way is to simply include important feature such as ETOPO as it is used in a model in the file, so that the mesher can refer to it later on. For the GF database with _globe for example, I store the bathymetry/topography that is used for _globe in the database files, because I need the topography for the locating the source in terms of depth. I don't think we have to be worried about bloating the files either because the topographies are peanuts with respect to model itself.

We had some discussion on the optimal way to drop the regional model into a background model, e.g., by providing perturbations to a 1D/3D background model (like PREM or s20rts), or by extrapolating the regional model into any gaps in the background, but this is as far as we got to my understanding.

I know this is very common in seismology to provide models in perturbations like that, but it becomes super weird the moment you add a different topography to your model. Say, your old model that you refer to is computed with CRUST2.0 and you compute your new model with CRUST1.0 then because of the difference in MOHO topographies in the two models the perturbations around the Moho are just "silly". Don't get me wrong I think perturbations are important for the interpretation of model differences, but may be ill-suited for representation in a standardized 3D model format. I'm happy to be corrected here!

Just to clarify, are you also interested in bringing in external mesh attributes like topography and moho topography, or just point model values?

Yeah, so I think that ideally most decisions that were made should be stored in a single file especially when it comes to topographies -- impact can be large here. [but maybe flags like gravity, rotation, etc. as well] The question is how would we translate this to the mesher. Here, we definitely need to discuss with Daniel since we don't want to reinvent the wheel.

I think a focused discussion would be great to see if we have the same goals and could work together on this. I am also willing to dedicate some time to this as I can see it being a very useful feature for SPECFEM moving forward.
Perhaps we can make this a discussion topic at the next monthly Zoom, and go from there.

Sounds perfect.

[bch0w]

This is all great! Thanks for the details and happy to advance this further during the Zoom next week. The choice for netCDF is primarily motivated by the availability at the EMC, and having it act as a central repository/ common file format for storing regional models, both SEM-derived and not. I will have to talk to @thurinj about the intricacies of the FORTRAN code.

Just wanted to throw in that I share the opinion that the perturbation is not the ideal way forward, I think our discussions were mainly focused on "how is specfem globe handling this right now?", which is how perturbations came up. We would definitely prefer to move away from this requirement to a more flexible framework for dropping in 3D models.

Now excuse me while I revisit my childhood and re-watch The Road to El Dorado

[lsawade]

Now excuse me while I revisit my childhood and re-watch The Road to El Dorado

A must.

[thurinj]

Here is the bare-bones piece of code that @bch0w mentioned:
specfem_nc_tomo

It is not much, but I hope it might kickstart things.

[SeismoFelix]

Thanks, everyone for addressing this. @thurinj the link you shared seems not to work. Probably I am well behind you on this but if this work it would be a huge leap in my career plan. So, more than happy to put my modest hands to help you with this.

[thurinj]

Kudos for pointing this out @SeismoFelix ! And now thanks to you the repo is even public haha 😅!

[danielpeter]

as discussed in out last zoom developer meeting, i committed these two PRs that add an initial version to read in (isotropic) IRIS EMC models:

• updates ADIOS2 and LIBXSMM library calls; adds initial IRIS EMC model support #818_
• adds scaling for EMC models #819_

it's intended as a starting point on how to use such external models - and hopefully we can iterate more on this.

one particular problem with these EMC models is that it was unclear if they were given with reference to a spherical earth. that is, SPECFEM3D GLOBE is used to read in mantle models defined for a spherical earth. once those model values are read in, the mesh will stretch out to accomodate topography and ellipticity.

for these EMC models, i assumed the opposite in that they defined the grid points at actual depth positions - including topography, but without the ellipticity factor. ellipticity is probably less of an issue for local models, but spherical reference would be important to know. unfortunately, based on the EMC format alone, one doesn't obtain that info. so for now, any EMC model is read in for GLL mesh point positions that are stretched out already to accommodate surface topography in case that is turned on in the 'Par_file'. and, ellipticity is added afterwards if selected.

[lsawade]

@danielpeter

do you think it is possible to read an EMC model different model values for the different regions. Say, the model has a parameter using_crust_version = 1.0, and arrays vpv_crust, vph_crust, vpv_mantle, vph_mantle, vsv_crust, vsh_crust, vsv_mantle, vsh_mantle.

How hard would it be to implement meshing with crust[1,2] and topography and the external model? And, of course the population of the different regions with the model values?

[danielpeter]

hi Lucas,

having an EMC model with velocity arrays separated by crust/mantle regions would be possible to implement with probably a few changes. at the moment, the EMC_model implementation is really for simple models only, i.e., only supports isotropic models.

the last question I didn't understand fully. you mean how to change the mesh with a Crust1.0 or Crust2.0 as crustal reference? at the moment, you can just add an ending _crust1.0 to switch from the default Crust2.0 to a Crust1.0 mesh. for these EMC models, this would mean setting the model as

MODEL               = EMC_model_crust1.0

for now, there seems to be no EMC model with different velocity arrays defined for different crust/mantle regions. all EMC models come with some single arrays, defined for the whole model range:
http://ds.iris.edu/ds/products/emc-earthmodels/

I couldn't find a model with a metadata entry like using_crust_version either. this would need to be added as well in our EMC_model implementation to support it.

[lsawade]

Hi Daniel,

I was mainly thinking about what the best idea is to implement Moho topography from external models, and I forgot that meshing is just done with _crust<1.0,2.0>, so forget about the second question. So, having separate parameters for vpv_crust and vpv_mantle is then probably not crazy to implement. There would have to be an extensive parameter grid check then, unless the user has to add an _transversely_isotropic to the model which would trigger transversely isotropic PREM to be used to populate the mesh, and in model_EMC we would look for vpv, vph, etc. What I am asking is ,does something like this make sense:

IF model == model_EMC

• no crust mesh
• isotropic PREM
• look for vp, vs, rho

ELSEIF model == model_EMC_crust<1.0,2.0>

• crust mesh
• isotropic PREM
• look for vp_crust, vp_mantle, vs_crust, vs_mantle, rho_crust, rho_mantle

ELSEIF model == model_EMC_transversely_isotropic

• no crust mesh
• transversely isotropic PREM
• look for vpv, vph, vsv, vsh, rho, or kappa, betav, betah?

ELSEIF model == model_EMC_transversely_isotropic_crust<1.0,2.0>

• crust mesh
• transversely isotropic PREM
• look for vpv_<crust,mantle>, vph_<crust,mantle>, vsv_<crust,mantle>, vsh_<crust,mantle>, rho_<crust,mantle>, or kappa_<crust,mantle>, betav_<crust,mantle>, betah_<crust,mantle>

Also, just FYI, I have embedded/meshed @clairdd 's model CANVAS into PREM. It is not in the EMC just yet, because they are finishing the manuscript (preliminary poster). I have not yet tried waveform modeling because we want to try to avoid as much as we possible can boundaries to isotropic PREM. Hence, @clairedd gave me a larger area of her model to work with, I just haven't gotten to it. What I wrote above is all due to the fact that her model is, in fact, transversely isotropic and does use crust1.0 in the mesh, she just gave me the isotropic non crust version for now so that it works with the model_EMC module.
I'm willing to try implementing transverse_isotropy in the model_EMC.F90 unless you prefer doing it yourself. All I'd need is the "ok" to above on how to implement it.

[danielpeter]

hi Lucas, yes sure please go ahead! the model_EMC file is there to iterate on, so make any changes as you see fit.

for this CANVAS model, I would use a model name like

MODEL   =  EMC_model_CANVAS

CANVAS is using the SPiRaL model as an initial starting point. I'm not sure how they used this model in Salvus (topography, internal discontinuities, etc.?), but we have SPiRaL implemented as well, which comes with its own crust definition. so, to overimpose this new EMC model values onto the SPiRaL one, it might be easier to just use a new name like above and then have the same reference flags as SPiRaL set in the get_model_parameters_flags() routine.

regarding the separate crust/mantle arrays, does CANVAS really need those or could it provide just the vpv,vph,.. velocities as single arrays for the whole model?

separate arrays for crust/mantle only make sense, if you want to model a first-order discontinuity with a jump of velocities at an internal interface, like the Moho. at longer seismic periods (CANVAS goes down to ~20 s periods), this might also just be visible as a strong velocity gradient in the model, which can be specified by a single array. again, this might depend on how they implemented SPiRaL in Salvus and how it was updated for CANVAS - well, just see what representation is really needed.

you might find that for separated arrays, you will need an additional element info (if it is located above/below the interface, i.e., in the crust or mantle) before overimposing the model velocities. at the moment, this is not done in the model_EMC_crustmantle() routine.

let me know how it goes with the implementation, happy to help out in case needed :)

[lsawade]

Ok, will work on it. I have so far only 'played' around with the mesher, even output some vtk files that all seem fine. At runtime however, I encountered an allocation issue even though I'm using a regional cutoff:

At line 221 of file src/specfem3D/prepare_gpu.f90
Fortran runtime error: Allocatable actual argument 'abs_boundary_ispec_outer_core' is not allocated

cutoff is at 400km. Am I overlooking something?

[lsawade]

After looking a lot at seismograms, I started looking more at the mesh, And it's mildly weird. I'm looking into what's happening, but it somehow looks like certain values are being misassigned:

Note that this is very much exaggerated. I used the following line in the paraview calculator

( 20 * ( sqrt (coordsX^2 + coordsY^2 + coordsZ^2) - 1) + 1)  * coords

wait, were you seeing lines like this, @bch0w ?

[lsawade]

Just to make sure that it's not just the vtk low-res plotting output, I also plotted a few high-res slices with all GLL points. This is looking really wonky. I think this is different from what you showed at the recent specfem meet, right, @bch0w ?

[lsawade]

Just an update, it looks as this is a bug, when I make the mesh very fine compared to the regional extent.

[bch0w]

Woops, sorry @lsawade, I though I had responded to this! This does not look like what I had showed at the last meeting, though. Seems like you've got a handle on what is causing this?

[SeismoFelix]

HI! Thanks for all the effort you have put on implementing routines for reading NetCDF external 3D velocity models. I tried to run the Alaska example, but unfortunately is not working for me.

• I was able to compile the devel branch version including the flags for dealing with NetCDF models.
• I downloaded the Alaska model and I created the symbolic link model.nc in the directory DATA/IRIS_EMC.
• I did not change the Par_file or any other file provided in the example.
• The mesher run xmeshfem3D fails when tries to read the model.nc file.
  The output_mesher.txt seems to be the same that the one provided in the example, but it does not continue when reaches this point:

broadcast model: EMC model
file : ./DATA/IRIS_EMC/model.nc

The output was:

running script...


Fri Mar 15 20:35:46 -05 2024
starting MPI mesher on 4 processors

 simulation time step size:
   DT determined =   0.11500000000000000     
   Par_file: user overrides with specified DT =    5.0000000000000003E-002

 using local mesh layout:
   number of layers in crust  =            4
   number of layers in mantle =            8

   number of doubling layers  =            2
            doubling layer at =            2
            doubling layer at =            5
   fictitious moho at depth   :    40.0000000     (km)


Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:

Could not print backtrace: executable file is not an executable
#0  0x10dd0faae
#1  0x10dd0ec8d
#2  0x7ff807381c1c
#3  0x10e5d32f1
#4  0x10e573ff3
#5  0x10e5eca30
#6  0x10e39500a
#7  0x10e423378
#8  0x10da0bd43
#9  0x10da12b73
#10  0x10d9c433f
#11  0x10d9c5a6e
#12  0x10d9c596a
#13  0x10d28b92a
#14  0x10d32fa15
#15  0x10d0433b1
#16  0x10d04be6a
#17  0x10cf7b51a
#18  0x10cf7c9a4
#19  0x10d01900a
#20  0x10d04d168
#21  0x7ff80702430f
--------------------------------------------------------------------------
prterun noticed that process rank 0 with PID 5927 on node Felixs-MBP exited on
signal 8 (Floating point exception: 8).

• To double check that this was an error of the NetCDF example and not coming from Specfem devel branch instalation I run another example (the one for the Middle East) and worked as expected.

I would appreciate any comment you may have about the origin of this failure.

Thanks a lot!

[danielpeter]

Felix, could you re-compile the code with debug flags added, -g -fbacktrace if you use a gfortran compiler or something like -g -traceback for ifort, to get something readable for the backtrace?

also, the latest devel should fix an intel compiler issue with EMC models, make sure to have this latest devel version in case.

[SeismoFelix]

Thanks for the suggestions @danielpeter . I added the debug flags when running configure for making the makefile:
FCFLAGS="-g -fbacktrace" ./configure FC=gfortran CC=gcc MPIFC=mpif90 --with-netcdf NETCDF_INC=/usr/local/include NETCDF_LIBS="-L/usr/local/lib -lnetcdff"

The make file was made without any issue including the debug flags:

FCFLAGS = -g -fbacktrace -I/usr/local/include

and the make all instruction worked as well.

However, the output when running the EMC example shows the same error I already posted. It does not show any additional information when crashes and I still get the same message:

Backtrace for this error: Could not print backtrace: executable file is not an executable

I also downloaded the latest devel version. I got this directory: specfem3d_globe-devel

Thanks for your help.

[danielpeter]

do you have more infos to reproduce the error, like your OS & compiler versions?

regarding the debugging options, you could also add the flag --enable-debug to the ./configure .. command to switch on more debugging checks:

./configure --enable-debug FC=gfortran CC=gcc MPIFC=mpif90 --with-netcdf NETCDF_INC=.. NETCDF_LIBS=..

maybe that triggers something readable.

[lsawade]

Hi @SeismoFelix, hi @ammcpherson, hi @danielpeter

So I tried the recent updates and it seems to work for me for the Alaska model.

@SeismoFelix, I have now tried again to input CANVAS (@clairedd's models), and got a SIGFPE when the _FillValue was a NaN, because the specfem fill value expects a float. So, replacing the _FillValue with a 9999.0 fixed that. That being said the traceback revealed that it happened much later, than what you are experiencing. For a bit I struggled with ieee_arithmic fortran module not being available on the cluster that I was working on, but that should appear at compilation and not runtime. Switching systems worked.

@danielpeter, and @ammcpherson, a few notes and questions,

• --enable-debug has a hard time compiling due to a bunch of conversion errors in src/shared/fft.f90. At least with the GCC compiler that I have been working on.
• When trying to run all of this with ADIOS2, the stacey boundary condition is failing to write. I did not look into it, because I just wanted to run some tests. I do not know whether this is on my end, or whether there is no implementation for stacey in ADIOS2.
• Thanks for implementing the manual crustal structure settings. When I started trying to get CANVAS into SPECFEM, I edited the part of specfem that setups the layers manually. This is much better ❤️
• Also, the model looks great already, awesome work, so thanks, @danielpeter, and thanks, @ammcpherson! There are just a few places where we (Claire and I) really need to be exact in terms of what their simulation is doing and what ours is doing.
• Attenuation and reference frequency. This morning I chatted with JT about all of this, and he made a good point in terms of dispersion. I assume that PREM's values are expected to be with respect to 1s, but input velocity model values are probably not. So, the velocity models values are with respect to -- most likely -- the highest resolvable frequency in the model, and these are then corrected to another velocity using the standard linear solids Q. Do you think there is a smart way to include reference frequency here in model_EMC? @clairedd is also checking with Michael Afansiev, to see how the the velocities are "picked" from the model.

Thanks already! I think I'm starting to get somewhere! Waveforms still look pretty different, but I ran many tests now, and just turning on and off some of the flags in specfem (topography, rotation, attenuation, ocean, etc.) has such an enormous effect on the waveforms at 20-50s, some better some worse, that now it's a question of converging (hopefully).

Cheers,

Lucas