max frequency resolution of the visp-roten-Diehl mesh

Hi @mafanasiev , just wanted to clarify a few things here. Hope it is not too trivial, again!

1. I am wondering why I get this error when increasing the frequency resolution to >1hz. If my suspicion is right, it may have to do with having too many thin or tightly-packed layers, or high vertical or horizontal complexity (e.g., topography, internal structure) or even contradictory refinement rules (e.g., tripling + doubling + constant in one layer)? Is it harder to smoothly coarsen/refine between layers — hence the error below when going to higher frequencies?

SNIPPET:

# ## Generate the mesh
reference_frequency = 2.0
rf = reference_frequency
evname = "Vallorcine"


mesh = lm.mesh_from_domain(
    domain=d,
    model=lm.MeshingProtocol(
        layered_model,
        intralayer_coarsening_policy=[
            lm.IntralayerVerticalRefine(refinement_type="doubling"),
            lm.IntralayerVerticalRefine(refinement_type="tripling"),
            lm.IntralayerConstant(),
        ],
        interlayer_coarsening_policy=[
            lm.InterlayerConstant(),
            lm.InterlayerDoubling(),
        ],
    ),
    mesh_resolution=sn.MeshResolution(reference_frequency=rf),
)



meshname = "mesh_{}Hz_blend_Rhone_models_28325_1030.h5".format(int(rf))
mesh.write_h5(meshname)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File salvus/mesh/layered_meshing/meshing_protocol/coarsening_policy/intralayer_vertical_refine.py:449, in salvus.mesh.layered_meshing.meshing_protocol.coarsening_policy.intralayer_vertical_refine.IntralayerVerticalRefine.basic_block()

File salvus/mesh/data_structures/mesh_block/mesh_block_utils.py:419, in salvus.mesh.data_structures.mesh_block.mesh_block_utils._smooth_three_dimensional_refinements()

ValueError: Smoothing refinements failed

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
Cell In[38], line 7
      3 rf = reference_frequency
      4 evname = "Vallorcine"
----> 7 mesh = lm.mesh_from_domain(
      8     domain=d,
      9     model=lm.MeshingProtocol(
     10         layered_model,
     11         intralayer_coarsening_policy=[
     12             lm.IntralayerVerticalRefine(refinement_type="doubling"),
     13             lm.IntralayerVerticalRefine(refinement_type="tripling"),
     14             lm.IntralayerConstant(),
     15         ],
     16         interlayer_coarsening_policy=[
     17             lm.InterlayerConstant(),
     18             lm.InterlayerDoubling(),
     19         ],
     20     ),
     21     mesh_resolution=sn.MeshResolution(reference_frequency=rf),
     22 )
     26 meshname = "mesh_{}Hz_blend_Rhone_models_28325_1030.h5".format(int(rf))
     27 mesh.write_h5(meshname)

File salvus/mesh/layered_meshing/detail/mesh_from_domain.py:1018, in salvus.mesh.layered_meshing.detail.mesh_from_domain.mesh_from_domain()

File salvus/mesh/layered_meshing/detail/mesh_from_domain.py:744, in salvus.mesh.layered_meshing.detail.mesh_from_domain.create_layered_mesh()

File salvus/mesh/layered_meshing/detail/mesh_from_domain.py:725, in genexpr()

File salvus/mesh/layered_meshing/detail/mesh_from_domain.py:688, in genexpr()

File salvus/mesh/layered_meshing/meshing_protocol/coarsening_policy/intralayer_detail.py:510, in salvus.mesh.layered_meshing.meshing_protocol.coarsening_policy.intralayer_detail.construct_mesh_block()

File salvus/mesh/layered_meshing/meshing_protocol/coarsening_policy/intralayer_vertical_refine.py:456, in salvus.mesh.layered_meshing.meshing_protocol.coarsening_policy.intralayer_vertical_refine.IntralayerVerticalRefine.basic_block()

TypeError: isinstance() argument 2 cannot be a parameterized generic

2. Another thing I do not fully get is the deeper crust-- below the blended model/topography. It seems that is is reading the Vs, Vp from Tobias' 3D crustal model (defined here as : bedrock), but in the mesh viz I don't see it very well. It is using Tobias' model, correct?

SNIPPET where I think Tobias model is read and included to the model--> then to the mesh:

bedrock = lm.material.elastic.Velocity.from_params(
    rho=2000,
    vp=visp_tools.read_velomod("/Users/maria/swiss_wavefields-main/data/velomod_P.out"),
    vs=visp_tools.read_velomod("/Users/maria/swiss_wavefields-main/data/velomod_S.out"),
)

3. It runs at 1Hz without problems, however the mesh on paraview doesn't look "sharp" enough, but maybe this is a plotting issue.

The version of your original nb.py that i am currently using: https://gitlab.seismo.ethz.ch/deep/visp/-/blob/main/rhone_visp_diehl_3D.ipynb

Any feedback appreciated as always :) Thank you!