Allow multiple external tomography files

bch0w commented 1 year ago

What does this PR do?

Allows Users to specify an arbitrary number of external tomography models
Introduces a new Par_file parameter: NFILES_TOMO
NFILES_TOMO tells SPECFEM to look for an integer number of tomography files in TOMOGRAPHY_PATH (valid range [1, inf])
Files require explicit naming as: tomography_model_01.xyz, tomography_model_02.xyz, ..., tomographymodel.xyz
Note the leading 0 (assuming <99 external models)
Adjusts SPECFEM machinery to accommodate the new parameter

Why was it initiated?

Related Issue: #1306

When importing external tomography models, fine spatial sampling at depth leads to unnecessarily large ASCII files. So we (@carltape et al.) split up our domains by depth regions (e.g., shallow, crust, mantle).

SPECFEM in it's current state can only accommodate one tomography model, so to get around this @rmodrak wrote up a small hardcoded edit to the source code to allow for 3 tomography models. Originally posted by @rmodrak in https://github.com/SPECFEM/specfem3d/issues/1306#issuecomment-435985634

Following @rmodrak I made this edit more generalized by allowing for a User-defined parameter that allow for an arbitrary number of tomography files to be included. By default this parameter is set to 1, which is the same behavior as before.

Likely this will need some code review. I'm not sure how far-reaching these changes are for other model types or for CUBIT meshes. Although I did test with my own 3-file external tomography model using the internal mesher and was able to reproduce correct results that I was getting when using the edited source file from @rmodrak

I'm open to other ways of implementing this that might not require a new Par_file parameter. Examples will need to be updated to deal with the new Par_file.

Here's relevant lines from output_generate_databases.txt:

   ...determining velocity model                                                 

      number of tomographic models       =            3                          
      maximum number of data records     =     23402520                          
      size of required tomography arrays =    1071.282     MB per process        

      material id:           -1                                                  
      file       : DATA/tomo_files/tomography_model_01.xyz                       
      data format: #x #y #z #vp #vs #density #Q_p #Q_s                           
      number of grid points = NX*NY*NZ:     1656720                              

      material id:           -2                                                  
      file       : DATA/tomo_files/tomography_model_02.xyz                       
      data format: #x #y #z #vp #vs #density #Q_p #Q_s                           
      number of grid points = NX*NY*NZ:    23402520                              

      material id:           -3                                                  
      file       : DATA/tomo_files/tomography_model_03.xyz                       
      data format: #x #y #z #vp #vs #density #Q_p #Q_s                           
      number of grid points = NX*NY*NZ:    14157080

And a figure confirming that this model was partitioned correctly to my mesh (New Zealand context).

nzatom_north_vs

output_generate_databases.txt

buildbot-princeton commented 1 year ago

The files changed require a buildbot test. Buildbot has started.

buildbot-princeton commented 1 year ago

The files changed require a buildbot test. Buildbot has started.

buildbot-princeton commented 1 year ago

Your changes failed the buildbot test. See build 934. Please fix the problems and push your changes back to your branch bch0w/specfem3d:feature-multiple_tomo_files. This pull-request will be reopened automatically. You don't need to create a new pull request

buildbot-princeton commented 1 year ago

Your changes failed the buildbot test. See build 935. Please fix the problems and push your changes back to your branch bch0w/specfem3d:feature-multiple_tomo_files. This pull-request will be reopened automatically. You don't need to create a new pull request

danielpeter commented 1 year ago

hi Bryant,

not sure if these additional options in Mesh_Par_file are really necessary.

NFILES_TOMO is more of a derived quantity, based on the input materials given. the current model_tomography.f90 routines work for multiple resolution tomography files, when defined in nummaterial_velocity_file using external meshes (CUBIT/Trelis,..). in nummaterial_velocity_file, you can have multiple lines, like:

2  -1 tomography elastic tomography_model_1.xyz 1
2  -2 tomography elastic tomography_model_2.xyz 2
2  -3 tomography elastic tomography_model_3.xyz 3
..

and corresponding material ids in the mesh defined to use different tomography files with different resolution.

in your case, you want to make this work for meshes generated by the in-house mesher xmeshfem3D? there, we defined the tomography materials as:

# number of materials
NMATERIALS                      = 1
# define the different materials in the model as :
# #material_id  #rho  #vp  #vs  #Q_Kappa  #Q_mu  #anisotropy_flag  #domain_id
#     Q_Kappa          : Q_Kappa attenuation quality factor
#     Q_mu             : Q_mu attenuation quality factor
#     anisotropy_flag  : 0 = no anisotropy / 1,2,... check the implementation in file aniso_model.f90
#     domain_id        : 1 = acoustic / 2 = elastic
# for tomography:
#   #material_id #type-keyword #domain-name #tomo-filename #tomo_id #domain_id
#   example:
#   -1 tomography elastic tomography_model.xyz 0 2
#    
-1 tomography elastic tomography_model.xyz 0 2

this one seems to work for only a single tomo file so far. we could just add multiple materials like this, together with the corresponding entries in the Domains section:

# number of materials
NMATERIALS                      = 3
# define the different materials in the model 
-1 tomography elastic tomography_model_1.xyz 0 2
-2 tomography elastic tomography_model_2.xyz 1 2
-3 tomography elastic tomography_model_3.xyz 2 2

# number of regions
NREGIONS                        = 3
# define the different regions of the model as :
#NEX_XI_BEGIN  #NEX_XI_END  #NEX_ETA_BEGIN  #NEX_ETA_END  #NZ_BEGIN #NZ_END  #material_id
1              24            1               24             1         4       -1
1              24            1               24             5         8       -2
1              24            1               24             9         12      -3

this would cover the handling with multiple resolution tomography files for the in-house mesher, without changing the input file format. what do you think?

if we can agree on this handling, then let me update the mesher to support this quickly.

bch0w commented 1 year ago

Thanks for your input @danielpeter! Yes, the idea would be to allow multiple external models with the internal mesher, similar to how external mesher can do it. I definitely kludged a few things to get this PR working so I'm glad to have your oversight.

I think your proposed idea is great, allowing for arbitrary file naming, and an intermixing of external models and internally defined properties (related to @carltape comment (via personal email), "Can you mix a tomography model with a user-defined homogeneous property?").

If it's easy enough for you to implement quickly that would be simplest. I'd also be happy to take a crack at it to improve my SPECFEM/FORTRAN skills; I'd just need to know which files to look at and could open a new PR for further code review. Either way, we'll be glad to have this feature moving forward!

danielpeter commented 1 year ago

will close this one and submit a new PR addressing this feature.

SPECFEM / specfem3d

Allow multiple external tomography files #1558

What does this PR do?

Why was it initiated?