SPECFEM / specfem3d

SPECFEM3D_Cartesian simulates acoustic (fluid), elastic (solid), coupled acoustic/elastic, poroelastic or seismic wave propagation in any type of conforming mesh of hexahedra (structured or not).
GNU General Public License v3.0
392 stars 223 forks source link

Question about creating external models with cubit, can't generate absorption boundary files #1710

Closed zuoxuan-de closed 2 weeks ago

zuoxuan-de commented 3 weeks ago

Hello guys, I'm creating a three-layer model containing an alleyway, where the alleyway is located in the middle layer, in addition, I manually set his free surface and other absorbing surfaces. But when I run cubit2SPECFEM3D, I have a problem with my absorbing_surface_file_bottom,free_or_absorbing_surface_file_zmax,absorbing_surface_file_ymin and other similar absorbing surface files are all empty. I saw your replies related to this in my previous question, I hope you can be generous with your suggestions, thanks a lot! @casarotti @homnath image This is the specific case of the lane in the second layer, which is empty in the center, without any medium, and has all free surfaces. image This is the free surface of the lane I set manually. image In order to set the absorbing layer on the upper surface, I set both the upper and lower interfaces (zmax,zmin) to face_bottom again.

Next is my specific python script:

#!/usr/bin/env python

import cubit
import cubit2specfem3d

import os
import sys

SEMoutput='MESH03'
CUBIToutput='MESH_GEOCUBIT03'

os.system('mkdir -p '+ SEMoutput)
os.system('mkdir -p '+ CUBIToutput)

# Simulating a coal seam with a roadway
# This geological structure is mainly divided into three layers, the upper and lower peripheral rocks and the coal seam,
# which contains two roadways, of which the roadway portion is not defined by any medium and is all free surface
cubit.cmd('reset')
# underlying rock
cubit.cmd('brick x 50 y 40 z 10')
cubit.cmd('vol 1 move x 25 y 20 z 5')

#################################################  intermediate seam   ##############################################

#Leftmost part of the coal seam in the y-direction, also the left wall of the left roadway
cubit.cmd('brick x 50 y 15 z 5')
cubit.cmd('vol 2 move x 25 y 7.5 z -2.5')

# #Rightmost part of the coal seam in the y-direction, also the right wall of the right roadway
cubit.cmd('brick x 50 y 15 z 5')
cubit.cmd('volume 3 move x 25 y 32.5 z -2.5')

# The top of the lane in the y-direction
cubit.cmd('brick x 5 y 10 z 5')
cubit.cmd('volume 4 move x 2.5 y 20 z -2.5')

# The end of the lane in the y-direction
cubit.cmd('brick x 5 y 10 z 5')
cubit.cmd('volume 5 move x 47.5 y 20 z -2.5')

####################################################  end ####################################################
# top layer of surrounding rock
cubit.cmd('brick x 50 y 40 z 10')
cubit.cmd('volume 6 move x 25 y 20 z -10')

cubit.cmd('merge all')
cubit.cmd('imprint all')

# Meshing the volumes
elementsize = 1.0

cubit.cmd('volume all size 1')
cubit.cmd('mesh volume all')

cubit.cmd('block 1001 add surface 24 28 38 42 46 49 ')
cubit.cmd('block 1001 name "face_topo"')

cubit.cmd('block 1002 add surface 1 32 ')
cubit.cmd('block 1002 name "face_bottom"')

cubit.cmd('block 1003 add surface 4 10 16 22 34 ')
cubit.cmd('block 1003 name "face_abs_xmin"')

cubit.cmd('block 1004 add surface 3 9 33')
cubit.cmd('block 1004 name "face_abs_ymin"')

cubit.cmd('block 1005 add surface 6 12 18 36  30 ')
cubit.cmd('block 1005 name "face_abs_xmax"')

cubit.cmd('block 1006 add surface 5 17 35')
cubit.cmd('block 1006 name "face_abs_ymax"')

from geocubitlib import cubit2specfem3d

#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
cubit.cmd('block 1 hex in volume 1')
cubit.cmd('block 1 name "elastic 1" ')        # elastic material region
cubit.cmd('block 1 attribute count 7')
cubit.cmd('block 1 attribute index 1 1')      # flag for material: 1 for 1. material
cubit.cmd('block 1 attribute index 2 3500')   # vp
cubit.cmd('block 1 attribute index 3 2000')   # vs
cubit.cmd('block 1 attribute index 4 2400')   # rho
cubit.cmd('block 1 attribute index 5 9999.0') # Q_kappa
cubit.cmd('block 1 attribute index 6 9999.0')  # Qmu
cubit.cmd('block 1 attribute index 7 0 ')      # anisotropy_flag

cubit.cmd('block 2 hex in volume 2')
cubit.cmd('block 2 name "elastic 2" ')        # elastic material region
cubit.cmd('block 2 attribute count 7')
cubit.cmd('block 2 attribute index 1 2')      # flag for material: 1 for 1. material
cubit.cmd('block 2 attribute index 2 1900')   # vp
cubit.cmd('block 2 attribute index 3 1100')   # vs
cubit.cmd('block 2 attribute index 4 1300')   # rho
cubit.cmd('block 2 attribute index 5 9999.0')  # Q_mu
cubit.cmd('block 2 attribute index 6 9999.0')  # Qmu
cubit.cmd('block 2 attribute index 7 0 ')      # anisotropy_flag

cubit.cmd('block 3 hex in volume 3')
cubit.cmd('block 3 name "elastic 3"')
cubit.cmd('block 3 attribute count 7')
cubit.cmd('block 3 attribute index 1 3 ')
cubit.cmd('block 3 attribute index 2  1900 ')
cubit.cmd('block 3 attribute index 3  1100 ')
cubit.cmd('block 3 attribute index 4  1300 ')
cubit.cmd('block 3 attribute index 5 9999')
cubit.cmd('block 3 attribute index 6 9999')
cubit.cmd('block 3 attribute index 7 0')

cubit.cmd('block 4 hex in volume 4')
cubit.cmd('block 4 name "elastic 4"')
cubit.cmd('block 4 attribute count 7')
cubit.cmd('block 4 attribute index 1 4 ')
cubit.cmd('block 4 attribute index 2 1900 ')
cubit.cmd('block 4 attribute index 3 1100 ')
cubit.cmd('block 4 attribute index 4 1300 ')
cubit.cmd('block 4 attribute index 5 9999')
cubit.cmd('block 4 attribute index 6 9999')
cubit.cmd('block 4 attribute index 7 0')

cubit.cmd('block 5 hex in volume 5')
cubit.cmd('block 5 name "elastic 5"')
cubit.cmd('block 5 attribute count 7')
cubit.cmd('block 5 attribute index 1 5 ')
cubit.cmd('block 5 attribute index 2 1900 ')
cubit.cmd('block 5 attribute index 3 1100 ')
cubit.cmd('block 5 attribute index 4 1300 ')
cubit.cmd('block 5 attribute index 5 9999')
cubit.cmd('block 5 attribute index 6 9999')
cubit.cmd('block 5 attribute index 7 0')

cubit.cmd('block 6 hex in volume 6')
cubit.cmd('block 6 name "elastic 6"')
cubit.cmd('block 6 attribute count 7')
cubit.cmd('block 6 attribute index 1 6 ')
cubit.cmd('block 6 attribute index 2 3500 ')
cubit.cmd('block 6 attribute index 3 2000 ')
cubit.cmd('block 6 attribute index 4 2400 ')
cubit.cmd('block 6 attribute index 5 9999')
cubit.cmd('block 6 attribute index 6 9999')
cubit.cmd('block 6 attribute index 7 0')

cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')

#### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT

os.system('mkdir -p MESH0703')
cubit2specfem3d.export2SPECFEM3D('MESH0703')

# all files needed by SCOTCH are now in directory MESH

One last question, can cubit determine the dimensions of the three dimensions of a hexahedral cell, similar to 110.5, instead of a square with side lengths of 1. Thank you again and I hope to hear back from you all sooner rather than later!

homnath commented 3 weeks ago

I manually set the material blocks and surface files (See the Cubit file) and used MeshAssist to convert the Exodus file to SPECFEM3D files. Here I define three blocks, because you have three layers. You still need to create nummaterial_velocity_file. For the top surface, you have option within the Par_file whether to set it free or absorbing surface.

Best, Hom Nath


Hom Nath Gharti, PhD Assistant Professor | Digital Earth Scientist https://www.digitalearthscience.com/ Department of Geological Sciences and Geological Engineering Miller Hall 314, 36 Union St Queen’s University Kingston, ON K7L 3N6, Canada

Queen’s University is situated on traditional Anishinaabe and Haudenosaunee Territory.https://www.queensu.ca/encyclopedia/t/traditional-territories


From: zuoxuan-de @.> Sent: Thursday, July 4, 2024 4:32 AM To: SPECFEM/specfem3d @.> Cc: Hom Nath Gharti @.>; Mention @.> Subject: [SPECFEM/specfem3d] Question about creating external models with cubit, can't generate absorption boundary files (Issue #1710)

Hello guys, I'm creating a three-layer model containing an alleyway, where the alleyway is located in the middle layer, in addition, I manually set his free surface and other absorbing surfaces. But when I run cubit2SPECFEM3D, I have a problem with my absorbing_surface_file_bottom,free_or_absorbing_surface_file_zmax,absorbing_surface_file_ymin and other similar absorbing surface files are all empty. I saw your replies related to this in my previous question, I hope you can be generous with your suggestions, thanks a lot! @casarottihttps://github.com/casarotti @homnathhttps://github.com/homnath image.png (view on web)https://github.com/SPECFEM/specfem3d/assets/129966017/51e27655-b6ab-4dd6-b1d0-494a1dd681e2 This is the specific case of the lane in the second layer, which is empty in the center, without any medium, and has all free surfaces. image.png (view on web)https://github.com/SPECFEM/specfem3d/assets/129966017/16f66f0c-e2f5-4d8a-8304-b5e1cb31a673 This is the free surface of the lane I set manually. image.png (view on web)https://github.com/SPECFEM/specfem3d/assets/129966017/8324b787-823f-4efb-9919-db3e1a37c567 In order to set the absorbing layer on the upper surface, I set both the upper and lower interfaces (zmax,zmin) to face_bottom again.

Next is my specific python script:

!/usr/bin/env python

import cubit import cubit2specfem3d

import os import sys

SEMoutput='MESH03' CUBIToutput='MESH_GEOCUBIT03'

os.system('mkdir -p '+ SEMoutput) os.system('mkdir -p '+ CUBIToutput)

Simulating a coal seam with a roadway This geological structure is mainly divided into three layers, the upper and lower peripheral rocks and the coal seam, which contains two roadways, of which the roadway portion is not defined by any medium and is all free surface

cubit.cmd('reset')

underlying rock

cubit.cmd('brick x 50 y 40 z 10') cubit.cmd('vol 1 move x 25 y 20 z 5')

################################################# intermediate seam ##############################################

Leftmost part of the coal seam in the y-direction, also the left wall of the left roadway

cubit.cmd('brick x 50 y 15 z 5') cubit.cmd('vol 2 move x 25 y 7.5 z -2.5')

Rightmost part of the coal seam in the y-direction, also the right wall of the right roadway

cubit.cmd('brick x 50 y 15 z 5') cubit.cmd('volume 3 move x 25 y 32.5 z -2.5')

The top of the lane in the y-direction

cubit.cmd('brick x 5 y 10 z 5') cubit.cmd('volume 4 move x 2.5 y 20 z -2.5')

The end of the lane in the y-direction

cubit.cmd('brick x 5 y 10 z 5') cubit.cmd('volume 5 move x 47.5 y 20 z -2.5')

#################################################### end ####################################################

top layer of surrounding rock

cubit.cmd('brick x 50 y 40 z 10') cubit.cmd('volume 6 move x 25 y 20 z -10')

cubit.cmd('merge all') cubit.cmd('imprint all')

Meshing the volumes

elementsize = 1.0

cubit.cmd('volume all size 1') cubit.cmd('mesh volume all')

cubit.cmd('block 1001 add surface 24 28 38 42 46 49 ') cubit.cmd('block 1001 name "face_topo"')

cubit.cmd('block 1002 add surface 1 32 ') cubit.cmd('block 1002 name "face_bottom"')

cubit.cmd('block 1003 add surface 4 10 16 22 34 ') cubit.cmd('block 1003 name "face_abs_xmin"')

cubit.cmd('block 1004 add surface 3 9 33') cubit.cmd('block 1004 name "face_abs_ymin"')

cubit.cmd('block 1005 add surface 6 12 18 36 30 ') cubit.cmd('block 1005 name "face_abs_xmax"')

cubit.cmd('block 1006 add surface 5 17 35') cubit.cmd('block 1006 name "face_abs_ymax"')

from geocubitlib import cubit2specfem3d

Define material properties for the 3 volumes

cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################') cubit.cmd('block 1 hex in volume 1') cubit.cmd('block 1 name "elastic 1" ') # elastic material region cubit.cmd('block 1 attribute count 7') cubit.cmd('block 1 attribute index 1 1') # flag for material: 1 for 1. material cubit.cmd('block 1 attribute index 2 3500') # vp cubit.cmd('block 1 attribute index 3 2000') # vs cubit.cmd('block 1 attribute index 4 2400') # rho cubit.cmd('block 1 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 1 attribute index 6 9999.0') # Qmu cubit.cmd('block 1 attribute index 7 0 ') # anisotropy_flag

cubit.cmd('block 2 hex in volume 2') cubit.cmd('block 2 name "elastic 2" ') # elastic material region cubit.cmd('block 2 attribute count 7') cubit.cmd('block 2 attribute index 1 2') # flag for material: 1 for 1. material cubit.cmd('block 2 attribute index 2 1900') # vp cubit.cmd('block 2 attribute index 3 1100') # vs cubit.cmd('block 2 attribute index 4 1300') # rho cubit.cmd('block 2 attribute index 5 9999.0') # Q_mu cubit.cmd('block 2 attribute index 6 9999.0') # Qmu cubit.cmd('block 2 attribute index 7 0 ') # anisotropy_flag

cubit.cmd('block 3 hex in volume 3') cubit.cmd('block 3 name "elastic 3"') cubit.cmd('block 3 attribute count 7') cubit.cmd('block 3 attribute index 1 3 ') cubit.cmd('block 3 attribute index 2 1900 ') cubit.cmd('block 3 attribute index 3 1100 ') cubit.cmd('block 3 attribute index 4 1300 ') cubit.cmd('block 3 attribute index 5 9999') cubit.cmd('block 3 attribute index 6 9999') cubit.cmd('block 3 attribute index 7 0')

cubit.cmd('block 4 hex in volume 4') cubit.cmd('block 4 name "elastic 4"') cubit.cmd('block 4 attribute count 7') cubit.cmd('block 4 attribute index 1 4 ') cubit.cmd('block 4 attribute index 2 1900 ') cubit.cmd('block 4 attribute index 3 1100 ') cubit.cmd('block 4 attribute index 4 1300 ') cubit.cmd('block 4 attribute index 5 9999') cubit.cmd('block 4 attribute index 6 9999') cubit.cmd('block 4 attribute index 7 0')

cubit.cmd('block 5 hex in volume 5') cubit.cmd('block 5 name "elastic 5"') cubit.cmd('block 5 attribute count 7') cubit.cmd('block 5 attribute index 1 5 ') cubit.cmd('block 5 attribute index 2 1900 ') cubit.cmd('block 5 attribute index 3 1100 ') cubit.cmd('block 5 attribute index 4 1300 ') cubit.cmd('block 5 attribute index 5 9999') cubit.cmd('block 5 attribute index 6 9999') cubit.cmd('block 5 attribute index 7 0')

cubit.cmd('block 6 hex in volume 6') cubit.cmd('block 6 name "elastic 6"') cubit.cmd('block 6 attribute count 7') cubit.cmd('block 6 attribute index 1 6 ') cubit.cmd('block 6 attribute index 2 3500 ') cubit.cmd('block 6 attribute index 3 2000 ') cubit.cmd('block 6 attribute index 4 2400 ') cubit.cmd('block 6 attribute index 5 9999') cubit.cmd('block 6 attribute index 6 9999') cubit.cmd('block 6 attribute index 7 0')

cubit.cmd('export mesh "top.e" dimension 3 overwrite') cubit.cmd('save as "meshing.cub" overwrite')

Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT

os.system('mkdir -p MESH0703') cubit2specfem3d.export2SPECFEM3D('MESH0703')

all files needed by SCOTCH are now in directory MESH

One last question, can cubit determine the dimensions of the three dimensions of a hexahedral cell, similar to 110.5, instead of a square with side lengths of 1. Thank you again and I hope to hear back from you all sooner rather than later!

— Reply to this email directly, view it on GitHubhttps://github.com/SPECFEM/specfem3d/issues/1710, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ABMCQ4VHAPIOJJBJT65LRCTZKUCAFAVCNFSM6AAAAABKLB2Q4OVHI2DSMVQWIX3LMV43ASLTON2WKOZSGM4TAMRXGM3TINA. You are receiving this because you were mentioned.Message ID: @.***>

danielpeter commented 3 weeks ago

there are a few problems with your meshing script:

-- use an order of "imprint & merge": that is, use first imprint, then merge command. this will connect your volumes correctly, otherwise your volumes and mesh won't be conforming, but still be separated. thus, replace:

cubit.cmd('merge all')
cubit.cmd('imprint all')

with

cubit.cmd('imprint all')
cubit.cmd('merge all')

-- the blocks defined for the surfaces need to have faces, not surfaces added: thus, replace any occurrences of add surface with face in surface, like:

cubit.cmd('block 1002 add surface 1 32 ')

with

cubit.cmd('block 1002 face in surface 1 32 ')

this will then find all faces for those blocks.

-- replace your block name face_topo to something like face_free and use a different block number than 1001: the face_topo is used for a free surface that is also a topography surface at the top of a mesh. when outputting such topography surfaces, it will include a vertical normal check to make sure that the surface has a correct face point numbering.

since your "free" surface has vertical sides, you should avoid this check and define it as a more general free surface. thus, use a block name like face_free or free_surface (anything with a free added will do) and use a block number other than 1001 which is reserved for topography surfaces again. thus, replace:

cubit.cmd('block 1001 add surface 24 28 38 42 46 49 ')
cubit.cmd('block 1001 name "face_topo"')

with

cubit.cmd('block 1000 face in surface 24 28 38 42 46 49 ')
cubit.cmd('block 1000 name "face_free"')