diff --git a/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py b/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py index 263617f76..5ea4b34cf 100755 --- a/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py +++ b/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py @@ -652,14 +652,16 @@ def block_definition(self): self.bc = bc print('HEX Blocks:') for m, f in zip(self.block_mat, self.block_flag): - print('block ', m, 'material flag ', f) + print(' block ', m, 'material flag ', f) print('Absorbing Boundary Conditions:') for m, f in zip(self.block_bc, self.block_bc_flag): - print('bc ', m, 'bc flag ', f) - print('Topography (free surface)') - print(self.topography) - print('Free surface') - print(self.free) + print(' bc ', m, 'bc flag ', f) + if self.topography: + print('Topography (free surface):') + print(' ',self.topography) + if self.free: + print('Free surface:') + print(' ',self.free) except: print('****************************************') print('sorry, no blocks or blocks not properly defined') @@ -765,7 +767,7 @@ def mat_parameter(self, properties): else: raise RuntimeError('Error: material id must be strictly positive or negative, but not equal to 0') # info output - print("material: ",txt) + print(" material: ",txt) return txt def nummaterial_write(self, nummaterial_name, placeholder=True): @@ -853,8 +855,7 @@ def create_hexnode_string(self, hexa, hexnode_string=True): else: map(int, txt.split()) - def create_facenode_string(self, hexa, face, normal=None, cknormal=True, - facenode_string=True): + def create_facenode_string(self, hexa, face, normal=None, cknormal=True, facenode_string=True): nodes = self.get_face_connectivity(face) if cknormal: nodes_ok = self.normal_check(nodes[0:4], normal) @@ -886,7 +887,7 @@ def mesh_write(self, mesh_name): meshfile.write(str(num_elems) + '\n') for block, flag in zip(self.block_mat, self.block_flag): hexes = cubit.get_block_hexes(block) - print('block ', block, ' hexes ', len(hexes)) + print(' block ', block, ' hexes ', len(hexes)) for hexa in hexes: txt = self.create_hexnode_string(hexa) meshfile.write(txt) @@ -897,7 +898,7 @@ def material_write(self, mat_name): mat = open(mat_name, 'w') print('Writing ' + mat_name + '.....') for block, flag in zip(self.block_mat, self.block_flag): - print('block ', block, 'flag ', flag) + print(' block ', block, 'flag ', flag) hexes = cubit.get_block_hexes(block) for hexa in hexes: mat.write(('%10i %10i\n') % (hexa, flag)) @@ -949,8 +950,7 @@ def free_write(self, freename=None): for block, flag in zip(self.block_bc, self.block_bc_flag): if block == self.topography: name = cubit.get_exodus_entity_name('block', block) - print('free surface (topography) block name:', \ - name, 'id:', block) + print(' free surface (topography): block name:', name, 'id:', block) quads_all = cubit.get_block_faces(block) print(' number of faces = ', len(quads_all)) dic_quads_all = dict(zip(quads_all, quads_all)) @@ -968,13 +968,13 @@ def free_write(self, freename=None): faces = cubit.get_sub_elements('hex', h, 2) for f in faces: if f in s: - txt = self.create_facenode_string( - h, f, normal, cknormal=True) + # topography surface: checks face orientation against vertical normal (0,0,1) + txt = self.create_facenode_string(h, f, normal, cknormal=True) freehex.write(txt) print(' 100 %') elif block == self.free: name = cubit.get_exodus_entity_name('block', block) - print('free surface block name:', name, 'id:', block) + print(' free surface: block name:', name, 'id:', block) quads_all = cubit.get_block_faces(block) print(' number of faces = ', len(quads_all)) dic_quads_all = dict(zip(quads_all, quads_all)) @@ -992,8 +992,8 @@ def free_write(self, freename=None): faces = cubit.get_sub_elements('hex', h, 2) for f in faces: if f in s: - txt = self.create_facenode_string( - h, f, normal, cknormal=False) + # free surface: no face orientation check against vertical normal + txt = self.create_facenode_string(h, f, normal, cknormal=False) freehex.write(txt) print(' 100 %') freehex.close() @@ -1122,60 +1122,60 @@ def abs_write(self, absname=None): # loops through all block definitions list_hex = cubit.parse_cubit_list('hex', 'all') for block, flag in zip(self.block_bc, self.block_bc_flag): - if block != self.topography: + if block != self.topography and block != self.free: name = cubit.get_exodus_entity_name('block', block) print(' block name:', name, 'id:', block) cknormal = True abshex_local = False # opens file if re.search('xmin', name): - print("xmin") + print(" xmin") abshex_local = open(absname + '_xmin', 'w') normal = (-1, 0, 0) elif re.search('xmax', name): - print("xmax") + print(" xmax") abshex_local = open(absname + '_xmax', 'w') normal = (1, 0, 0) elif re.search('ymin', name): - print("ymin") + print(" ymin") abshex_local = open(absname + '_ymin', 'w') normal = (0, -1, 0) elif re.search('ymax', name): - print("ymax") + print(" ymax") abshex_local = open(absname + '_ymax', 'w') normal = (0, 1, 0) elif re.search('bottom', name): - print("bottom") + print(" bottom") abshex_local = open(absname + '_bottom', 'w') normal = (0, 0, -1) elif re.search('abs', name): - print("abs all - experimental, check the output") + print(" abs all - experimental, check the output") cknormal = False abshex_local = open(absname, 'w') normal = None else: if block == 1003: - print("xmin") + print(" xmin") abshex_local = open(absname + '_xmin', 'w') normal = (-1, 0, 0) elif block == 1004: - print("ymin") + print(" ymin") abshex_local = open(absname + '_ymin', 'w') normal = (0, -1, 0) elif block == 1005: - print("xmax") + print(" xmax") abshex_local = open(absname + '_xmax', 'w') normal = (1, 0, 0) elif block == 1006: - print("ymax") + print(" ymax") abshex_local = open(absname + '_ymax', 'w') normal = (0, 1, 0) elif block == 1002: - print("bottom") + print(" bottom") abshex_local = open(absname + '_bottom', 'w') normal = (0, 0, -1) elif block == 1000: - print("custom") + print(" custom") abshex_local = open(absname, 'w') cknormal = False normal = None @@ -1199,8 +1199,7 @@ def abs_write(self, absname=None): faces = cubit.get_sub_elements('hex', h, 2) for f in faces: if f in s: - txt = self.create_facenode_string( - h, f, normal=normal, cknormal=cknormal) + txt = self.create_facenode_string(h, f, normal=normal, cknormal=cknormal) abshex_local.write(txt) abshex_local.close() print(' 100 %') @@ -1215,19 +1214,23 @@ def surface_write(self, pathdir=None): # > block 10 name 'moho_surface' import re for block in self.block_bc: - if block != self.topography: + if block != self.topography and block != self.free: name = cubit.get_exodus_entity_name('block', block) - # skips block names like face_abs**, face_topo** + # skips block names like face_abs**, face_topo**, face_free** if re.search('abs', name): continue elif re.search('topo', name): continue + elif re.search('free', name): + continue elif re.search('surface', name): + # e.g. moho_surface filename = pathdir + name + '_file' else: continue # gets face elements - print(' surface block name: ', name, 'id: ', block) + print('Optional surface:') + print(' surface: block name: ', name, 'id: ', block) quads_all = cubit.get_block_faces(block) print(' face = ', len(quads_all)) if len(quads_all) == 0: @@ -1243,8 +1246,8 @@ def surface_write(self, pathdir=None): Nhex = len(list_hex) # Create a set to speed up "if f in dic_quads_all.keys()" big time s = set(dic_quads_all.keys()) - percentageDone = -1 - for h in list_hex: + lastDisplayed = -1 + for ih, h in enumerate(list_hex): percentageDone = float(ih)/Nhex*100.0 if int(percentageDone)%10 == 0 and lastDisplayed != int(percentageDone): print(' ', int(percentageDone), '%') @@ -1252,8 +1255,7 @@ def surface_write(self, pathdir=None): faces = cubit.get_sub_elements('hex', h, 2) for f in faces: if f in s: - txt = self.create_facenode_string( - h, f, cknormal=False) + txt = self.create_facenode_string(h, f, cknormal=False) surfhex_local.write(txt) # closes file surfhex_local.close() @@ -1388,7 +1390,7 @@ def _get_bc_flag(self, block): normal = (0, 0, -1) elif re.search('abs', name): label = "all" - print("abs all - experimental, check the output") + print(" abs all - experimental, check the output") cknormal = False else: if block == 1003: @@ -1424,9 +1426,7 @@ def _get_bc_flag(self, block): # faces = cubit.get_sub_elements('hex', h, 2) # for f in faces: # if dic_quads_all.has_key(f): -# abs_array.append(self.create_facenode_string( -# h, f, normal=normal, cknormal=cknormal, -# facenode_string=False)) +# abs_array.append(self.create_facenode_string(h, f, normal=normal, cknormal=cknormal, facenode_string=False)) # print('Ok') # cubit.cd('set info on') # cubit.cmd('set echo on') @@ -1474,9 +1474,7 @@ def _get_bc_flag(self, block): # for f in faces: # if dic_quads_all.has_key(f): # # print(f) -# free_array.append(self.create_facenode_string( -# h, f, normal, cknormal=True, -# facenode_string=False)) +# free_array.append(self.create_facenode_string(h, f, normal, cknormal=True, facenode_string=False)) # elif block == self.free: # name = cubit.get_exodus_entity_name('block', block) # print('free surface block name:', name, 'id:', block) @@ -1488,9 +1486,7 @@ def _get_bc_flag(self, block): # faces = cubit.get_sub_elements('hex', h, 2) # for f in faces: # if dic_quads_all.has_key(f): -# free_array.append(self.create_facenode_string( -# h, f, normal, cknormal=False, -# facenode_string=False)) +# free_array.append(self.create_facenode_string(h, f, normal, cknormal=False, facenode_string=False)) # print('Ok') # cubit.cmd('set info on') # cubit.cmd('set echo on') diff --git a/EXAMPLES/applications/Mount_StHelens/mesh_mount.py b/EXAMPLES/applications/Mount_StHelens/mesh_mount.py index 593ed25a7..26291371d 100755 --- a/EXAMPLES/applications/Mount_StHelens/mesh_mount.py +++ b/EXAMPLES/applications/Mount_StHelens/mesh_mount.py @@ -158,8 +158,8 @@ else: print("#" + time.strftime("%a, %d %b %Y %H:%M:%S", time.localtime())) print("# imprinting volume, this will take around 1 min, please be patience...") - cubit.cmd('merge all') cubit.cmd('imprint all') + cubit.cmd('merge all') # exports only surfaces which will create single volume cubit.cmd('export acis "topo_2.acis" surface 3 10 12 14 15 9 ascii overwrite') diff --git a/EXAMPLES/applications/Mount_StHelens/mesh_mount_stl.py b/EXAMPLES/applications/Mount_StHelens/mesh_mount_stl.py index 12713490b..825a020f4 100755 --- a/EXAMPLES/applications/Mount_StHelens/mesh_mount_stl.py +++ b/EXAMPLES/applications/Mount_StHelens/mesh_mount_stl.py @@ -132,7 +132,7 @@ # imprints topography surfaces into brick volume, creates connected surfaces # note: this is a develop feature cubit.cmd('imprint all') - +cubit.cmd('merge all') cubit.cmd('version') diff --git a/EXAMPLES/applications/fault_examples/tpv5/TPV5.py b/EXAMPLES/applications/fault_examples/tpv5/TPV5.py index ffceab079..7747d3472 100755 --- a/EXAMPLES/applications/fault_examples/tpv5/TPV5.py +++ b/EXAMPLES/applications/fault_examples/tpv5/TPV5.py @@ -131,6 +131,7 @@ cubit.cmd("imprint all") cubit.cmd("merge all") + cubit.cmd("surface 1 size "+str(elementsize)) cubit.cmd("volume 1 size "+str(elementsize)) cubit.cmd("surface 1 scheme pave") diff --git a/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py b/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py index ff549e299..9bb5c0d7f 100755 --- a/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py +++ b/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py @@ -94,8 +94,8 @@ cubit.cmd('delete volume 2 4') -cubit.cmd('merge all') cubit.cmd('imprint all') +cubit.cmd('merge all') cubit.cmd('volume 3 size '+str(elementsize)) cubit.cmd('mesh volume 3') diff --git a/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8.py b/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8.py index 044d7fdb6..5cbc766b2 100755 --- a/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8.py +++ b/EXAMPLES/applications/layered_halfspace/2lay_mesh_boundary_fig8.py @@ -94,8 +94,8 @@ cubit.cmd('delete volume 2 4') -cubit.cmd('merge all') cubit.cmd('imprint all') +cubit.cmd('merge all') # Meshing the volumes ## middle volume diff --git a/EXAMPLES/applications/waterlayered_halfspace/make_mesh_waterlayer-default.py b/EXAMPLES/applications/waterlayered_halfspace/make_mesh_waterlayer-default.py index 84095d03d..8571d2322 100755 --- a/EXAMPLES/applications/waterlayered_halfspace/make_mesh_waterlayer-default.py +++ b/EXAMPLES/applications/waterlayered_halfspace/make_mesh_waterlayer-default.py @@ -88,8 +88,8 @@ cubit.cmd('delete volume 2 4') -cubit.cmd('merge all') cubit.cmd('imprint all') +cubit.cmd('merge all') # Meshing the volumes cubit.cmd('volume 3 size '+str(elementsize)) diff --git a/EXAMPLES/applications/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py b/EXAMPLES/applications/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py index fe312b2e1..c9b52e079 100755 --- a/EXAMPLES/applications/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py +++ b/EXAMPLES/applications/waterlayered_halfspace/waterlayer_mesh_boundary_fig8.py @@ -68,8 +68,8 @@ cubit.cmd('delete volume 2 4') -cubit.cmd('merge all') cubit.cmd('imprint all') +cubit.cmd('merge all') # Meshing the volumes cubit.cmd('volume 3 size 3589.2') diff --git a/EXAMPLES/applications/waterlayered_halfspace/waterlayer_only-nodoubling.py b/EXAMPLES/applications/waterlayered_halfspace/waterlayer_only-nodoubling.py index c5ed17ce9..950c1bcbd 100755 --- a/EXAMPLES/applications/waterlayered_halfspace/waterlayer_only-nodoubling.py +++ b/EXAMPLES/applications/waterlayered_halfspace/waterlayer_only-nodoubling.py @@ -64,8 +64,8 @@ cubit.cmd('delete volume 2 4') -cubit.cmd('merge all') cubit.cmd('imprint all') +cubit.cmd('merge all') # Meshing the volumes #elementsize = 1196.4 #hi-resolution diff --git a/EXAMPLES/applications/waterlayered_poroelastic/create_mesh.py b/EXAMPLES/applications/waterlayered_poroelastic/create_mesh.py index 767420f07..dea735f17 100755 --- a/EXAMPLES/applications/waterlayered_poroelastic/create_mesh.py +++ b/EXAMPLES/applications/waterlayered_poroelastic/create_mesh.py @@ -95,8 +95,8 @@ cubit.cmd('delete volume 4') # merges surfaces -cubit.cmd('merge all') cubit.cmd('imprint all') +cubit.cmd('merge all') # resets numbering cubit.cmd('compress all') diff --git a/src/decompose_mesh/part_decompose_mesh.F90 b/src/decompose_mesh/part_decompose_mesh.F90 index e1468dffe..f85afb035 100644 --- a/src/decompose_mesh/part_decompose_mesh.F90 +++ b/src/decompose_mesh/part_decompose_mesh.F90 @@ -620,7 +620,7 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, & implicit none integer, intent(in) :: iproc integer, intent(in) :: nspec - integer, intent(in) :: nb_wd + integer, intent(in) :: nb_wd integer, dimension(:), pointer :: glob2loc_elmnts integer, dimension(1:nspec) :: part integer, dimension(1:nb_wd), intent(in) :: boundary_to_ispec_wd, side_wd @@ -632,7 +632,7 @@ subroutine write_wavefield_discontinuity_database(iproc, outputpath_name, & character(len=MAX_STRING_LEN) :: prname write(prname, "('proc',i6.6,'_')") iproc open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=trim(outputpath_name)//'/'//trim(prname)//& + file = trim(outputpath_name)//'/'//trim(prname)//& trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & form='unformatted', action='write') local_nb_wd = 0 diff --git a/src/decompose_mesh/read_mesh_files.F90 b/src/decompose_mesh/read_mesh_files.F90 index 389be55e2..706ecf290 100644 --- a/src/decompose_mesh/read_mesh_files.F90 +++ b/src/decompose_mesh/read_mesh_files.F90 @@ -39,7 +39,7 @@ subroutine read_mesh_files() use fault_scotch, only: ANY_FAULT,read_fault_files,save_nodes_coords,close_faults use constants, only: FNAME_WAVEFIELD_DISCONTINUITY_INTERFACE - + implicit none ! local parameters @@ -592,7 +592,7 @@ subroutine read_mesh_files() read(IIN_DB, *, iostat=ier) boundary_to_ispec_wd(ib), side_wd(ib) enddo close(IIN_DB) - endif + endif ! reads in absorbing boundary files open(unit=IIN_DB, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmin', & diff --git a/src/decompose_mesh/write_mesh_databases.F90 b/src/decompose_mesh/write_mesh_databases.F90 index 7c42244a2..3b3259bee 100644 --- a/src/decompose_mesh/write_mesh_databases.F90 +++ b/src/decompose_mesh/write_mesh_databases.F90 @@ -201,7 +201,7 @@ subroutine write_mesh_databases() enddo endif - + !! setting up wavefield discontinuity boundary if (IS_WAVEFIELD_DISCONTINUITY) then do ipart = 0, nparts-1 diff --git a/src/generate_databases/save_arrays_solver.F90 b/src/generate_databases/save_arrays_solver.F90 index 59cde021b..5fee37dc5 100644 --- a/src/generate_databases/save_arrays_solver.F90 +++ b/src/generate_databases/save_arrays_solver.F90 @@ -110,6 +110,7 @@ subroutine save_arrays_solver_mesh() if (myrank == 0) then write(IMAIN,*) ' using binary file format' write(IMAIN,*) ' database file (for rank 0): ',trim(filename) + write(IMAIN,*) call flush_IMAIN() endif @@ -585,7 +586,7 @@ subroutine save_arrays_solver_files() ! user output call synchronize_all() if (myrank == 0) then - write(IMAIN,*) ' saving additonal mesh files with surface/coupling points' + write(IMAIN,*) ' saving additional mesh files with surface/coupling points' call flush_IMAIN() endif diff --git a/src/generate_databases/wavefield_discontinuity_generate_databases.f90 b/src/generate_databases/wavefield_discontinuity_generate_databases.f90 index 12c3292d2..e531bbd32 100644 --- a/src/generate_databases/wavefield_discontinuity_generate_databases.f90 +++ b/src/generate_databases/wavefield_discontinuity_generate_databases.f90 @@ -1,3 +1,31 @@ +!===================================================================== +! +! S p e c f e m 3 D +! ----------------- +! +! Main historical authors: Dimitri Komatitsch and Jeroen Tromp +! CNRS, France +! and Princeton University, USA +! (there are currently many more authors!) +! (c) October 2017 +! +! This program is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License along +! with this program; if not, write to the Free Software Foundation, Inc., +! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +! +!===================================================================== + + module wavefield_discontinuity_generate_databases use constants, only: CUSTOM_REAL @@ -23,7 +51,7 @@ module wavefield_discontinuity_generate_databases !! written in solver database and used in solver integer, dimension(:), allocatable :: ispec_to_elem_wd - !! number of distinct gll points on the boundary + !! number of distinct GLL points on the boundary !! written in solver database and used in solver integer :: nglob_wd @@ -75,7 +103,7 @@ subroutine read_partition_files_wavefield_discontinuity() implicit none !integer :: ier open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & + file = trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & action='read', form='unformatted') read(IFILE_WAVEFIELD_DISCONTINUITY) nb_wd allocate(boundary_to_ispec_wd(nb_wd), side_wd(nb_wd)) @@ -90,7 +118,7 @@ subroutine save_arrays_solver_mesh_wavefield_discontinuity() use generate_databases_par, only: prname implicit none open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_DATABASE), & + file = trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_DATABASE), & action='write', form='unformatted') write(IFILE_WAVEFIELD_DISCONTINUITY) ispec_to_elem_wd write(IFILE_WAVEFIELD_DISCONTINUITY) nglob_wd @@ -112,7 +140,7 @@ end subroutine save_arrays_solver_mesh_wavefield_discontinuity subroutine setup_boundary_wavefield_discontinuity() use generate_databases_par, only: NDIM, NGLLX, NGLLY, NGLLZ, CUSTOM_REAL, & NGLLSQUARE, ibool, NSPEC_AB - + implicit none integer :: i1, i2, j1, j2, k1, k2, i, j, k integer :: ib, iside, ispec, iglob, ispec_wd, iglob_wd @@ -126,7 +154,7 @@ subroutine setup_boundary_wavefield_discontinuity() nfaces_wd = 0 ispec_to_elem_wd(:) = 0 elem_list_temp(:) = 0 - gllp_list_temp(:) = 0 + gllp_list_temp(:) = 0 do ib = 1, nb_wd ispec = boundary_to_ispec_wd(ib) iside = side_wd(ib) @@ -137,7 +165,7 @@ subroutine setup_boundary_wavefield_discontinuity() elem_list_temp(nspec_wd+1) = ispec nspec_wd = nspec_wd + 1 endif - do i=i1,i2; do j=j1,j2; do k=k1,k2 + do i = i1,i2; do j = j1,j2; do k = k1,k2 iglob = ibool(i,j,k,ispec) call find_point_wd(iglob, gllp_list_temp, nglob_wd, iglob_wd) if (iglob_wd == 0) then @@ -162,7 +190,7 @@ subroutine setup_boundary_wavefield_discontinuity() iside = side_wd(ib) call get_points_boundary_wd(iside, i1, i2, j1, j2, k1, k2, is_face) call find_point_wd(ispec, elem_list_temp, nspec_wd, ispec_wd) - do i=i1,i2; do j=j1,j2; do k=k1,k2 + do i = i1,i2; do j = j1,j2; do k = k1,k2 iglob = ibool(i,j,k,ispec) call find_point_wd(iglob, gllp_list_temp, nglob_wd, iglob_wd) ibool_wd(i,j,k,ispec_wd) = iglob_wd @@ -181,7 +209,7 @@ subroutine setup_boundary_wavefield_discontinuity() do ispec = 1, NSPEC_AB ispec_wd = ispec_to_elem_wd(ispec) if (ispec_wd > 0) then - do k=1,NGLLZ; do j=1,NGLLY; do i=1, NGLLX + do k = 1,NGLLZ; do j = 1,NGLLY; do i = 1, NGLLX iglob_wd= ibool_wd(i,j,k,ispec_wd) if (iglob_wd > 0) then call get_mass_wd(i, j, k, ispec, val_mass) @@ -194,7 +222,7 @@ subroutine setup_boundary_wavefield_discontinuity() end subroutine setup_boundary_wavefield_discontinuity subroutine write_discontinuity_surface_file() - use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,NGLLSQUARE,CUSTOM_REAL,& + use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,NGLLSQUARE,CUSTOM_REAL, & IFILE_WAVEFIELD_DISCONTINUITY use generate_databases_par, only: prname, myrank, LOCAL_PATH, & nodes_coords_ext_mesh, elmnts_ext_mesh, NGNOD @@ -215,7 +243,7 @@ subroutine write_discontinuity_surface_file() double precision, dimension(NGNOD) :: xelm,yelm,zelm call create_name_database(prname,myrank,LOCAL_PATH) open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=prname(1:len_trim(prname))//'wavefield_discontinuity_points',& + file = prname(1:len_trim(prname))//'wavefield_discontinuity_points', & action='write', form='formatted') do iglob_wd = 1, nglob_wd iglob = boundary_to_iglob_wd(iglob_wd) @@ -235,7 +263,7 @@ subroutine write_discontinuity_surface_file() call get_shape3D(shape3D_shrink, dershape3D_shrink, & xigll_shrink, yigll_shrink, zigll_shrink, NGNOD, NGLLX, NGLLY, NGLLZ) open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=prname(1:len_trim(prname))//'wavefield_discontinuity_faces',& + file = prname(1:len_trim(prname))//'wavefield_discontinuity_faces', & action='write', form='formatted') do iface_wd = 1, nfaces_wd ispec = face_ispec_wd(iface_wd) @@ -245,7 +273,7 @@ subroutine write_discontinuity_surface_file() yelm(ia) = nodes_coords_ext_mesh(2,iglob) zelm(ia) = nodes_coords_ext_mesh(3,iglob) enddo - call calc_coords(xstore_shrink(1,1,1), ystore_shrink(1,1,1),& + call calc_coords(xstore_shrink(1,1,1), ystore_shrink(1,1,1), & zstore_shrink(1,1,1), xelm,yelm,zelm,shape3D_shrink) do igll = 1, NGLLSQUARE i = face_ijk_wd(1,igll,iface_wd) diff --git a/src/meshfem3D/create_meshfem_mesh.f90 b/src/meshfem3D/create_meshfem_mesh.f90 index 5fe925c07..3b3c46e35 100644 --- a/src/meshfem3D/create_meshfem_mesh.f90 +++ b/src/meshfem3D/create_meshfem_mesh.f90 @@ -82,7 +82,7 @@ subroutine create_meshfem_mesh() ! HDF5 file i/o use shared_parameters, only: HDF5_ENABLED - + !! setting up wavefield discontinuity interface use shared_parameters, only: IS_WAVEFIELD_DISCONTINUITY diff --git a/src/meshfem3D/get_wavefield_discontinuity.f90 b/src/meshfem3D/get_wavefield_discontinuity.f90 index 7f3f26090..2c9e5c143 100644 --- a/src/meshfem3D/get_wavefield_discontinuity.f90 +++ b/src/meshfem3D/get_wavefield_discontinuity.f90 @@ -34,7 +34,7 @@ subroutine write_wavefield_discontinuity_database() implicit none ! local variables open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=prname(1:len_trim(prname))//& + file = prname(1:len_trim(prname))//& trim(FNAME_WAVEFIELD_DISCONTINUITY_MESH), & form='unformatted', action='write') write(IFILE_WAVEFIELD_DISCONTINUITY) nb_wd @@ -72,7 +72,6 @@ end subroutine write_wavefield_discontinuity_file subroutine find_wavefield_discontinuity_elements() ! read the wavefield_discontinuity_box file ! find the wavefield discontinuity interface - use constants, only: CUSTOM_REAL use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & FNAME_WAVEFIELD_DISCONTINUITY_BOX use meshfem_par, only: xstore,ystore,zstore,nspec @@ -85,7 +84,7 @@ subroutine find_wavefield_discontinuity_elements() IS_EXTRAPOLATION_MODE logical :: covered(26) double precision :: x_min, x_max, y_min, y_max, z_min, z_max - double precision :: dx, dy, dz, x_mid, y_mid, z_mid, ratio_small=1.0e-6 + double precision :: dx, dy, dz, x_mid, y_mid, z_mid, ratio_small = 1.0e-6 open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & file=trim(FNAME_WAVEFIELD_DISCONTINUITY_BOX), & form='formatted', action='read') diff --git a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 index 8d3627c0e..873cb01fb 100644 --- a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 +++ b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 @@ -40,7 +40,7 @@ subroutine compute_forces_viscoelastic_calling() use constants, only: FAULT_SYNCHRONIZE_DISPL_VELOC,FAULT_SYNCHRONIZE_ACCEL use fault_solver_dynamic, only: bc_dynflt_set3d_all,SIMULATION_TYPE_DYN,fault_output_synchronize_GPU,NT_RECORD_LENGTH use fault_solver_kinematic, only: bc_kinflt_set_all,SIMULATION_TYPE_KIN - + !! solving wavefield discontinuity problem with non-split-node scheme use wavefield_discontinuity_solver, only: & add_traction_discontinuity, read_wavefield_discontinuity_file diff --git a/src/specfem3D/wavefield_discontinuity_solver.f90 b/src/specfem3D/wavefield_discontinuity_solver.f90 index e13886ffa..2c4c94841 100644 --- a/src/specfem3D/wavefield_discontinuity_solver.f90 +++ b/src/specfem3D/wavefield_discontinuity_solver.f90 @@ -1,3 +1,30 @@ +!===================================================================== +! +! S p e c f e m 3 D +! ----------------- +! +! Main historical authors: Dimitri Komatitsch and Jeroen Tromp +! CNRS, France +! and Princeton University, USA +! (there are currently many more authors!) +! (c) October 2017 +! +! This program is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License along +! with this program; if not, write to the Free Software Foundation, Inc., +! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +! +!===================================================================== + !! Solving the wavefield discontinuity problem with a non-split-node !! scheme !! Tianshi Liu, 2023.5 @@ -9,7 +36,7 @@ module wavefield_discontinuity_solver !! read from solver database integer, dimension(:), allocatable :: ispec_to_elem_wd - !! number of distinct gll points on the boundary + !! number of distinct GLL points on the boundary !! read from solver database integer :: nglob_wd @@ -66,15 +93,15 @@ module wavefield_discontinuity_solver real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: traction_wd contains - + subroutine read_mesh_databases_wavefield_discontinuity() use constants, only: IFILE_WAVEFIELD_DISCONTINUITY, & FNAME_WAVEFIELD_DISCONTINUITY_DATABASE - use specfem_par, only: CUSTOM_REAL, prname, & + use specfem_par, only: prname, & NSPEC_AB, NGLLX, NGLLY, NGLLZ, NDIM, NGLLSQUARE implicit none open(unit=IFILE_WAVEFIELD_DISCONTINUITY, & - file=trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_DATABASE), & + file = trim(prname)//trim(FNAME_WAVEFIELD_DISCONTINUITY_DATABASE), & action='read', form='unformatted') allocate(ispec_to_elem_wd(NSPEC_AB)) @@ -137,9 +164,9 @@ subroutine add_displacement_discontinuity_element(ispec, dummyx_loc, & integer :: ispec_wd, i, j, k, iglob_wd ispec_wd = ispec_to_elem_wd(ispec) if (ispec_wd /= 0) then - do k=1,NGLLZ - do j=1,NGLLY - do i=1,NGLLX + do k = 1,NGLLZ + do j = 1,NGLLY + do i = 1,NGLLX iglob_wd = ibool_wd(i,j,k,ispec_wd) if (iglob_wd /= 0) then dummyx_loc(i,j,k) = dummyx_loc(i,j,k) + displ_wd(1, iglob_wd) @@ -153,7 +180,7 @@ subroutine add_displacement_discontinuity_element(ispec, dummyx_loc, & end subroutine add_displacement_discontinuity_element subroutine add_traction_discontinuity(accel, nglob) - use specfem_par, only: CUSTOM_REAL, NGLLX, NGLLY, NGLLZ, NGLLSQUARE, & + use specfem_par, only: CUSTOM_REAL, NGLLSQUARE, & ibool, NDIM !use specfem_par_elastic, only: accel implicit none diff --git a/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py b/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py index c977bb21d..b5a12b242 100644 --- a/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py +++ b/utils/dynamic_rupture/dipping_fault_planar_dip10_kink.py @@ -140,6 +140,7 @@ def main(parameter): ##### MESHING fault #### cubit.cmd("imprint all") cubit.cmd("merge all") + cubit.cmd("vol all size "+str(h_size*sf)) cubit.cmd("mesh volume all")