adapted irregular mesh with lasif

[yjgao-gfz]

Hello, here I ran into one problem when using the irregular mesh from salvus smoothsem,
but get an error which is related to the domain identification:

/tmp/ipykernel_30620/1694013832.py in
1 events=lasif.api.list_events(lasif_root, output=True)
----> 2 submit_forward('1','buildmesh/japan_hindu.h5',events)

/tmp/ipykernel_30620/1399908083.py in submit_forward(iteration, modelfile, events)
41 event=event,
42 iteration=iteration,
---> 43 side_set="r1"
44 )
45 if i%5==0 and i!=0 :

~/Salvus/salvus_lasif_yjgao-main/LASIF_2.0/lasif/salvus_utils.py in create_salvus_forward_simulation(comm, event, iteration, mesh, side_set)
95 mesh, sim_set["absorbing_boundaries_in_km"]
96 )
---> 97 if not domain.is_global_domain():
98 absorbing = sc.boundary.Absorbing(
99 width_in_meters=comm.project.simulation_settings[

~/Salvus/salvus_lasif_yjgao-main/LASIF_2.0/lasif/domain.py in is_global_domain(self)
503 def is_global_domain(self):
504 if not self.is_read:
--> 505 self._read()
506
507 if self.is_global_mesh:

~/Salvus/salvus_lasif_yjgao-main/LASIF_2.0/lasif/domain.py in _read(self)
188 # the points on the boundary are arranged in a way that a proper
189 # polygon will be drawn.
--> 190 sorted_indices = self.get_sorted_edge_coords()
191 x, y, z = self.domain_edge_coords[np.append(sorted_indices, 0)].T
192 lats, lons, _ = xyz_to_lat_lon_radius(x[0], y[0], z[0])

~/Salvus/salvus_lasif_yjgao-main/LASIF_2.0/lasif/domain.py in get_sorted_edge_coords(self)
478 # start sorting with the first node
479 print(indices_sorted)
--> 480 indices_sorted[0] = 0
481 for i in range(num_edge_points)[1:]:
482 prev_idx = indices_sorted[i - 1]

IndexError: index 0 is out of bounds for axis 0 with size 0

Below is the mesh building:

sm = SmoothieSEM()
sm.basic.model = "csem"
sm.basic.min_period_in_seconds = 30.0
sm.basic.elements_per_wavelength = 2.0
sm.basic.number_of_lateral_elements = 24
sm.advanced.tensor_order = 1

sm.source.latitude = 40
sm.source.longitude = 140
sm.chunk.max_colatitude = 80.0
sm.spherical.min_radius = 4000.0
sm.topography.topography_file = "topography_earth2014_egm2008_lmax_10800.nc"
sm.topography.topography_varname = (
"topography_earth2014_egm2008_lmax_10800_lmax_4096"
)
sm.topography.moho_topography_file = "moho_topography_crust_1_0_egm2008.nc"
sm.topography.moho_topography_varname = (
"moho_topography_crust_1_0_egm2008_lmax_128"
)
sm.refinement.lateral_refinements = [
{
"theta_min": 40.0,
"theta_max": 70.0,
"r_min": 4500.0,
"phi_min": -15.0,
"phi_max": 15.0,
}
]

sm.source.azimuth = 290.0

from obspy.clients.fdsn import Client

source = sn.simple_config.source.seismology.SideSetMomentTensorPoint3D(
latitude=40,
longitude=140,
depth_in_m=120,
side_set_name="r1",
mrr=5.47e15,
mtt=-4.11e16,
mpp=3.56e16,
mrt=2.26e16,
mrp=-2.25e16,
mtp=1.92e16,
)

get USarray stations from iris

inv = Client("gfz").get_stations(
network="*", level="station", format="text",minlatitude=34.782,
maxlatitude=46.5322,minlongitude=66.3029, maxlongitude=88.2316)

inv2 = Client("iris").get_stations(
network="*", level="station", format="text",minlatitude=34.782,
maxlatitude=46.5322,minlongitude=66.3029, maxlongitude=88.2316)

receivers = sn.simple_config.receiver.seismology.parse(
inv, dimensions=3, fields=["displacement"]
)

receivers2 = sn.simple_config.receiver.seismology.parse(
inv2, dimensions=3, fields=["displacement"]
)

newreceivers=receivers+receivers2

event_collection = sn.EventCollection.from_sources(
sources=[source], receivers=newreceivers
)
sm.spherical.ellipticity = 0.0033528106647474805
from salvus.mesh.mask_generators import RayMaskGenerator
rmg = RayMaskGenerator(
event_collection,
phases=["P","pP","S","sS"],
number_of_points_per_ray=100,
distance_in_km=1000.0,
)

m= sm.create_mesh(mesh_processing_callback=rmg)
m.find_surface()
#print(m.side_sets)
m.get_element_nodes().shape
m.side_sets['r1']
m.map_nodal_fields_to_element_nodal()
m.attach_global_variable("reference_frame", "spherical")
m.attach_global_variable("rotation_convention", "Spherical")
m.write_h5('japan_hindu.h5')