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
390 stars 223 forks source link

models may not match #1609

Open liaohaiyang1534 opened 1 year ago

liaohaiyang1534 commented 1 year ago

Hi,

i use a python scipt to make a model (given at the end), and it looks like this in the software CoreForm Cubit 2022.11:

)EP3~DHA8G8%TPP{AN03@{6

when I set

# save AVS or OpenDX mesh files to check the mesh
SAVE_MESH_FILES                 = .true.

and check part_array.vtk from OUTPUT_FILES/DATABASES_MPI, it looks like this:

image

and I think this is alright, it matches the size of my model.

but when I check proc000010_vs.vtk:

image

It looks a little strange.

Should I ignore it? Or is there some potential bug?

the python script that I use:


> `
> from __future__ import print_function
> 
> import os
> import sys
> 
> elementsize = 4
> 
> # default directories
> SEMoutput='MESH'
> CUBIToutput='MESH_GEOCUBIT'
> 
> os.system('mkdir -p '+ SEMoutput)
> os.system('mkdir -p '+ CUBIToutput)
> 
> import cubit
> try:
>     #cubit.init([""])
>     cubit.init(["-noecho","-nojournal"])
> except:
>     pass
> 
> version = cubit.get_version()
> version_major = int(version.split(".")[0])
> version_minor = int(version.split(".")[1])
> print("cubit version: ",version)
> 
> cubit.cmd('reset')
> cubit.cmd('brick x 300 y 150 z 100')
> 
> # This seems to conflict with boundary_definition.py
> 
> cubit.cmd('volume 1 move x 150 y 75 z -50')
> 
> # create vertices for discontinuity
> cubit.cmd('split curve 9  distance 20')
> cubit.cmd('split curve 10  distance 20')
> cubit.cmd('split curve 11  distance 20')
> cubit.cmd('split curve 12  distance 20')
> 
> # create surface for interface
> cubit.cmd('create surface vertex 9 10 12 11')
> 
> cubit.cmd('section volume 1 with surface 7 keep normal')
> cubit.cmd('section volume 1 with surface 7 reverse')
> 
> # create vertices for auxiliary interface to allow for refinement
> cubit.cmd('split curve 29  distance 30')
> cubit.cmd('split curve 31  distance 30')
> cubit.cmd('split curve 32  distance 30')
> cubit.cmd('split curve 36  distance 30')
> 
> # create surface for buffer interface to refine BELOW the discontinuity
> cubit.cmd('create surface vertex 25 26 28 27')
> 
> cubit.cmd('section volume 3 with surface 19 keep normal')
> cubit.cmd('section volume 3 with surface 19 reverse')
> 
> cubit.cmd('delete volume 2 4')
> 
> cubit.cmd('merge all')
> cubit.cmd('imprint all')
> 
> # Meshing the volumes
> cubit.cmd('volume 3 size '+str(elementsize))
> cubit.cmd('mesh volume 3')
> 
> cubit.cmd('volume 1 size '+str(elementsize))
> cubit.cmd('mesh volume 1')
> 
> cubit.cmd('volume 5 size '+str(elementsize))
> cubit.cmd('mesh volume 5')
> 
> #### 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)
> 
> #### Define material properties for the 3 volumes ################
> cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
> 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  ')     # material 1
> cubit.cmd('block 1 attribute index 2 1200 ')  # vp
> cubit.cmd('block 1 attribute index 3 400 ')      # vs
> cubit.cmd('block 1 attribute index 4 2000 ')  # rho 
> cubit.cmd('block 2 attribute index 5 9999.0')      # Q_kappa
> cubit.cmd('block 2 attribute index 6 9999.0')      # Q_mu
> cubit.cmd('block 2 attribute index 7 0 ')     # anisotropy_flag
> 
> 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  ')     # material 2
> cubit.cmd('block 2 attribute index 2 2000 ')  # vp
> cubit.cmd('block 2 attribute index 3 700 ')  # vs
> cubit.cmd('block 2 attribute index 4 3200 ')  # rho
> cubit.cmd('block 2 attribute index 5 9999.0')      # Q_kappa
> cubit.cmd('block 2 attribute index 6 9999.0')      # Q_mu
> cubit.cmd('block 2 attribute index 7 0 ')     # anisotropy_flag
> 
> cubit.cmd('block 3 name "elastic 3" ')        # elastic material region
> cubit.cmd('block 3 attribute count 7')
> cubit.cmd('block 3 attribute index 1 3  ')     
> cubit.cmd('block 3 attribute index 2 2500 ')  # vp
> cubit.cmd('block 3 attribute index 3 800 ')  # vs
> cubit.cmd('block 3 attribute index 4 3600 ')  # rho
> cubit.cmd('block 3 attribute index 5 9999.0')      # Q_kappa
> cubit.cmd('block 3 attribute index 6 9999.0')      # Q_mu
> cubit.cmd('block 3 attribute index 7 0 ')     # anisotropy_flag
> 
> cubit.cmd('export mesh "' + CUBIToutput + '/top.e" dimension 3 overwrite')
> cubit.cmd('save as "' + CUBIToutput + '/top.cub" overwrite')
> 
> #### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT
> 
> cubit2specfem3d.export2SPECFEM3D(SEMoutput)
> 
> # screen shot
> # (could crash version < 16.4)
> if version_major >= 16 and version_minor >= 4:
>     cubit.cmd('view top')
>     cubit.cmd('hardcopy "' + CUBIToutput + '/waterlayered.png" png')
> 
> # all files needed by SCOTCH are now in directory MESH`
> 

I'm very appreciated to have any advice and suggestions for this. (If you need any more clarifications, please let me know)

Regards

Haiyang Liao haiyangliao@smail.nju.edu.cn

homnath commented 1 year ago

Hi Haiyang Liao,

How does it look like when you visualize all processor files, proc000000_vs.vtk, proc000001_vs.vtk, proc000002_vs.vtk, etc., altogether? The individual mesh may look irregular due to partitioning but altogether, they should match your original model.

Best, Hom Nath


From: liaohaiyang1534 @.> Sent: Saturday, May 27, 2023 11:15 PM To: SPECFEM/specfem3d @.> Cc: Subscribed @.***> Subject: [SPECFEM/specfem3d] models may not match (Issue #1609)

Hi,

i use a python scipt to make a model (given at the end), and it looks like this in the software CoreForm Cubit 2022.11:

[)EP3~DHA8G8%TPP{AN03@{6]https://user-images.githubusercontent.com/104192749/241474009-81f180e1-06e3-4825-8374-0acb5c092df9.png

when I set

save AVS or OpenDX mesh files to check the mesh

SAVE_MESH_FILES = .true.

and check part_array.vtk from OUTPUT_FILES/DATABASES_MPI, it looks like this:

[image]https://user-images.githubusercontent.com/104192749/241476203-a3c7f16f-a220-43b3-8e37-65cd37947672.png

and I think this is alright, it matches the size of my model.

but when I check proc000010_vs.vtk:

[image]https://user-images.githubusercontent.com/104192749/241476548-8c2f60ac-5cca-47b2-9b18-f12a0b4afcfc.png

It looks a little strange.

Should I ignore it? Or is there some potential bug?

the python script that I use:

` from future import print_function

import os import sys

elementsize = 4

default directories

SEMoutput='MESH' CUBIToutput='MESH_GEOCUBIT'

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

import cubit try:

cubit.init([""])

cubit.init(["-noecho","-nojournal"])

except: pass

version = cubit.get_version() version_major = int(version.split(".")[0]) version_minor = int(version.split(".")[1]) print("cubit version: ",version)

cubit.cmd('reset') cubit.cmd('brick x 300 y 150 z 100')

This seems to conflict with boundary_definition.py

cubit.cmd('volume 1 move x 150 y 75 z -50')

create vertices for discontinuity

cubit.cmd('split curve 9 distance 20') cubit.cmd('split curve 10 distance 20') cubit.cmd('split curve 11 distance 20') cubit.cmd('split curve 12 distance 20')

create surface for interface

cubit.cmd('create surface vertex 9 10 12 11')

cubit.cmd('section volume 1 with surface 7 keep normal') cubit.cmd('section volume 1 with surface 7 reverse')

create vertices for auxiliary interface to allow for refinement

cubit.cmd('split curve 29 distance 30') cubit.cmd('split curve 31 distance 30') cubit.cmd('split curve 32 distance 30') cubit.cmd('split curve 36 distance 30')

create surface for buffer interface to refine BELOW the discontinuity

cubit.cmd('create surface vertex 25 26 28 27')

cubit.cmd('section volume 3 with surface 19 keep normal') cubit.cmd('section volume 3 with surface 19 reverse')

cubit.cmd('delete volume 2 4')

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

Meshing the volumes

cubit.cmd('volume 3 size '+str(elementsize)) cubit.cmd('mesh volume 3')

cubit.cmd('volume 1 size '+str(elementsize)) cubit.cmd('mesh volume 1')

cubit.cmd('volume 5 size '+str(elementsize)) cubit.cmd('mesh volume 5')

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)

Define material properties for the 3 volumes

cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################') 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 ') # material 1 cubit.cmd('block 1 attribute index 2 1200 ') # vp cubit.cmd('block 1 attribute index 3 400 ') # vs cubit.cmd('block 1 attribute index 4 2000 ') # rho cubit.cmd('block 2 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 2 attribute index 6 9999.0') # Q_mu cubit.cmd('block 2 attribute index 7 0 ') # anisotropy_flag

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 ') # material 2 cubit.cmd('block 2 attribute index 2 2000 ') # vp cubit.cmd('block 2 attribute index 3 700 ') # vs cubit.cmd('block 2 attribute index 4 3200 ') # rho cubit.cmd('block 2 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 2 attribute index 6 9999.0') # Q_mu cubit.cmd('block 2 attribute index 7 0 ') # anisotropy_flag

cubit.cmd('block 3 name "elastic 3" ') # elastic material region cubit.cmd('block 3 attribute count 7') cubit.cmd('block 3 attribute index 1 3 ') cubit.cmd('block 3 attribute index 2 2500 ') # vp cubit.cmd('block 3 attribute index 3 800 ') # vs cubit.cmd('block 3 attribute index 4 3600 ') # rho cubit.cmd('block 3 attribute index 5 9999.0') # Q_kappa cubit.cmd('block 3 attribute index 6 9999.0') # Q_mu cubit.cmd('block 3 attribute index 7 0 ') # anisotropy_flag

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

Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT

cubit2specfem3d.export2SPECFEM3D(SEMoutput)

screen shot

(could crash version < 16.4)

if version_major >= 16 and version_minor >= 4: cubit.cmd('view top') cubit.cmd('hardcopy "' + CUBIToutput + '/waterlayered.png" png')

all files needed by SCOTCH are now in directory MESH`

I'm very appreciated to have any advice and suggestions for this. (If you need any more clarifications, please let me know)

Regards

Haiyang Liao @.**@.>

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