search

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 & unstructured).
https://specfem.org
GNU General Public License v3.0
416 stars 231 forks source link

Problems using cubit to generate cpml model files for forward simulation #1770

Open zuoxuan-de opened 1 day ago

zuoxuan-de commented 1 day ago

Hi guys, I've recently been trying to use cubit to generate cpml model files for forward simulation, but some times it works, but when I change the model it doesn't, there's a specific model creation process below, but no matter how much I try it, it doesn't work, the specific error report is as follows, I hope you can give me a little bit of advice, thank you very much! I created a three-layer model in which there are two excavated lanes in the middle layer, and in addition, a trap column runs through all three layers of the medium.

!/usr/bin/env python

import os

########################################################################### ###########################################################################

SEMoutput='collapse_1201_ini_MESH' CUBIToutput='collapse_1201_ini_GEOCUBIT'

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

import cubit try: cubit.init([""]) except: pass

cubit.cmd('reset') cubit.cmd('brick x 100 y 75 z 25') cubit.cmd('move vol 1 x 50 y 37.5 z -12.5')

cubit.cmd('brick x 100 y 5 z 5') cubit.cmd('move vol 2 x 50 y 7.5 z -12.5')

cubit.cmd('brick x 100 y 5 z 5') cubit.cmd('move vol 3 x 50 y 67.5 z -12.5')

cubit.cmd('frustum z 21 radius 8 top 8') cubit.cmd('move vol 4 x 40 y 37.5 z -12.5')

SUBTRACT VOL B FROM VOL A

cubit.cmd('subtract vol 2 from vol 1') cubit.cmd('subtract vol 3 from vol 1') cubit.cmd('subtract vol 4 from vol 1')

cubit.cmd('frustum z 21 radius 8 top 8') cubit.cmd('move vol 6 x 40 y 37.5 z -12.5')

cubit.cmd('webcut vol all plane zplane offset -10 ') cubit.cmd('webcut vol all plane zplane offset -15 ')

cubit.cmd('delete volume 9 11')

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

Meshing the volumes

cubit.cmd('volume all size 1')

cubit.cmd('curve 48 49 61 interval 10') cubit.cmd('curve 48 49 61 scheme equal')

cubit.cmd('mesh volume all')

cubit.cmd('draw volume all')

cubit.cmd('save as "' + CUBIToutput + '/meshing.cub" overwrite')

End of meshing

obsolete:

import boundary_definition

import cubit2specfem3d

new:

from geocubitlib import boundary_definition

from geocubitlib import cubit2specfem3d

This is boundary_definition.py of GEOCUBIT

..... which extracts the bounding faces and defines them into blocks

boundary_definition.entities=['face']

boundary_definition.define_bc(boundary_definition.entities,parallel=True)

cubit.cmd('block 1001 face in surface 1 ') cubit.cmd('block 1001 name "face_topo"')

cubit.cmd('block 1002 face in surface 2 ') cubit.cmd('block 1002 name "face_bottom"')

cubit.cmd('block 1003 face in surface 42 58 67 62 73') cubit.cmd('block 1003 name "face_abs_xmin"')

cubit.cmd('block 1004 face in surface 43 63 72') cubit.cmd('block 1004 name "face_abs_ymin"')

cubit.cmd('block 1005 face in surface 44 64 59 66 70') cubit.cmd('block 1005 name "face_abs_xmax"')

cubit.cmd('block 1006 face in surface 41 60 69') cubit.cmd('block 1006 name "face_abs_ymax"')

boundary

boundary_definition.define_bc(parallel=True)

Define material properties for the 3 volumes

cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')

cubit.cmd('block 1 hex in volume 5 7') cubit.cmd('block 1 name "elastic 1" ') # material region cubit.cmd('block 1 attribute count 7') cubit.cmd('block 1 attribute index 1 1 ') # volume 2 cubit.cmd('block 1 attribute index 2 2700 ') # vp cubit.cmd('block 1 attribute index 3 1500 ')# vs cubit.cmd('block 1 attribute index 4 2200 ') # rho cubit.cmd('block 1 attribute index 5 9999.0 ')# Qkappa 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 9 10 11') cubit.cmd('block 2 name "elastic 2" ') # material region cubit.cmd('block 2 attribute count 7') cubit.cmd('block 2 attribute index 1 2 ') # volume 2 cubit.cmd('block 2 attribute index 2 2000 ') # vp cubit.cmd('block 2 attribute index 3 1150 ')# vs cubit.cmd('block 2 attribute index 4 1400 ') # rho cubit.cmd('block 2 attribute index 5 9999.0 ')# Qkappa 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 6 8 12') cubit.cmd('block 3 name "elastic 3" ') # material region cubit.cmd('block 3 attribute count 7') cubit.cmd('block 3 attribute index 1 3 ') # volume 2 cubit.cmd('block 3 attribute index 2 2500 ') # vp cubit.cmd('block 3 attribute index 3 1400 ')# vs cubit.cmd('block 3 attribute index 4 2150 ') # rho cubit.cmd('block 3 attribute index 5 9999.0 ')# Qkappa cubit.cmd('block 3 attribute index 6 9999.0 ')# Qmu cubit.cmd('block 3 attribute index 7 0 ') # anisotropy_flag

file export

cubit2specfem3d.export2SPECFEM3D(SEMoutput,cpml=True, cpml_size=2, top_absorbing=True ) image Also, I created a similar model, except without the trap columns, but it works fine. The specific error reported in the runtime here is as follows: error.txt

And, I made sure that the parameters used here regarding stability are correct. If you have any further questions, please let me know.

zuoxuan-de commented 1 day ago

And, this model with collapse works fine when I use Stacey boundary conditions.