[cig-commits] r20950 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: . START_UP src
surendra at geodynamics.org
surendra at geodynamics.org
Fri Oct 26 20:15:38 PDT 2012
Author: surendra
Date: 2012-10-26 20:15:37 -0700 (Fri, 26 Oct 2012)
New Revision: 20950
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_generate_databases.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90
Removed:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
seismo/3D/FAULT_SOURCE/branches/new_fault_db/START_UP/Makefile.in
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic_noDev.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
Log:
renamed the fault_solver module to fault_solver_dynamic. renamed fault_object module to fault_generate_databases. Added documentation about mesh generation with split nodes. Added documentation about sign convention of fault quantities
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/README_SPECFEM3D_FAULT 2012-10-27 03:15:37 UTC (rev 20950)
@@ -2,10 +2,12 @@
This package contains a version of the 3D spectral element code SPECFEM3D-SESAME (Komatitsch, Tromp et al)
modified by J.-P. Ampuero (Caltech Seismolab) and P. Galvez (ETH Zurich)
to model dynamic earthquake rupture on non-planar faults with slip-weakening friction.
-The main modifications are encapsulated in three new modules:
+The main modifications are encapsulated in the following new modules:
decompose_mesh_SCOTCH/fault_scotch.f90.
- src/fault_object.f90
- src/fault_solver.f90
+ src/fault_generate_databases.f90
+ src/fault_solver_dynamic.f90
+ src/fault_solver_kinematic.f90
+ src/fault_solver_common.f90
We also include examples and Matlab functions for post-processing and visualization.
This is a preliminary version, still under development and testing.
@@ -454,3 +456,46 @@
-N simulation_name
-m bae -M your at email.address
+
+
+
+
+VII. MESH GENERATION WITH SPLIT NODES
+---------------------------------------
+Vertical strike-slip embedded fault. For more complicated meshes, please refer to EXAMPLES
+Note : There has to be at least one webcut honoring the fault plane
+
+reset
+brick x 10 y 10 z 10
+webcut volume all with plane xplane
+webcut volume all with plane yplane
+webcut volume all with plane xplane offset 3
+webcut volume all with plane zplane offset 3
+webcut volume all with plane zplane offset -3
+imprint all
+merge all
+unmerge surf 160
+mesh vol all
+
+set node constraint off
+node in surf 168 move X 0 Y 0.01 Z 0
+node in surf 160 move X 0 Y -0.01 Z 0
+
+
+
+
+VIII. SIGN CONVENTION FOR FAULT QUANTITIES
+---------------------------------------
+Fault in X-Z plane
+Assuming side with y-coord +ve is given as first side during mesh generation with CUBIT,
+ strike : pointing to x<0
+ dip : pointing z<0
+ normal : pointing y<0
+
+slip = u(x,0-,z) - u(x,0+,z)
+slip-rate = v(x,0-,z) - v(x,0+,z)
+Traction = sigma . n
+
+sigma (normal stress) is negative in compression
+
+
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/START_UP/Makefile.in
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/START_UP/Makefile.in 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/START_UP/Makefile.in 2012-10-27 03:15:37 UTC (rev 20950)
@@ -72,7 +72,7 @@
$O/detect_surface.o \
$O/exit_mpi.o \
$O/get_absorbing_boundary.o \
- $O/fault_object.o\
+ $O/fault_generate_databases.o\
$O/get_coupling_surfaces.o \
$O/get_model.o \
$O/get_MPI.o \
@@ -155,7 +155,7 @@
$O/compute_gradient.o \
$O/compute_interpolated_dva.o \
$O/fault_solver_common.o \
- $O/fault_solver.o \
+ $O/fault_solver_dynamic.o \
$O/fault_solver_kinematic.o \
$O/initialize_simulation.o \
$O/read_mesh_databases.o \
@@ -300,7 +300,7 @@
$O/specfem3D.o: constants.h specfem3D.f90
${FCCOMPILE_NO_CHECK} -c -o $O/specfem3D.o specfem3D.f90
-$O/compute_forces_elastic_noDev.o: constants.h compute_forces_elastic_noDev.f90
+$O/compute_forces_elastic_noDev.o: constants.h compute_forces_elastic_noDev.f90 $O/fault_solver_dynamic.o
${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_elastic_noDev.o compute_forces_elastic_noDev.f90
$O/compute_forces_elastic_Dev.o: constants.h compute_forces_elastic_Dev.f90
@@ -309,13 +309,13 @@
$O/fault_solver_common.o: fault_solver_common.f90 constants.h $O/specfem3D_par.o
${FCCOMPILE_NO_CHECK} -c -o $O/fault_solver_common.o fault_solver_common.f90
-$O/fault_solver.o: fault_solver.f90 constants.h $O/specfem3D_par.o $O/fault_solver_common.o
- ${FCCOMPILE_NO_CHECK} -c -o $O/fault_solver.o fault_solver.f90
+$O/fault_solver_dynamic.o: fault_solver_dyanmic.f90 constants.h $O/specfem3D_par.o $O/fault_solver_common.o
+ ${FCCOMPILE_NO_CHECK} -c -o $O/fault_solver_dynamic.o fault_solver_dynamic.f90
$O/fault_solver_kinematic.o: fault_solver_kinematic.f90 constants.h $O/specfem3D_par.o $O/fault_solver_common.o
${FCCOMPILE_NO_CHECK} -c -o $O/fault_solver_kinematic.o fault_solver_kinematic.f90
-$O/compute_forces_elastic.o: constants.h compute_forces_elastic.f90 $O/fault_solver.o $O/fault_solver_kinematic.o
+$O/compute_forces_elastic.o: constants.h compute_forces_elastic.f90 $O/fault_solver_dynamic.o $O/fault_solver_kinematic.o
${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_elastic.o compute_forces_elastic.f90
$O/compute_forces_acoustic.o: constants.h compute_forces_acoustic.f90
@@ -372,7 +372,7 @@
$O/setup_sources_receivers.o: constants.h setup_sources_receivers.f90
${FCCOMPILE_NO_CHECK} -c -o $O/setup_sources_receivers.o setup_sources_receivers.f90
-$O/prepare_timerun.o: constants.h prepare_timerun.f90 $O/fault_solver.o $O/fault_solver_kinematic.o
+$O/prepare_timerun.o: constants.h prepare_timerun.f90 $O/fault_solver_dynamic.o $O/fault_solver_kinematic.o
${FCCOMPILE_NO_CHECK} -c -o $O/prepare_timerun.o prepare_timerun.f90
$O/iterate_time.o: constants.h iterate_time.f90
@@ -561,10 +561,10 @@
$O/create_regions_mesh_ext_par.o: create_regions_mesh_ext_par.f90
${FCCOMPILE_CHECK} -c -o $O/create_regions_mesh_ext_par.o create_regions_mesh_ext_par.f90
-$O/fault_object.o: fault_object.f90 $O/create_regions_mesh_ext_par.o $O/generate_databases.o
- ${FCCOMPILE_CHECK} -c -o $O/fault_object.o fault_object.f90
+$O/fault_generate_databases.o: fault_generate_databases.f90 $O/create_regions_mesh_ext_par.o $O/generate_databases.o
+ ${FCCOMPILE_CHECK} -c -o $O/fault_generate_databases.o fault_generate_databases.f90
-$O/create_regions_mesh.o: constants.h create_regions_mesh.f90 $O/fault_object.o $O/create_regions_mesh_ext_par.o
+$O/create_regions_mesh.o: constants.h create_regions_mesh.f90 $O/fault_generate_databases.o $O/create_regions_mesh_ext_par.o
${FCCOMPILE_CHECK} -c -o $O/create_regions_mesh.o create_regions_mesh.f90
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic.f90 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -32,7 +32,7 @@
use specfem_par_acoustic
use specfem_par_elastic
use specfem_par_poroelastic
- use fault_solver, only : bc_dynflt_set3d_all,SIMULATION_TYPE_DYN
+ use fault_solver_dynamic, only : bc_dynflt_set3d_all,SIMULATION_TYPE_DYN
use fault_solver_kinematic, only : bc_kinflt_set_all,SIMULATION_TYPE_KIN
implicit none
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic_noDev.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic_noDev.f90 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/compute_forces_elastic_noDev.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -59,7 +59,7 @@
use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM, &
NUM_REGIONS_ATTENUATION,N_SLS,SAVE_MOHO_MESH, &
ONE_THIRD,FOUR_THIRDS
- use fault_solver, only : Kelvin_Voigt_eta
+ use fault_solver_dynamic, only : Kelvin_Voigt_eta
implicit none
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/create_regions_mesh.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -53,7 +53,7 @@
! create the different regions of the mesh
use create_regions_mesh_ext_par
- use fault_object, only: fault_read_input,fault_setup, &
+ use fault_generate_databases, only: fault_read_input,fault_setup, &
fault_save_arrays,fault_save_arrays_test,&
nnodes_coords_open,nodes_coords_open,ANY_FAULT_IN_THIS_PROC,&
ANY_FAULT, PARALLEL_FAULT
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_generate_databases.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_generate_databases.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_generate_databases.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -0,0 +1,669 @@
+!=====================================================================
+!
+! s p e c f e m 3 d v e r s i o n 1 . 4
+! ---------------------------------------
+!
+! dimitri komatitsch and jeroen tromp
+! seismological laboratory - california institute of technology
+! (c) california institute of technology september 2006
+!
+! 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 2 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.
+!
+!===============================================================================================================
+
+! This module was written by:
+! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
+
+! Splitting fault nodes is done in CUBIT.
+
+module fault_generate_databases
+
+ use create_regions_mesh_ext_par, only: NGLLX,NGLLY,NGLLZ,NGLLSQUARE,NGNOD2D,NDIM,CUSTOM_REAL
+! these variables are defined in 'constants.h', which is included in create_regions_mesh_ext_par
+
+ implicit none
+ private
+
+ type fault_db_type
+ private
+ integer :: nspec=0,nglob=0
+ real(kind=CUSTOM_REAL) :: eta
+ integer, dimension(:), pointer:: ispec1, ispec2, ibulk1, ibulk2, iface1, iface2
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoordbulk1,ycoordbulk1,zcoordbulk1,xcoordbulk2,ycoordbulk2,zcoordbulk2
+ integer, dimension(:,:), pointer :: ibool1, ibool2, inodes1, inodes2
+ integer, dimension(:,:,:), pointer :: ijk1, ijk2
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer:: jacobian2Dw
+ real(kind=CUSTOM_REAL), dimension(:,:,:), pointer:: normal
+ end type fault_db_type
+
+ type(fault_db_type), allocatable, save :: fault_db(:)
+ ! fault_db(i) is the database of the i-th fault in the mesh
+ real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
+ double precision, allocatable, save :: nodes_coords_open(:,:)
+ integer, save :: nnodes_coords_open
+ logical, save :: ANY_FAULT_IN_THIS_PROC = .false.
+ logical, save :: ANY_FAULT = .false.
+ logical, parameter :: PARALLEL_FAULT = .false.
+
+ ! corners indices of reference cube faces
+ integer,dimension(3,4),parameter :: iface1_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
+ integer,dimension(3,4),parameter :: iface2_corner_ijk = &
+ reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
+ integer,dimension(3,4),parameter :: iface3_corner_ijk = &
+ reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
+ integer,dimension(3,4),parameter :: iface4_corner_ijk = &
+ reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
+ integer,dimension(3,4),parameter :: iface5_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
+ integer,dimension(3,4),parameter :: iface6_corner_ijk = &
+ reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
+ integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
+ reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
+ iface3_corner_ijk,iface4_corner_ijk, &
+ iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
+
+ public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, &
+ nodes_coords_open, nnodes_coords_open, ANY_FAULT_IN_THIS_PROC, ANY_FAULT, PARALLEL_FAULT
+
+contains
+
+!=================================================================================================================
+subroutine fault_read_input(prname,NDIM)
+
+ character(len=256), intent(in) :: prname
+ integer, intent(in) :: NDIM
+
+ integer :: nb,i,iflt,ier,nspec,dummy_node
+ integer, parameter :: IIN = 100
+
+ ! read fault input file
+ nb = 0
+ open(unit=IIN,file='DATA/FAULT/Par_file_faults',status='old',action='read',iostat=ier)
+ if (ier==0) then
+ read(IIN,*) nb
+ else
+ write(6,*) 'No faults in the domain'
+ write(6,*) 'Par_file_faults does not exist '
+ close(IIN)
+ end if
+
+ ANY_FAULT = (nb>0)
+ if (.not. ANY_FAULT) return
+
+ allocate(fault_db(nb))
+ do i=1,nb
+ read(IIN,*) fault_db(i)%eta
+ enddo
+ close(IIN)
+
+ ! read fault database file
+ open(unit=IIN,file=prname(1:len_trim(prname))//'Database_fault', &
+ status='old',action='read',form='formatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(IIN,*) 'error opening file: ',prname(1:len_trim(prname))//'Database_fault'
+ write(IIN,*) 'make sure file exists'
+ stop
+ endif
+
+ do iflt=1,size(fault_db)
+
+ read(IIN,*) nspec
+ fault_db(iflt)%nspec = nspec
+
+ if (nspec == 0) cycle
+
+ ANY_FAULT_IN_THIS_PROC = .true.
+
+ allocate(fault_db(iflt)%ispec1(nspec))
+ allocate(fault_db(iflt)%inodes1(4,nspec))
+ allocate(fault_db(iflt)%ispec2(nspec))
+ allocate(fault_db(iflt)%inodes2(4,nspec))
+
+ do i=1,nspec
+ read(IIN,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%inodes1(:,i)
+ enddo
+ do i=1,nspec
+ read(IIN,*) fault_db(iflt)%ispec2(i), fault_db(iflt)%inodes2(:,i)
+ enddo
+
+ ! loading ispec1 ispec2 iface1 iface2 of fault elements.
+! allocate(fault_db(iflt)%iface1(nspec))
+! allocate(fault_db(iflt)%iface2(nspec))
+! do i=1,fault_db(iflt)%nspec
+! read(IIN,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%ispec2(i), &
+! fault_db(iflt)%iface1(i), fault_db(iflt)%iface2(i)
+! enddo
+
+ enddo
+
+ ! read nodes_coords with fault crack open.
+ read(IIN,*) nnodes_coords_open
+ allocate(nodes_coords_open(NDIM,nnodes_coords_open))
+ do i = 1, nnodes_coords_open
+ read(IIN,*) dummy_node, nodes_coords_open(:,i)
+ enddo
+
+ close(IIN)
+
+end subroutine fault_read_input
+
+!==================================================================================================================
+subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
+ xstore,ystore,zstore,nspec,nglob,myrank)
+
+ integer, intent(in) :: nspec
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(inout) :: xstore,ystore,zstore
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: ibool
+ integer, intent(in) :: myrank
+ integer, intent(in) :: nglob
+ integer, intent(in) :: nnodes_ext_mesh
+ double precision, dimension(NDIM,nnodes_ext_mesh),intent(in) :: nodes_coords_ext_mesh
+
+ integer :: iflt
+
+ if (.not. ANY_FAULT_IN_THIS_PROC ) return
+
+ do iflt=1,size(fault_db)
+ !NOTE: the small fault opening in *_dummy does not affect this subroutine (see get_element_face_id)
+ call setup_iface(fault_db(iflt),nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool)
+
+ ! saving gll indices for each fault face, needed for ibulks
+ call setup_ijk(fault_db(iflt))
+
+ ! ibools = mapping from local indices on the fault (GLL index, element index)
+ ! to global indices on the fault
+ call setup_ibools(fault_db(iflt),xstore,ystore,zstore,nspec,fault_db(iflt)%nspec*NGLLSQUARE)
+
+ ! ibulks = mapping global indices of fault nodes
+ ! from global indices on the fault to global indices on the bulk
+ call setup_ibulks(fault_db(iflt),ibool,nspec)
+
+ ! close the fault in (xyz)store_dummy
+ call close_fault(fault_db(iflt))
+
+ call setup_Kelvin_Voigt_eta(fault_db(iflt),nspec)
+
+ call save_fault_xyzcoord_ibulk(fault_db(iflt))
+
+ call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob,myrank)
+
+ enddo
+
+end subroutine fault_setup
+
+
+!=============================================================================================================
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ ! sets face id of reference element associated with this face
+subroutine setup_iface(fdb,nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
+
+ type(fault_db_type), intent(inout) :: fdb
+ integer, intent(in) :: nnodes_ext_mesh,nspec,nglob
+ double precision, dimension(NDIM,nnodes_ext_mesh), intent(in) :: nodes_coords_ext_mesh
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+ real(kind=CUSTOM_REAL), dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ integer :: icorner,e
+
+ allocate(fdb%iface1(fdb%nspec))
+ allocate(fdb%iface2(fdb%nspec))
+ do e=1,fdb%nspec
+ ! side 1
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,fdb%inodes1(icorner,e))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,fdb%inodes1(icorner,e))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,fdb%inodes1(icorner,e))
+ enddo
+ call get_element_face_id(fdb%ispec1(e),xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ fdb%iface1(e))
+ ! side 2
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,fdb%inodes2(icorner,e))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,fdb%inodes2(icorner,e))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,fdb%inodes2(icorner,e))
+ enddo
+ call get_element_face_id(fdb%ispec2(e),xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ fdb%iface2(e))
+ enddo
+
+end subroutine setup_iface
+
+!=============================================================================================================
+subroutine setup_ijk(fdb)
+
+ type(fault_db_type), intent(inout) :: fdb
+
+ integer :: e,i,j,igll
+ integer :: ijk_face1(3,NGLLX,NGLLY), ijk_face2(3,NGLLX,NGLLY)
+
+ allocate(fdb%ijk1(3,NGLLX*NGLLY,fdb%nspec))
+ allocate(fdb%ijk2(3,NGLLX*NGLLY,fdb%nspec))
+
+ do e=1,fdb%nspec
+ call get_element_face_gll_indices(fdb%iface1(e),ijk_face1,NGLLX,NGLLY)
+ call get_element_face_gll_indices(fdb%iface2(e),ijk_face2,NGLLX,NGLLY)
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ igll = igll + 1
+ fdb%ijk1(:,igll,e)=ijk_face1(:,i,j)
+ fdb%ijk2(:,igll,e)=ijk_face2(:,i,j)
+ enddo
+ enddo
+ enddo
+
+end subroutine setup_ijk
+
+!=============================================================================================================
+ subroutine setup_Kelvin_Voigt_eta(fdb,nspec)
+
+ type(fault_db_type), intent(in) :: fdb
+ integer, intent(in) :: nspec ! number of spectral elements in each block
+
+ if (fdb%eta > 0.0_CUSTOM_REAL) then
+ if (.not.allocated(Kelvin_Voigt_eta)) then
+ allocate(Kelvin_Voigt_eta(nspec))
+ Kelvin_Voigt_eta(:) = 0.0_CUSTOM_REAL
+ endif
+ Kelvin_Voigt_eta(fdb%ispec1) = fdb%eta
+ Kelvin_Voigt_eta(fdb%ispec2) = fdb%eta
+ endif
+
+ end subroutine
+
+!===============================================================================================================
+! The lexicographic ordering of node coordinates
+! guarantees that the fault nodes are
+! consistently ordered on both sides of the fault,
+! such that the K-th node of side 1 is facing the K-th node of side 2
+
+subroutine setup_ibools(fdb,xstore,ystore,zstore,nspec,npointot)
+
+ use generate_databases_par, only: nodes_coords_ext_mesh
+
+ type(fault_db_type), intent(inout) :: fdb
+ integer, intent(in) :: nspec,npointot
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: xstore,ystore,zstore
+
+ double precision :: xp(npointot),yp(npointot),zp(npointot),xmin,xmax
+ integer :: loc(npointot)
+ logical :: ifseg(npointot)
+ integer :: ispec,k,igll,ie,je,ke,e
+
+ xmin = minval(nodes_coords_ext_mesh(1,:))
+ xmax = maxval(nodes_coords_ext_mesh(1,:))
+
+ k = 0
+ do e = 1,fdb%nspec
+ ispec = fdb%ispec1(e)
+ do igll=1,NGLLSQUARE
+ ie=fdb%ijk1(1,igll,e)
+ je=fdb%ijk1(2,igll,e)
+ ke=fdb%ijk1(3,igll,e)
+ k = k+1
+ xp(k) = xstore(ie,je,ke,ispec)
+ yp(k) = ystore(ie,je,ke,ispec)
+ zp(k) = zstore(ie,je,ke,ispec)
+ enddo
+ enddo
+ allocate( fdb%ibool1(NGLLSQUARE,fdb%nspec) )
+ call get_global(fdb%nspec,xp,yp,zp,fdb%ibool1,loc,ifseg,fdb%nglob,npointot,xmin,xmax)
+
+! xp,yp,zp need to be recomputed on side 2
+! because they are generally not in the same order as on side 1,
+! because ispec1(e) is not necessarily facing ispec2(e).
+
+ k = 0
+ do e = 1,fdb%nspec
+ ispec = fdb%ispec2(e)
+ do igll=1,NGLLSQUARE
+ ie=fdb%ijk2(1,igll,e)
+ je=fdb%ijk2(2,igll,e)
+ ke=fdb%ijk2(3,igll,e)
+ k = k+1
+ xp(k) = xstore(ie,je,ke,ispec)
+ yp(k) = ystore(ie,je,ke,ispec)
+ zp(k) = zstore(ie,je,ke,ispec)
+ enddo
+ enddo
+ allocate( fdb%ibool2(NGLLSQUARE,fdb%nspec) )
+ call get_global(fdb%nspec,xp,yp,zp,fdb%ibool2,loc,ifseg,fdb%nglob,npointot,xmin,xmax)
+
+end subroutine setup_ibools
+
+
+!=================================================================================
+
+subroutine setup_ibulks(fdb,ibool,nspec)
+
+ type(fault_db_type), intent(inout) :: fdb
+ integer, intent(in) :: nspec, ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer :: e,k, K1, K2, ie,je,ke
+
+ allocate( fdb%ibulk1(fdb%nglob) )
+ allocate( fdb%ibulk2(fdb%nglob) )
+
+ do e=1, fdb%nspec
+ do k=1, NGLLSQUARE
+
+ ie=fdb%ijk1(1,k,e)
+ je=fdb%ijk1(2,k,e)
+ ke=fdb%ijk1(3,k,e)
+ K1= fdb%ibool1(k,e)
+ fdb%ibulk1(K1)=ibool(ie,je,ke,fdb%ispec1(e))
+
+ ie=fdb%ijk2(1,k,e)
+ je=fdb%ijk2(2,k,e)
+ ke=fdb%ijk2(3,k,e)
+ K2= fdb%ibool2(k,e)
+ fdb%ibulk2(K2)=ibool(ie,je,ke,fdb%ispec2(e))
+
+ enddo
+ enddo
+
+end subroutine setup_ibulks
+
+!=============================================================================================================
+! We only close *store_dummy. *store is close already in create_regions_mesh.f90
+! Fortunately only *store_dummy is needed to compute jacobians and normals
+
+subroutine close_fault(fdb)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
+
+ type(fault_db_type), intent(inout) :: fdb
+
+ integer :: i,K1,K2
+
+ do i=1,fdb%nglob
+ K1 = fdb%ibulk1(i)
+ K2 = fdb%ibulk2(i)
+ xstore_dummy(K1) = 0.5d0*( xstore_dummy(K1) + xstore_dummy(K2) )
+ xstore_dummy(K2) = xstore_dummy(K1)
+ ystore_dummy(K1) = 0.5d0*( ystore_dummy(K1) + ystore_dummy(K2) )
+ ystore_dummy(K2) = ystore_dummy(K1)
+ zstore_dummy(K1) = 0.5d0*( zstore_dummy(K1) + zstore_dummy(K2) )
+ zstore_dummy(K2) = zstore_dummy(K1)
+ enddo
+
+end subroutine close_fault
+
+!=================================================================================
+
+ subroutine save_fault_xyzcoord_ibulk(fdb)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
+
+ type(fault_db_type), intent(inout) :: fdb
+
+ integer :: K1, K2, i
+
+ allocate( fdb%xcoordbulk1(fdb%nglob) )
+ allocate( fdb%ycoordbulk1(fdb%nglob) )
+ allocate( fdb%zcoordbulk1(fdb%nglob) )
+ allocate( fdb%xcoordbulk2(fdb%nglob) )
+ allocate( fdb%ycoordbulk2(fdb%nglob) )
+ allocate( fdb%zcoordbulk2(fdb%nglob) )
+
+ do i=1, fdb%nglob
+ K1 =fdb%ibulk1(i)
+ K2 =fdb%ibulk2(i)
+ fdb%xcoordbulk1(i) = xstore_dummy(K1)
+ fdb%ycoordbulk1(i) = ystore_dummy(K1)
+ fdb%zcoordbulk1(i) = zstore_dummy(K1)
+
+ fdb%xcoordbulk2(i) = xstore_dummy(K2)
+ fdb%ycoordbulk2(i) = ystore_dummy(K2)
+ fdb%zcoordbulk2(i) = zstore_dummy(K2)
+ enddo
+
+ end subroutine save_fault_xyzcoord_ibulk
+
+
+!=================================================================================
+
+ subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob,myrank)
+
+ use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz
+
+ type(fault_db_type), intent(inout) :: fdb
+ integer, intent(in) :: nspec,nglob, ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer, intent(in) :: myrank
+
+ ! (assumes NGLLX=NGLLY=NGLLZ)
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
+ integer,dimension(NGNOD2D) :: iglob_corners_ref
+ integer :: ispec_flt,ispec,i,j,k,igll
+ integer :: iface_ref,icorner
+
+ allocate(fdb%normal(NDIM,NGLLSQUARE,fdb%nspec))
+ allocate(fdb%jacobian2Dw(NGLLSQUARE,fdb%nspec))
+
+ do ispec_flt=1,fdb%nspec
+
+ iface_ref= fdb%iface1(ispec_flt)
+ ispec = fdb%ispec1(ispec_flt)
+
+ ! takes indices of corners of reference face
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface_ref)
+ j = iface_all_corner_ijk(2,icorner,iface_ref)
+ k = iface_all_corner_ijk(3,icorner,iface_ref)
+
+ ! global reference indices
+ iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
+
+ ! reference corner coordinates
+ xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
+ ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
+ zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
+ enddo
+
+ ! gets face GLL 2Djacobian, weighted from element face
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
+
+ ! normal convention: points away from domain1, reference element.
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! directs normals such that they point outwards of element
+ call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! stores informations about this face
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! adds all gll points on that face
+ igll = igll + 1
+ ! stores weighted jacobian and normals
+ fdb%jacobian2Dw(igll,ispec_flt) = jacobian2Dw_face(i,j)
+ fdb%normal(:,igll,ispec_flt) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ enddo ! ispec_flt
+
+end subroutine setup_normal_jacobian
+
+!====================================================================================
+! saves all fault data in ASCII files for verification
+subroutine fault_save_arrays_test(prname)
+
+ character(len=256), intent(in) :: prname ! 'proc***'
+
+ integer, parameter :: IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+ integer :: nbfaults,iflt,ier
+ character(len=256) :: filename
+
+ if (.not.ANY_FAULT) return
+
+! saves mesh file proc***_fault_db.txt
+ filename = prname(1:len_trim(prname))//'fault_db.txt'
+ open(unit=IOUT,file=trim(filename),status='unknown',action='write',iostat=ier)
+ if( ier /= 0 ) stop 'error opening database proc######_external_mesh.bin'
+
+ nbfaults = size(fault_db)
+ write(IOUT,*) 'NBFAULTS = ',nbfaults
+ do iflt=1,nbfaults
+ write(IOUT,*) 'BEGIN FAULT # ',iflt
+ call save_one_fault_test(fault_db(iflt),IOUT)
+ write(IOUT,*) 'END FAULT # ',iflt
+ enddo
+ close(IOUT)
+
+end subroutine fault_save_arrays_test
+
+!-------------------------------------------------------------------------------------
+
+subroutine save_one_fault_test(f,IOUT)
+
+ type(fault_db_type), intent(in) :: f
+ integer, intent(in) :: IOUT
+
+ integer :: e,k
+ character(15) :: fmt1,fmt2
+
+ write(fmt1,'("(a,",I0,"(x,I7))")') NGLLSQUARE+1 ! fmt = (a,(NGLL^2+1)(x,I7))
+ write(fmt2,'("(a,",I0,"(x,F0.4))")') NGLLSQUARE+1 ! fmt = (a,(NGLL^2+1)(x,F0.16))
+
+ write(IOUT,*) 'NSPEC NGLOB NGLL = ',f%nspec,f%nglob,NGLLX
+ if (f%nspec==0) return
+ do e=1,f%nspec
+ write(IOUT,*) 'FLT_ELEM = ',e
+ write(IOUT,*) 'ISPEC1 ISPEC2 = ',f%ispec1(e),f%ispec2(e)
+ write(IOUT,fmt1) 'IBOOL1 = ',f%ibool1(:,e)
+ write(IOUT,fmt1) 'IBOOL2 = ',f%ibool2(:,e)
+ write(IOUT,fmt1) 'I1 = ',f%ijk1(1,:,e)
+ write(IOUT,fmt1) 'J1 = ',f%ijk1(2,:,e)
+ write(IOUT,fmt1) 'K1 = ',f%ijk1(3,:,e)
+ write(IOUT,fmt1) 'I2 = ',f%ijk2(1,:,e)
+ write(IOUT,fmt1) 'J2 = ',f%ijk2(2,:,e)
+ write(IOUT,fmt1) 'K2 = ',f%ijk2(3,:,e)
+ write(IOUT,fmt2) 'JAC2DW = ',f%jacobian2Dw(:,e)
+ write(IOUT,fmt2) 'N1 = ',f%normal(1,:,e)
+ write(IOUT,fmt2) 'N2 = ',f%normal(2,:,e)
+ write(IOUT,fmt2) 'N3 = ',f%normal(3,:,e)
+ enddo
+
+ write(IOUT,*) 'FLT_NODE IBULK1 IBULK2'
+ do k=1,f%nglob
+ write(IOUT,*) k,f%ibulk1(k),f%ibulk2(k)
+ enddo
+
+ write(IOUT,*) 'FLT_NODE xcoordbulk ycoordbulk zcoordbulk'
+ do k=1,f%nglob
+ write(IOUT,*) f%ibulk1(k),f%xcoordbulk1(k),f%ycoordbulk1(k),f%zcoordbulk1(k)
+ write(IOUT,*) f%ibulk2(k),f%xcoordbulk2(k),f%ycoordbulk2(k),f%zcoordbulk2(k)
+ enddo
+
+end subroutine save_one_fault_test
+
+!=================================================================================
+! saves fault data needed by the solver in binary files
+subroutine fault_save_arrays(prname)
+
+ character(len=256), intent(in) :: prname ! 'proc***'
+
+ integer, parameter :: IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+ integer :: nbfaults,iflt,ier
+ character(len=256) :: filename
+ integer :: size_Kelvin_Voigt
+
+ if (.not.ANY_FAULT) return
+! opening Kelvin_voig_eta.bin (Necessary for all processors , if number of fault elements = 0
+! then the file will be empty with
+ filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
+ open(unit=IOUT,file=trim(filename),status='unknown',action='write',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(IOUT) 'error opening file ',trim(filename)
+ ! JPA instead of IOUT this error message should go to an error log file
+ stop
+ endif
+
+! saves mesh file proc***_Kelvin_voigt_eta.bin
+ if (allocated(Kelvin_Voigt_eta)) then
+ size_Kelvin_Voigt = size(Kelvin_Voigt_eta)
+ else
+ size_Kelvin_Voigt = 0
+ endif
+ write(IOUT) size_Kelvin_Voigt
+ if (size_Kelvin_Voigt /= 0) Write(IOUT) Kelvin_Voigt_eta
+ Close(IOUT)
+
+! saves mesh file proc***_fault_db.bin
+ filename = prname(1:len_trim(prname))//'fault_db.bin'
+ open(unit=IOUT,file=trim(filename),status='unknown',action='write',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(IOUT,*) 'error opening file ',trim(filename)
+ ! JPA instead of IOUT this error message should go to an error log file
+ stop
+ endif
+
+ nbfaults = size(fault_db)
+ write(IOUT) nbfaults
+ do iflt=1,nbfaults
+ call save_one_fault_bin(fault_db(iflt),IOUT)
+ enddo
+ close(IOUT)
+
+
+end subroutine fault_save_arrays
+
+!----------------------------------------------
+
+subroutine save_one_fault_bin(f,IOUT)
+
+ type(fault_db_type), intent(in) :: f
+ integer, intent(in) :: IOUT
+
+ write(IOUT) f%nspec,f%nglob
+ if (f%nspec==0) return
+ write(IOUT) f%ibool1
+ write(IOUT) f%jacobian2Dw
+ write(IOUT) f%normal
+ write(IOUT) f%ibulk1
+ write(IOUT) f%ibulk2
+ write(IOUT) f%xcoordbulk1
+ write(IOUT) f%ycoordbulk1
+ write(IOUT) f%zcoordbulk1
+
+end subroutine save_one_fault_bin
+
+!------------------------------------------------
+
+
+end module fault_generate_databases
Deleted: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -1,669 +0,0 @@
-!=====================================================================
-!
-! s p e c f e m 3 d v e r s i o n 1 . 4
-! ---------------------------------------
-!
-! dimitri komatitsch and jeroen tromp
-! seismological laboratory - california institute of technology
-! (c) california institute of technology september 2006
-!
-! 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 2 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.
-!
-!===============================================================================================================
-
-! This module was written by:
-! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
-
-! Splitting fault nodes is done in CUBIT.
-
-module fault_object
-
- use create_regions_mesh_ext_par, only: NGLLX,NGLLY,NGLLZ,NGLLSQUARE,NGNOD2D,NDIM,CUSTOM_REAL
-! these variables are defined in 'constants.h', which is included in create_regions_mesh_ext_par
-
- implicit none
- private
-
- type fault_db_type
- private
- integer :: nspec=0,nglob=0
- real(kind=CUSTOM_REAL) :: eta
- integer, dimension(:), pointer:: ispec1, ispec2, ibulk1, ibulk2, iface1, iface2
- real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoordbulk1,ycoordbulk1,zcoordbulk1,xcoordbulk2,ycoordbulk2,zcoordbulk2
- integer, dimension(:,:), pointer :: ibool1, ibool2, inodes1, inodes2
- integer, dimension(:,:,:), pointer :: ijk1, ijk2
- real(kind=CUSTOM_REAL), dimension(:,:), pointer:: jacobian2Dw
- real(kind=CUSTOM_REAL), dimension(:,:,:), pointer:: normal
- end type fault_db_type
-
- type(fault_db_type), allocatable, save :: fault_db(:)
- ! fault_db(i) is the database of the i-th fault in the mesh
- real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
- double precision, allocatable, save :: nodes_coords_open(:,:)
- integer, save :: nnodes_coords_open
- logical, save :: ANY_FAULT_IN_THIS_PROC = .false.
- logical, save :: ANY_FAULT = .false.
- logical, parameter :: PARALLEL_FAULT = .false.
-
- ! corners indices of reference cube faces
- integer,dimension(3,4),parameter :: iface1_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
- integer,dimension(3,4),parameter :: iface2_corner_ijk = &
- reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
- integer,dimension(3,4),parameter :: iface3_corner_ijk = &
- reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
- integer,dimension(3,4),parameter :: iface4_corner_ijk = &
- reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
- integer,dimension(3,4),parameter :: iface5_corner_ijk = &
- reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
- integer,dimension(3,4),parameter :: iface6_corner_ijk = &
- reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
- integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
- reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
- iface3_corner_ijk,iface4_corner_ijk, &
- iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
-
- public :: fault_read_input, fault_setup, fault_save_arrays_test, fault_save_arrays, &
- nodes_coords_open, nnodes_coords_open, ANY_FAULT_IN_THIS_PROC, ANY_FAULT, PARALLEL_FAULT
-
-contains
-
-!=================================================================================================================
-subroutine fault_read_input(prname,NDIM)
-
- character(len=256), intent(in) :: prname
- integer, intent(in) :: NDIM
-
- integer :: nb,i,iflt,ier,nspec,dummy_node
- integer, parameter :: IIN = 100
-
- ! read fault input file
- nb = 0
- open(unit=IIN,file='DATA/FAULT/Par_file_faults',status='old',action='read',iostat=ier)
- if (ier==0) then
- read(IIN,*) nb
- else
- write(6,*) 'No faults in the domain'
- write(6,*) 'Par_file_faults does not exist '
- close(IIN)
- end if
-
- ANY_FAULT = (nb>0)
- if (.not. ANY_FAULT) return
-
- allocate(fault_db(nb))
- do i=1,nb
- read(IIN,*) fault_db(i)%eta
- enddo
- close(IIN)
-
- ! read fault database file
- open(unit=IIN,file=prname(1:len_trim(prname))//'Database_fault', &
- status='old',action='read',form='formatted',iostat=ier)
- if( ier /= 0 ) then
- write(IIN,*) 'error opening file: ',prname(1:len_trim(prname))//'Database_fault'
- write(IIN,*) 'make sure file exists'
- stop
- endif
-
- do iflt=1,size(fault_db)
-
- read(IIN,*) nspec
- fault_db(iflt)%nspec = nspec
-
- if (nspec == 0) cycle
-
- ANY_FAULT_IN_THIS_PROC = .true.
-
- allocate(fault_db(iflt)%ispec1(nspec))
- allocate(fault_db(iflt)%inodes1(4,nspec))
- allocate(fault_db(iflt)%ispec2(nspec))
- allocate(fault_db(iflt)%inodes2(4,nspec))
-
- do i=1,nspec
- read(IIN,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%inodes1(:,i)
- enddo
- do i=1,nspec
- read(IIN,*) fault_db(iflt)%ispec2(i), fault_db(iflt)%inodes2(:,i)
- enddo
-
- ! loading ispec1 ispec2 iface1 iface2 of fault elements.
-! allocate(fault_db(iflt)%iface1(nspec))
-! allocate(fault_db(iflt)%iface2(nspec))
-! do i=1,fault_db(iflt)%nspec
-! read(IIN,*) fault_db(iflt)%ispec1(i), fault_db(iflt)%ispec2(i), &
-! fault_db(iflt)%iface1(i), fault_db(iflt)%iface2(i)
-! enddo
-
- enddo
-
- ! read nodes_coords with fault crack open.
- read(IIN,*) nnodes_coords_open
- allocate(nodes_coords_open(NDIM,nnodes_coords_open))
- do i = 1, nnodes_coords_open
- read(IIN,*) dummy_node, nodes_coords_open(:,i)
- enddo
-
- close(IIN)
-
-end subroutine fault_read_input
-
-!==================================================================================================================
-subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, &
- xstore,ystore,zstore,nspec,nglob,myrank)
-
- integer, intent(in) :: nspec
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(inout) :: xstore,ystore,zstore
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: ibool
- integer, intent(in) :: myrank
- integer, intent(in) :: nglob
- integer, intent(in) :: nnodes_ext_mesh
- double precision, dimension(NDIM,nnodes_ext_mesh),intent(in) :: nodes_coords_ext_mesh
-
- integer :: iflt
-
- if (.not. ANY_FAULT_IN_THIS_PROC ) return
-
- do iflt=1,size(fault_db)
- !NOTE: the small fault opening in *_dummy does not affect this subroutine (see get_element_face_id)
- call setup_iface(fault_db(iflt),nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool)
-
- ! saving gll indices for each fault face, needed for ibulks
- call setup_ijk(fault_db(iflt))
-
- ! ibools = mapping from local indices on the fault (GLL index, element index)
- ! to global indices on the fault
- call setup_ibools(fault_db(iflt),xstore,ystore,zstore,nspec,fault_db(iflt)%nspec*NGLLSQUARE)
-
- ! ibulks = mapping global indices of fault nodes
- ! from global indices on the fault to global indices on the bulk
- call setup_ibulks(fault_db(iflt),ibool,nspec)
-
- ! close the fault in (xyz)store_dummy
- call close_fault(fault_db(iflt))
-
- call setup_Kelvin_Voigt_eta(fault_db(iflt),nspec)
-
- call save_fault_xyzcoord_ibulk(fault_db(iflt))
-
- call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob,myrank)
-
- enddo
-
-end subroutine fault_setup
-
-
-!=============================================================================================================
- ! looks for i,j,k indices of GLL points on boundary face
- ! determines element face by given CUBIT corners
- ! sets face id of reference element associated with this face
-subroutine setup_iface(fdb,nnodes_ext_mesh,nodes_coords_ext_mesh,nspec,nglob,ibool)
-
- use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
-
- type(fault_db_type), intent(inout) :: fdb
- integer, intent(in) :: nnodes_ext_mesh,nspec,nglob
- double precision, dimension(NDIM,nnodes_ext_mesh), intent(in) :: nodes_coords_ext_mesh
- integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
-
- real(kind=CUSTOM_REAL), dimension(NGNOD2D) :: xcoord,ycoord,zcoord
- integer :: icorner,e
-
- allocate(fdb%iface1(fdb%nspec))
- allocate(fdb%iface2(fdb%nspec))
- do e=1,fdb%nspec
- ! side 1
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,fdb%inodes1(icorner,e))
- ycoord(icorner) = nodes_coords_ext_mesh(2,fdb%inodes1(icorner,e))
- zcoord(icorner) = nodes_coords_ext_mesh(3,fdb%inodes1(icorner,e))
- enddo
- call get_element_face_id(fdb%ispec1(e),xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- fdb%iface1(e))
- ! side 2
- do icorner=1,NGNOD2D
- xcoord(icorner) = nodes_coords_ext_mesh(1,fdb%inodes2(icorner,e))
- ycoord(icorner) = nodes_coords_ext_mesh(2,fdb%inodes2(icorner,e))
- zcoord(icorner) = nodes_coords_ext_mesh(3,fdb%inodes2(icorner,e))
- enddo
- call get_element_face_id(fdb%ispec2(e),xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- fdb%iface2(e))
- enddo
-
-end subroutine setup_iface
-
-!=============================================================================================================
-subroutine setup_ijk(fdb)
-
- type(fault_db_type), intent(inout) :: fdb
-
- integer :: e,i,j,igll
- integer :: ijk_face1(3,NGLLX,NGLLY), ijk_face2(3,NGLLX,NGLLY)
-
- allocate(fdb%ijk1(3,NGLLX*NGLLY,fdb%nspec))
- allocate(fdb%ijk2(3,NGLLX*NGLLY,fdb%nspec))
-
- do e=1,fdb%nspec
- call get_element_face_gll_indices(fdb%iface1(e),ijk_face1,NGLLX,NGLLY)
- call get_element_face_gll_indices(fdb%iface2(e),ijk_face2,NGLLX,NGLLY)
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- igll = igll + 1
- fdb%ijk1(:,igll,e)=ijk_face1(:,i,j)
- fdb%ijk2(:,igll,e)=ijk_face2(:,i,j)
- enddo
- enddo
- enddo
-
-end subroutine setup_ijk
-
-!=============================================================================================================
- subroutine setup_Kelvin_Voigt_eta(fdb,nspec)
-
- type(fault_db_type), intent(in) :: fdb
- integer, intent(in) :: nspec ! number of spectral elements in each block
-
- if (fdb%eta > 0.0_CUSTOM_REAL) then
- if (.not.allocated(Kelvin_Voigt_eta)) then
- allocate(Kelvin_Voigt_eta(nspec))
- Kelvin_Voigt_eta(:) = 0.0_CUSTOM_REAL
- endif
- Kelvin_Voigt_eta(fdb%ispec1) = fdb%eta
- Kelvin_Voigt_eta(fdb%ispec2) = fdb%eta
- endif
-
- end subroutine
-
-!===============================================================================================================
-! The lexicographic ordering of node coordinates
-! guarantees that the fault nodes are
-! consistently ordered on both sides of the fault,
-! such that the K-th node of side 1 is facing the K-th node of side 2
-
-subroutine setup_ibools(fdb,xstore,ystore,zstore,nspec,npointot)
-
- use generate_databases_par, only: nodes_coords_ext_mesh
-
- type(fault_db_type), intent(inout) :: fdb
- integer, intent(in) :: nspec,npointot
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: xstore,ystore,zstore
-
- double precision :: xp(npointot),yp(npointot),zp(npointot),xmin,xmax
- integer :: loc(npointot)
- logical :: ifseg(npointot)
- integer :: ispec,k,igll,ie,je,ke,e
-
- xmin = minval(nodes_coords_ext_mesh(1,:))
- xmax = maxval(nodes_coords_ext_mesh(1,:))
-
- k = 0
- do e = 1,fdb%nspec
- ispec = fdb%ispec1(e)
- do igll=1,NGLLSQUARE
- ie=fdb%ijk1(1,igll,e)
- je=fdb%ijk1(2,igll,e)
- ke=fdb%ijk1(3,igll,e)
- k = k+1
- xp(k) = xstore(ie,je,ke,ispec)
- yp(k) = ystore(ie,je,ke,ispec)
- zp(k) = zstore(ie,je,ke,ispec)
- enddo
- enddo
- allocate( fdb%ibool1(NGLLSQUARE,fdb%nspec) )
- call get_global(fdb%nspec,xp,yp,zp,fdb%ibool1,loc,ifseg,fdb%nglob,npointot,xmin,xmax)
-
-! xp,yp,zp need to be recomputed on side 2
-! because they are generally not in the same order as on side 1,
-! because ispec1(e) is not necessarily facing ispec2(e).
-
- k = 0
- do e = 1,fdb%nspec
- ispec = fdb%ispec2(e)
- do igll=1,NGLLSQUARE
- ie=fdb%ijk2(1,igll,e)
- je=fdb%ijk2(2,igll,e)
- ke=fdb%ijk2(3,igll,e)
- k = k+1
- xp(k) = xstore(ie,je,ke,ispec)
- yp(k) = ystore(ie,je,ke,ispec)
- zp(k) = zstore(ie,je,ke,ispec)
- enddo
- enddo
- allocate( fdb%ibool2(NGLLSQUARE,fdb%nspec) )
- call get_global(fdb%nspec,xp,yp,zp,fdb%ibool2,loc,ifseg,fdb%nglob,npointot,xmin,xmax)
-
-end subroutine setup_ibools
-
-
-!=================================================================================
-
-subroutine setup_ibulks(fdb,ibool,nspec)
-
- type(fault_db_type), intent(inout) :: fdb
- integer, intent(in) :: nspec, ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
- integer :: e,k, K1, K2, ie,je,ke
-
- allocate( fdb%ibulk1(fdb%nglob) )
- allocate( fdb%ibulk2(fdb%nglob) )
-
- do e=1, fdb%nspec
- do k=1, NGLLSQUARE
-
- ie=fdb%ijk1(1,k,e)
- je=fdb%ijk1(2,k,e)
- ke=fdb%ijk1(3,k,e)
- K1= fdb%ibool1(k,e)
- fdb%ibulk1(K1)=ibool(ie,je,ke,fdb%ispec1(e))
-
- ie=fdb%ijk2(1,k,e)
- je=fdb%ijk2(2,k,e)
- ke=fdb%ijk2(3,k,e)
- K2= fdb%ibool2(k,e)
- fdb%ibulk2(K2)=ibool(ie,je,ke,fdb%ispec2(e))
-
- enddo
- enddo
-
-end subroutine setup_ibulks
-
-!=============================================================================================================
-! We only close *store_dummy. *store is close already in create_regions_mesh.f90
-! Fortunately only *store_dummy is needed to compute jacobians and normals
-
-subroutine close_fault(fdb)
-
- use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
-
- type(fault_db_type), intent(inout) :: fdb
-
- integer :: i,K1,K2
-
- do i=1,fdb%nglob
- K1 = fdb%ibulk1(i)
- K2 = fdb%ibulk2(i)
- xstore_dummy(K1) = 0.5d0*( xstore_dummy(K1) + xstore_dummy(K2) )
- xstore_dummy(K2) = xstore_dummy(K1)
- ystore_dummy(K1) = 0.5d0*( ystore_dummy(K1) + ystore_dummy(K2) )
- ystore_dummy(K2) = ystore_dummy(K1)
- zstore_dummy(K1) = 0.5d0*( zstore_dummy(K1) + zstore_dummy(K2) )
- zstore_dummy(K2) = zstore_dummy(K1)
- enddo
-
-end subroutine close_fault
-
-!=================================================================================
-
- subroutine save_fault_xyzcoord_ibulk(fdb)
-
- use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy
-
- type(fault_db_type), intent(inout) :: fdb
-
- integer :: K1, K2, i
-
- allocate( fdb%xcoordbulk1(fdb%nglob) )
- allocate( fdb%ycoordbulk1(fdb%nglob) )
- allocate( fdb%zcoordbulk1(fdb%nglob) )
- allocate( fdb%xcoordbulk2(fdb%nglob) )
- allocate( fdb%ycoordbulk2(fdb%nglob) )
- allocate( fdb%zcoordbulk2(fdb%nglob) )
-
- do i=1, fdb%nglob
- K1 =fdb%ibulk1(i)
- K2 =fdb%ibulk2(i)
- fdb%xcoordbulk1(i) = xstore_dummy(K1)
- fdb%ycoordbulk1(i) = ystore_dummy(K1)
- fdb%zcoordbulk1(i) = zstore_dummy(K1)
-
- fdb%xcoordbulk2(i) = xstore_dummy(K2)
- fdb%ycoordbulk2(i) = ystore_dummy(K2)
- fdb%zcoordbulk2(i) = zstore_dummy(K2)
- enddo
-
- end subroutine save_fault_xyzcoord_ibulk
-
-
-!=================================================================================
-
- subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob,myrank)
-
- use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz
-
- type(fault_db_type), intent(inout) :: fdb
- integer, intent(in) :: nspec,nglob, ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
- integer, intent(in) :: myrank
-
- ! (assumes NGLLX=NGLLY=NGLLZ)
- real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
- real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
- real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
- integer,dimension(NGNOD2D) :: iglob_corners_ref
- integer :: ispec_flt,ispec,i,j,k,igll
- integer :: iface_ref,icorner
-
- allocate(fdb%normal(NDIM,NGLLSQUARE,fdb%nspec))
- allocate(fdb%jacobian2Dw(NGLLSQUARE,fdb%nspec))
-
- do ispec_flt=1,fdb%nspec
-
- iface_ref= fdb%iface1(ispec_flt)
- ispec = fdb%ispec1(ispec_flt)
-
- ! takes indices of corners of reference face
- do icorner = 1,NGNOD2D
- i = iface_all_corner_ijk(1,icorner,iface_ref)
- j = iface_all_corner_ijk(2,icorner,iface_ref)
- k = iface_all_corner_ijk(3,icorner,iface_ref)
-
- ! global reference indices
- iglob_corners_ref(icorner) = ibool(i,j,k,ispec)
-
- ! reference corner coordinates
- xcoord(icorner) = xstore_dummy(iglob_corners_ref(icorner))
- ycoord(icorner) = ystore_dummy(iglob_corners_ref(icorner))
- zcoord(icorner) = zstore_dummy(iglob_corners_ref(icorner))
- enddo
-
- ! gets face GLL 2Djacobian, weighted from element face
- call get_jacobian_boundary_face(myrank,nspec, &
- xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY)
-
- ! normal convention: points away from domain1, reference element.
- do j=1,NGLLY
- do i=1,NGLLX
- ! directs normals such that they point outwards of element
- call get_element_face_normal(ispec,iface_ref,xcoord,ycoord,zcoord, &
- ibool,nspec,nglob, &
- xstore_dummy,ystore_dummy,zstore_dummy, &
- normal_face(:,i,j) )
- enddo
- enddo
-
- ! stores informations about this face
- igll = 0
- do j=1,NGLLY
- do i=1,NGLLX
- ! adds all gll points on that face
- igll = igll + 1
- ! stores weighted jacobian and normals
- fdb%jacobian2Dw(igll,ispec_flt) = jacobian2Dw_face(i,j)
- fdb%normal(:,igll,ispec_flt) = normal_face(:,i,j)
- enddo
- enddo
-
- enddo ! ispec_flt
-
-end subroutine setup_normal_jacobian
-
-!====================================================================================
-! saves all fault data in ASCII files for verification
-subroutine fault_save_arrays_test(prname)
-
- character(len=256), intent(in) :: prname ! 'proc***'
-
- integer, parameter :: IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
- integer :: nbfaults,iflt,ier
- character(len=256) :: filename
-
- if (.not.ANY_FAULT) return
-
-! saves mesh file proc***_fault_db.txt
- filename = prname(1:len_trim(prname))//'fault_db.txt'
- open(unit=IOUT,file=trim(filename),status='unknown',action='write',iostat=ier)
- if( ier /= 0 ) stop 'error opening database proc######_external_mesh.bin'
-
- nbfaults = size(fault_db)
- write(IOUT,*) 'NBFAULTS = ',nbfaults
- do iflt=1,nbfaults
- write(IOUT,*) 'BEGIN FAULT # ',iflt
- call save_one_fault_test(fault_db(iflt),IOUT)
- write(IOUT,*) 'END FAULT # ',iflt
- enddo
- close(IOUT)
-
-end subroutine fault_save_arrays_test
-
-!-------------------------------------------------------------------------------------
-
-subroutine save_one_fault_test(f,IOUT)
-
- type(fault_db_type), intent(in) :: f
- integer, intent(in) :: IOUT
-
- integer :: e,k
- character(15) :: fmt1,fmt2
-
- write(fmt1,'("(a,",I0,"(x,I7))")') NGLLSQUARE+1 ! fmt = (a,(NGLL^2+1)(x,I7))
- write(fmt2,'("(a,",I0,"(x,F0.4))")') NGLLSQUARE+1 ! fmt = (a,(NGLL^2+1)(x,F0.16))
-
- write(IOUT,*) 'NSPEC NGLOB NGLL = ',f%nspec,f%nglob,NGLLX
- if (f%nspec==0) return
- do e=1,f%nspec
- write(IOUT,*) 'FLT_ELEM = ',e
- write(IOUT,*) 'ISPEC1 ISPEC2 = ',f%ispec1(e),f%ispec2(e)
- write(IOUT,fmt1) 'IBOOL1 = ',f%ibool1(:,e)
- write(IOUT,fmt1) 'IBOOL2 = ',f%ibool2(:,e)
- write(IOUT,fmt1) 'I1 = ',f%ijk1(1,:,e)
- write(IOUT,fmt1) 'J1 = ',f%ijk1(2,:,e)
- write(IOUT,fmt1) 'K1 = ',f%ijk1(3,:,e)
- write(IOUT,fmt1) 'I2 = ',f%ijk2(1,:,e)
- write(IOUT,fmt1) 'J2 = ',f%ijk2(2,:,e)
- write(IOUT,fmt1) 'K2 = ',f%ijk2(3,:,e)
- write(IOUT,fmt2) 'JAC2DW = ',f%jacobian2Dw(:,e)
- write(IOUT,fmt2) 'N1 = ',f%normal(1,:,e)
- write(IOUT,fmt2) 'N2 = ',f%normal(2,:,e)
- write(IOUT,fmt2) 'N3 = ',f%normal(3,:,e)
- enddo
-
- write(IOUT,*) 'FLT_NODE IBULK1 IBULK2'
- do k=1,f%nglob
- write(IOUT,*) k,f%ibulk1(k),f%ibulk2(k)
- enddo
-
- write(IOUT,*) 'FLT_NODE xcoordbulk ycoordbulk zcoordbulk'
- do k=1,f%nglob
- write(IOUT,*) f%ibulk1(k),f%xcoordbulk1(k),f%ycoordbulk1(k),f%zcoordbulk1(k)
- write(IOUT,*) f%ibulk2(k),f%xcoordbulk2(k),f%ycoordbulk2(k),f%zcoordbulk2(k)
- enddo
-
-end subroutine save_one_fault_test
-
-!=================================================================================
-! saves fault data needed by the solver in binary files
-subroutine fault_save_arrays(prname)
-
- character(len=256), intent(in) :: prname ! 'proc***'
-
- integer, parameter :: IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
- integer :: nbfaults,iflt,ier
- character(len=256) :: filename
- integer :: size_Kelvin_Voigt
-
- if (.not.ANY_FAULT) return
-! opening Kelvin_voig_eta.bin (Necessary for all processors , if number of fault elements = 0
-! then the file will be empty with
- filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
- open(unit=IOUT,file=trim(filename),status='unknown',action='write',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(IOUT) 'error opening file ',trim(filename)
- ! JPA instead of IOUT this error message should go to an error log file
- stop
- endif
-
-! saves mesh file proc***_Kelvin_voigt_eta.bin
- if (allocated(Kelvin_Voigt_eta)) then
- size_Kelvin_Voigt = size(Kelvin_Voigt_eta)
- else
- size_Kelvin_Voigt = 0
- endif
- write(IOUT) size_Kelvin_Voigt
- if (size_Kelvin_Voigt /= 0) Write(IOUT) Kelvin_Voigt_eta
- Close(IOUT)
-
-! saves mesh file proc***_fault_db.bin
- filename = prname(1:len_trim(prname))//'fault_db.bin'
- open(unit=IOUT,file=trim(filename),status='unknown',action='write',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(IOUT,*) 'error opening file ',trim(filename)
- ! JPA instead of IOUT this error message should go to an error log file
- stop
- endif
-
- nbfaults = size(fault_db)
- write(IOUT) nbfaults
- do iflt=1,nbfaults
- call save_one_fault_bin(fault_db(iflt),IOUT)
- enddo
- close(IOUT)
-
-
-end subroutine fault_save_arrays
-
-!----------------------------------------------
-
-subroutine save_one_fault_bin(f,IOUT)
-
- type(fault_db_type), intent(in) :: f
- integer, intent(in) :: IOUT
-
- write(IOUT) f%nspec,f%nglob
- if (f%nspec==0) return
- write(IOUT) f%ibool1
- write(IOUT) f%jacobian2Dw
- write(IOUT) f%normal
- write(IOUT) f%ibulk1
- write(IOUT) f%ibulk2
- write(IOUT) f%xcoordbulk1
- write(IOUT) f%ycoordbulk1
- write(IOUT) f%zcoordbulk1
-
-end subroutine save_one_fault_bin
-
-!------------------------------------------------
-
-
-end module fault_object
Deleted: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-27 02:06:40 UTC (rev 20949)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -1,1457 +0,0 @@
-!=====================================================================
-!
-! s p e c f e m 3 d v e r s i o n 1 . 4
-! ---------------------------------------
-!
-! dimitri komatitsch and jeroen tromp
-! seismological laboratory - california institute of technology
-! (c) california institute of technology september 2006
-!
-! 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 2 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.
-!
-!===============================================================================================================
-
-! This module was written by:
-! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
-! Surendra Nadh Somala : Added Heterogenous initial stress capabilities (based on TPV16)
-! Surendra Nadh Somala : Added Rate and State Friction
-! Somala/Ampuero : fault parallelization
-
-module fault_solver
-
- use fault_solver_common
-
- implicit none
-
- include 'constants.h'
-
- private
-
- ! outputs on selected fault nodes at every time step:
- ! slip, slip velocity, fault stresses
- type dataT_type
- integer :: npoin=0,npoin_local=0
- integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
- integer, dimension(:,:), pointer :: iglob_all=>null()
- integer, dimension(:), pointer :: islice=>null(),glob_indx=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: dist=>null()
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: dist_all=>null()
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1=>null(),v1=>null(),t1=>null(), &
- d2=>null(),v2=>null(),t2=>null(), &
- t3=>null(),theta=>null()
- character(len=70), dimension(:), pointer :: name=>null()
- end type dataT_type
-
-
- ! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
- ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
- ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
- ! process zone time = first time when slip = Dc
- type dataXZ_type
- real(kind=CUSTOM_REAL), dimension(:), pointer :: stg=>null(), sta=>null(), d1=>null(), d2=>null(), v1=>null(), v2=>null(), &
- t1=>null(), t2=>null(), t3=>null(), tRUP=>null(), tPZ=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord=>null(), ycoord=>null(), zcoord=>null()
- integer :: npoin=0
- end type dataXZ_type
-
- type swf_type
- private
- integer :: kind
- logical :: healing = .false.
- real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), &
- theta=>null(), T=>null(), C=>null()
- end type swf_type
-
- type rsf_type
- private
- integer :: StateLaw = 1 ! 1=ageing law, 2=slip law
- real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), &
- V_init=>null(), &
- a=>null(), b=>null(), theta=>null(), &
- T=>null(), C=>null(), &
- fw=>null(), Vw=>null()
- end type rsf_type
-
- type, extends (fault_type) :: bc_dynflt_type
- private
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null(), Fload=>null()
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ,dataXZ_all
- type(swf_type), pointer :: swf => null()
- type(rsf_type), pointer :: rsf => null()
- logical :: allow_opening = .false. ! default : do not allow opening
- end type bc_dynflt_type
-
- type(bc_dynflt_type), allocatable, save :: faults(:)
-
- !slip velocity threshold for healing
- !WARNING: not very robust
- real(kind=CUSTOM_REAL), save :: V_HEALING
-
- !slip velocity threshold for definition of rupture front
- real(kind=CUSTOM_REAL), save :: V_RUPT
-
- !Number of time steps defined by the user : NTOUT
- integer, save :: NTOUT,NSNAP
-
- logical, save :: SIMULATION_TYPE_DYN = .false.
-
- logical, save :: TPV16 = .false.
-
- logical, save :: RATE_AND_STATE = .false.
-
- real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
-
- public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- SIMULATION_TYPE_DYN,RATE_AND_STATE
-
-
-contains
-
-!=====================================================================
-! BC_DYNFLT_init initializes dynamic faults
-!
-! prname fault database is read from file prname_fault_db.bin
-! Minv inverse mass matrix
-! dt global time step
-!
-subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,vel)
-
- character(len=256), intent(in) :: prname ! 'proc***'
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- double precision, intent(in) :: DTglobal
- integer, intent(in) :: nt
- real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
-
- real(kind=CUSTOM_REAL) :: dt
- integer :: iflt,ier,dummy_idfault
- integer :: nbfaults
- integer :: size_Kelvin_Voigt
- integer :: SIMULATION_TYPE
- character(len=256) :: filename
- integer, parameter :: IIN_PAR =151
- integer, parameter :: IIN_BIN =170
-
- NAMELIST / BEGIN_FAULT / dummy_idfault
-
- dummy_idfault = 0
-
- open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File Par_file_faults not found: assume no faults'
- close(IIN_PAR)
- return
- endif
-
- read(IIN_PAR,*) nbfaults
- if (nbfaults==0) return
-
- filename = prname(1:len_trim(prname))//'fault_db.bin'
- open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File ',trim(filename),' not found. Abort'
- stop
- endif
- ! WARNING TO DO: should be an MPI abort
-
- ! Reading etas of each fault
- do iflt = 1,nbfaults
- read(IIN_PAR,*) ! etas
- enddo
- read(IIN_PAR,*) SIMULATION_TYPE
- if ( SIMULATION_TYPE /= 1 ) then
- close(IIN_BIN)
- close(IIN_PAR)
- return
- endif
- SIMULATION_TYPE_DYN = .true.
- read(IIN_PAR,*) NTOUT
- read(IIN_PAR,*) NSNAP
- read(IIN_PAR,*) V_HEALING
- read(IIN_PAR,*) V_RUPT
-
- read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
- allocate( faults(nbfaults) )
- dt = real(DTglobal)
- do iflt=1,nbfaults
- read(IIN_PAR,nml=BEGIN_FAULT,end=100)
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,vel,iflt)
- enddo
- close(IIN_BIN)
- close(IIN_PAR)
-
- filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
- open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(6,*) 'File ',trim(filename),' not found. Abort'
- stop
- endif
- read(IIN_BIN) size_Kelvin_Voigt
- if (size_Kelvin_Voigt > 0) then
- allocate(Kelvin_Voigt_eta(size_Kelvin_Voigt))
- read(IIN_BIN) Kelvin_Voigt_eta
- endif
- close(IIN_BIN)
- return
-100 stop 'Did not find BEGIN_FAULT block #'
- ! WARNING TO DO: should be an MPI abort
-
-end subroutine BC_DYNFLT_init
-
-!---------------------------------------------------------------------
-
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,vel,iflt)
-
- real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
-
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: dt
-
- real(kind=CUSTOM_REAL) :: S1,S2,S3
- integer :: n1,n2,n3
-
- NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
-
- call initialize_fault(bc,IIN_BIN,Minv,dt)
-
- if (bc%nspec>0) then
-
- allocate(bc%T(3,bc%nglob))
- allocate(bc%D(3,bc%nglob))
- allocate(bc%V(3,bc%nglob))
- bc%T = 0e0_CUSTOM_REAL
- bc%D = 0e0_CUSTOM_REAL
- bc%V = 0e0_CUSTOM_REAL
-
- ! Set initial fault stresses
- allocate(bc%T0(3,bc%nglob))
- S1 = 0e0_CUSTOM_REAL
- S2 = 0e0_CUSTOM_REAL
- S3 = 0e0_CUSTOM_REAL
- n1=0
- n2=0
- n3=0
- read(IIN_PAR, nml=INIT_STRESS)
- bc%T0(1,:) = S1
- bc%T0(2,:) = S2
- bc%T0(3,:) = S3
- call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
- call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
- call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
-
- !WARNING : Quick and dirty free surface condition at z=0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
- ! end do
-
- ! Set friction parameters and initialize friction variables
- allocate(bc%MU(bc%nglob))
- if (RATE_AND_STATE) then
- allocate(bc%rsf)
- call rsf_init(bc%rsf,bc%T0,bc%V,bc%Fload,bc%coord,IIN_PAR)
- else
- allocate(bc%swf)
- call swf_init(bc%swf,bc%MU,bc%coord,IIN_PAR)
- if (TPV16) call TPV16_init() !WARNING: ad hoc, initializes T0 and swf
- endif
-
- endif
-
-! call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc)
-
-!--------------------------------------------------------
-contains
-
-subroutine TPV16_init
-
- integer :: ier, ipar
- integer, parameter :: IIN_NUC =270 ! WARNING: not safe, should look for an available unit
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
- dyn_fc,swcd,cohes,tim_forcedRup
- integer, dimension(bc%nglob) :: inp_nx,inp_nz
- real(kind=CUSTOM_REAL) :: minX, siz_str,siz_dip, hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- integer :: relz_num,sub_relz_num, num_cell_str,num_cell_dip, hypo_cell_str,hypo_cell_dip
- integer :: i
-
- open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
- read(IIN_NUC,*) relz_num,sub_relz_num
- read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
- read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- do ipar=1,bc%nglob
- read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
- Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
- enddo
- close(IIN_NUC)
-
- minX = minval(bc%coord(1,:))
-
- do i=1,bc%nglob
-
- ! WARNING: nearest neighbor interpolation
- ipar = minloc( (minX+loc_str(:)-bc%coord(1,i))**2 + (-loc_dip(:)-bc%coord(3,i))**2 , 1)
- !loc_dip is negative of Z-coord
-
- bc%T0(3,i) = -sigma0(ipar)
- bc%T0(1,i) = tau0_str(ipar)
- bc%T0(2,i) = tau0_dip(ipar)
-
- bc%swf%mus(i) = static_fc(ipar)
- bc%swf%mud(i) = dyn_fc(ipar)
- bc%swf%Dc(i) = swcd(ipar)
- bc%swf%C(i) = cohes(ipar)
- bc%swf%T(i) = tim_forcedRup(ipar)
-
- enddo
-
-end subroutine TPV16_init
-
-end subroutine init_one_fault
-
-!---------------------------------------------------------------------
-! REPLACES the value of a fault parameter inside an area with prescribed shape
-subroutine init_2d_distribution(a,coord,iin,n)
-!JPA refactor: background value should be an argument
-
- real(kind=CUSTOM_REAL), intent(inout) :: a(:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- integer, intent(in) :: iin,n
-
- real(kind=CUSTOM_REAL) :: b(size(a))
- character(len=20) :: shape
- real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
- real(kind=CUSTOM_REAL) :: r1(size(a))
- integer :: i
-
- NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
-
- if (n==0) return
-
- do i=1,n
- shape = ''
- xc = 0e0_CUSTOM_REAL
- yc = 0e0_CUSTOM_REAL
- zc = 0e0_CUSTOM_REAL
- r = 0e0_CUSTOM_REAL
- l = 0e0_CUSTOM_REAL
- lx = 0e0_CUSTOM_REAL
- ly = 0e0_CUSTOM_REAL
- lz = 0e0_CUSTOM_REAL
- valh = 0e0_CUSTOM_REAL
-
- read(iin,DIST2D)
- select case(shape)
- case ('circle')
- b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
- case ('circle-exp')
- r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
- where(r1<r)
- b =exp(r1**2/(r1**2 - r**2) )
- elsewhere
- b =0._CUSTOM_REAL
- endwhere
- case ('ellipse')
- b = heaviside( 1e0_CUSTOM_REAL - sqrt( (coord(1,:)-xc)**2/lx**2 + (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2 ) )
- case ('square')
- b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('cylinder')
- b = heaviside(r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2)) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('rectangle')
- b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('rectangle-taper')
- b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
- case default
- stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
- end select
-
- ! REPLACE the value inside the prescribed area
- where (b /= 0e0_CUSTOM_REAL) a = b
- enddo
-
-end subroutine init_2d_distribution
-
-!---------------------------------------------------------------------
-elemental function heaviside(x)
-
- real(kind=CUSTOM_REAL), intent(in) :: x
- real(kind=CUSTOM_REAL) :: heaviside
-
- if (x>=0e0_CUSTOM_REAL) then
- heaviside = 1e0_CUSTOM_REAL
- else
- heaviside = 0e0_CUSTOM_REAL
- endif
-
-end function heaviside
-
-!=====================================================================
-! adds boundary term Bt into Force array for each fault.
-!
-! NOTE: On non-split nodes at fault edges, dD=dV=dA=0
-! and the net contribution of B*T is =0
-!
-subroutine bc_dynflt_set3d_all(F,V,D)
-
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V,D
- real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
-
- integer :: i
-
- if (.not. allocated(faults)) return
- do i=1,size(faults)
- call BC_DYNFLT_set3d(faults(i),F,V,D,i)
- enddo
-
-end subroutine bc_dynflt_set3d_all
-
-!---------------------------------------------------------------------
-subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
-
- use specfem_par, only: it,NSTEP,myrank
-
- real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
- integer, intent(in) :: iflt
-
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T,dD,dV,dA
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength,tStick,tnew, &
- theta_old, theta_new, dc, &
- ta,Vf_old,Vf_new,TxExt
- real(kind=CUSTOM_REAL) :: half_dt,TLoad,DTau0,GLoad,time
- integer :: i
-
- if (bc%nspec > 0) then !Surendra : for parallel faults
-
- half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- Vf_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
-
- ! get predicted values
- dD = get_jump(bc,D) ! dD_predictor
- dV = get_jump(bc,V) ! dV_predictor
- dA = get_weighted_jump(bc,MxA) ! dA_free
-
- ! rotate to fault frame (tangent,normal)
- ! component 3 is normal to the fault
- dD = rotate(bc,dD,1)
- dV = rotate(bc,dV,1)
- dA = rotate(bc,dA,1)
-
- ! T_stick
- T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
- T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
- T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
-
- !Warning : dirty particular free surface condition z = 0.
- ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
- ! end do
-
- ! add initial stress
- T = T + bc%T0
-
- ! Solve for normal stress (negative is compressive)
- ! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
-
- ! smooth loading within nucleation patch
- !WARNING : ad hoc for SCEC benchmark TPV10x
- if (RATE_AND_STATE) then
- TxExt = 0._CUSTOM_REAL
- TLoad = 1.0_CUSTOM_REAL
- DTau0 = 25e6_CUSTOM_REAL
- time = it*bc%dt !time will never be zero. it starts from 1
- if (time <= TLoad) then
- GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
- else
- GLoad = 1.0_CUSTOM_REAL
- endif
- TxExt = DTau0 * bc%Fload * GLoad
- T(1,:) = T(1,:) + TxExt
- endif
-
- tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
-
- if (.not. RATE_AND_STATE) then ! Update slip weakening friction:
- ! Update slip state variable
- ! WARNING: during opening the friction state variable should not evolve
- theta_old = bc%swf%theta
- call swf_update_state(bc%D,dD,bc%V,bc%swf)
-
- ! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
-
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
- if (TPV16) then
- where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
- endif
-
- ! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
-
- ! Solve for shear stress
- tnew = min(tStick,strength)
-
- else ! Update rate and state friction:
- !JPA the solver below can be refactored into a loop with two passes
-
- ! first pass
- theta_old = bc%rsf%theta
- call rsf_update_state(Vf_old,bc%dt,bc%rsf)
- do i=1,bc%nglob
- Vf_new(i)=rtsafe(funcd,0.0,Vf_old(i)+5.0,1e-5,tStick(i),-T(3,i),bc%Z(i),bc%rsf%f0(i), &
- bc%rsf%V0(i),bc%rsf%a(i),bc%rsf%b(i),bc%rsf%L(i),bc%rsf%theta(i),bc%rsf%StateLaw)
- enddo
-
- ! second pass
- bc%rsf%theta = theta_old
- call rsf_update_state(0.5e0_CUSTOM_REAL*(Vf_old + Vf_new),bc%dt,bc%rsf)
- do i=1,bc%nglob
- Vf_new(i)=rtsafe(funcd,0.0,Vf_old(i)+5.0,1e-5,tStick(i),-T(3,i),bc%Z(i),bc%rsf%f0(i), &
- bc%rsf%V0(i),bc%rsf%a(i),bc%rsf%b(i),bc%rsf%L(i),bc%rsf%theta(i),bc%rsf%StateLaw)
- enddo
-
- tnew = tStick - bc%Z*Vf_new
-
- endif
-
- tStick = max(tStick,1e0_CUSTOM_REAL) ! to avoid division by zero
- T(1,:) = tnew * T(1,:)/tStick
- T(2,:) = tnew * T(2,:)/tStick
-
- ! Save total tractions
- bc%T = T
-
- ! Subtract initial stress
- T = T - bc%T0
-
- if (RATE_AND_STATE) T(1,:) = T(1,:) - TxExt
- !JPA: this eliminates the effect of TxExt on the equations of motion. Why is it needed?
-
- ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
- dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
- dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
- dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
-
- ! Update slip and slip rate, in fault frame
- bc%D = dD
- bc%V = dV + half_dt*dA
-
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
-
- ! Add boundary term B*T to M*a
- call add_BT(bc,MxA,T)
-
- !-- intermediate storage of outputs --
- Vf_new = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- if(.not. RATE_AND_STATE) then
- theta_new = bc%swf%theta
- dc = bc%swf%dc
- else
- theta_new = bc%rsf%theta
- dc = bc%rsf%L
- endif
-
- call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
- Vf_old, Vf_new, it*bc%dt,bc%dt)
-! call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
-
- !-- outputs --
- ! write dataT every NTOUT time step or at the end of simulation
-! if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it,bc%rsf%StateLaw)
-
- endif
-
- ! write dataXZ every NSNAP time step
- if ( mod(it,NSNAP) == 0) then
- if (.NOT. PARALLEL_FAULT) then
- if (bc%nspec > 0) call write_dataXZ(bc%dataXZ,it,iflt)
- else
- call gather_dataXZ(bc)
- if (myrank==0) call write_dataXZ(bc%dataXZ_all,it,iflt)
- endif
- endif
-
- if ( it == NSTEP) then
- if (.NOT. PARALLEL_FAULT) then
- call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
- else
- if (myrank==0) call SCEC_Write_RuptureTime(bc%dataXZ_all,bc%dt,NSTEP,iflt)
- endif
- endif
-
-end subroutine BC_DYNFLT_set3d
-
-!===============================================================
-
-subroutine swf_init(f,mu,coord,IIN_PAR)
-
- type(swf_type), intent(out) :: f
- real(kind=CUSTOM_REAL), intent(out) :: mu(:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- integer, intent(in) :: IIN_PAR
-
- integer :: nglob
- real(kind=CUSTOM_REAL) :: mus,mud,dc,C,T
- integer :: nmus,nmud,ndc,nC,nForcedRup
-
- NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
-
- nglob = size(mu)
- allocate( f%mus(nglob) )
- allocate( f%mud(nglob) )
- allocate( f%Dc(nglob) )
- allocate( f%theta(nglob) )
- allocate( f%C(nglob) )
- allocate( f%T(nglob) )
-
- ! WARNING: if V_HEALING is negative we turn off healing
- f%healing = (V_HEALING > 0e0_CUSTOM_REAL)
-
- mus = 0.6e0_CUSTOM_REAL
- mud = 0.1e0_CUSTOM_REAL
- dc = 1e0_CUSTOM_REAL
- C = 0._CUSTOM_REAL
- T = HUGEVAL
-
- nmus = 0
- nmud = 0
- ndc = 0
- nC = 0
- nForcedRup = 0
-
- read(IIN_PAR, nml=SWF)
-
- f%mus = mus
- f%mud = mud
- f%Dc = dc
- f%C = C
- f%T = T
- call init_2d_distribution(f%mus,coord,IIN_PAR,nmus)
- call init_2d_distribution(f%mud,coord,IIN_PAR,nmud)
- call init_2d_distribution(f%Dc ,coord,IIN_PAR,ndc)
- call init_2d_distribution(f%C ,coord,IIN_PAR,nC)
- call init_2d_distribution(f%T ,coord,IIN_PAR,nForcedRup)
-
- f%theta = 0e0_CUSTOM_REAL
-
- mu = swf_mu(f)
-
-end subroutine swf_init
-
-!---------------------------------------------------------------------
-subroutine swf_update_state(dold,dnew,vold,f)
-
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
- type(swf_type), intent(inout) :: f
-
- real(kind=CUSTOM_REAL) :: vnorm
- integer :: k,npoin
-
- f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
-
- if (f%healing) then
- npoin = size(vold,2)
- do k=1,npoin
- vnorm = sqrt(vold(1,k)**2 + vold(2,k)**2)
- if (vnorm<V_HEALING) f%theta(k) = 0e0_CUSTOM_REAL
- enddo
- endif
-end subroutine swf_update_state
-
-!---------------------------------------------------------------------
-function swf_mu(f) result(mu)
-
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
-
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- mu = max( mu, f%mud)
-
-end function swf_mu
-
-!=====================================================================
-
-subroutine rsf_init(f,T0,V,nucFload,coord,IIN_PAR)
-
- type(rsf_type), intent(out) :: f
- real(kind=CUSTOM_REAL), intent(in) :: T0(:,:)
- real(kind=CUSTOM_REAL), intent(inout) :: V(:,:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- real(kind=CUSTOM_REAL), pointer :: nucFload(:)
- integer, intent(in) :: IIN_PAR
-
- real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init,fw,Vw, C,T
- integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init,nfw,nVw, nC,nForcedRup
- real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
- real(kind=CUSTOM_REAL) :: x,z
- logical :: c1,c2,c3,c4
- real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
- integer :: i !,nglob_bulk
- real(kind=CUSTOM_REAL) :: Fload
- integer :: nFload
-! real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: init_vel
- integer :: nglob
-
- NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup,Vw,fw,nVw,nfw
- NAMELIST / ASP / Fload,nFload
-
- nglob = size(coord,2)
-
- allocate( f%V0(nglob) )
- allocate( f%f0(nglob) )
- allocate( f%a(nglob) )
- allocate( f%b(nglob) )
- allocate( f%L(nglob) )
- allocate( f%V_init(nglob) )
- allocate( f%theta(nglob) )
- allocate( f%C(nglob) )
- allocate( f%T(nglob) )
- allocate( f%fw(nglob) )
- allocate( f%Vw(nglob) )
-
- V0 =1.e-6_CUSTOM_REAL
- f0 =0.6_CUSTOM_REAL
- a =0.0080_CUSTOM_REAL !0.0080_CUSTOM_REAL
- b =0.0040_CUSTOM_REAL !0.0120_CUSTOM_REAL
- L =0.0135_CUSTOM_REAL
- V_init =1.e-12_CUSTOM_REAL
- theta_init =1.084207680000000e+09_CUSTOM_REAL
- C = 0._CUSTOM_REAL
- T = HUGEVAL
- fw = 0.2_CUSTOM_REAL
- Vw = 0.1_CUSTOM_REAL
-
- nV0 =0
- nf0 =0
- na =0
- nb =0
- nL =0
- nV_init =0
- ntheta_init =0
- nC = 0
- nForcedRup = 0
- nfw =0
- nVw =0
-
- read(IIN_PAR, nml=RSF)
-
- f%V0 = V0
- f%f0 = f0
- f%a = a
- f%b = b
- f%L = L
- f%V_init = V_init
- f%theta = theta_init
- f%C = C
- f%T = T
- f%fw = fw
- f%Vw = Vw
-
- call init_2d_distribution(f%V0,coord,IIN_PAR,nV0)
- call init_2d_distribution(f%f0,coord,IIN_PAR,nf0)
- call init_2d_distribution(f%a,coord,IIN_PAR,na)
- call init_2d_distribution(f%b,coord,IIN_PAR,nb)
- call init_2d_distribution(f%L,coord,IIN_PAR,nL)
- call init_2d_distribution(f%V_init,coord,IIN_PAR,nV_init)
- call init_2d_distribution(f%theta,coord,IIN_PAR,ntheta_init)
- call init_2d_distribution(f%C,coord,IIN_PAR,nC)
- call init_2d_distribution(f%T,coord,IIN_PAR,nForcedRup)
- call init_2d_distribution(f%fw,coord,IIN_PAR,nfw)
- call init_2d_distribution(f%Vw,coord,IIN_PAR,nVw)
-
-!!$ ! WARNING : Not general enough
-!!$ nglob_bulk = size(vel,2)
-!!$ allocate(init_vel(3,nglob_bulk))
-!!$ init_vel = 0._CUSTOM_REAL
-!!$ init_vel(1,bc%ibulk1) = -f%V_init/2._CUSTOM_REAL
-!!$ init_vel(1,bc%ibulk2) = f%V_init/2._CUSTOM_REAL
-!!$ where(ystore > 0) init_vel(1,:) = -V_init/2._CUSTOM_REAL
-!!$ where(ystore < 0) init_vel(1,:) = V_init/2._CUSTOM_REAL
-!!$ !init_vel = rotate(bc,init_vel,-1) ! directly assigned in global coordinates here as it is a simplified case
-!!$ vel = vel + init_vel
-
- ! WARNING: below is an ad-hoc setting of a(x,z) for some SCEC benchmark
- ! This should be instead an option in init_2d_distribution
- W1=15000._CUSTOM_REAL
- W2=7500._CUSTOM_REAL
- w=3000._CUSTOM_REAL
- hypo_z = -7500._CUSTOM_REAL
- do i=1,nglob
- x=coord(1,i)
- z=coord(3,i)
- c1=abs(x)<W1+w
- c2=abs(x)>W1
- c3=abs(z-hypo_z)<W2+w
- c4=abs(z-hypo_z)>W2
- if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
-
- if (c1 .and. c2) then
- b11 = w/(abs(x)-W1-w)
- b12 = w/(abs(x)-W1)
- B1 = HALF * (ONE + tanh(b11 + b12))
- elseif(abs(x)<=W1) then
- B1 = 1._CUSTOM_REAL
- else
- B1 = 0._CUSTOM_REAL
- endif
-
- if (c3 .and. c4) then
- b21 = w/(abs(z-hypo_z)-W2-w)
- b22 = w/(abs(z-hypo_z)-W2)
- B2 = HALF * (ONE + tanh(b21 + b22))
- elseif(abs(z-hypo_z)<=W2) then
- B2 = 1._CUSTOM_REAL
- else
- B2 = 0._CUSTOM_REAL
- endif
-
- f%a(i) = 0.008 + 0.008 * (ONE - B1*B2)
- f%Vw(i) = 0.1 + 0.9 * (ONE - B1*B2)
-
- elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
- f%a(i) = 0.008
- f%Vw(i) = 0.1_CUSTOM_REAL
- else
- f%a(i) = 0.016
- f%Vw(i) = 1.0_CUSTOM_REAL
- endif
-
- enddo
-
- ! WARNING: The line below scratches an earlier initialization of theta through theta_init
- ! We should implement it as an option for the user
- if(f%stateLaw == 1) then
- f%theta = f%L/f%V0 &
- * exp( ( f%a * log(TWO*sinh(-sqrt(T0(1,:)**2+T0(2,:)**2)/T0(3,:)/f%a)) &
- - f%f0 - f%a*log(f%V_init/f%V0) ) &
- / f%b )
- else
- f%theta = f%a * log(TWO*f%V0/f%V_init * sinh(-sqrt(T0(1,:)**2+T0(2,:)**2)/T0(3,:)/f%a))
- endif
-
- ! WARNING : ad hoc for SCEC benchmark TPV10x
- allocate( nucFload(nglob) )
- Fload = 0.e0_CUSTOM_REAL
- nFload = 0
- read(IIN_PAR, nml=ASP)
- nucFload = Fload
- call init_2d_distribution(nucFload,coord,IIN_PAR,nFload)
-
- ! WARNING: the line below is only valid for pure strike-slip faulting
- V(1,:) = f%V_init
-
-end subroutine rsf_init
-
-!---------------------------------------------------------------------
-! Rate and state friction coefficient
-function rsf_mu(f,V) result(mu)
-
- type(rsf_type), intent(in) :: f
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- real(kind=CUSTOM_REAL) :: mu(size(V))
- mu = f%a * asinh( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
-
-end function rsf_mu
-
-!---------------------------------------------------------------------
-subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
-
- real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
- double precision :: arg,fn,df,x
- integer :: statelaw
-
- if(statelaw == 1) then
- arg = exp((f0+dble(b)*log(V0*theta/L))/a)/TWO/V0
- else
- arg = exp(theta/a)/TWO/V0
- endif
- fn = tStick - Z*x - a*Seff*asinh(x*arg)
- df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
-end subroutine funcd
-
-!---------------------------------------------------------------------
-function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
-
- integer, parameter :: MAXIT=200
- real(kind=CUSTOM_REAL) :: x1,x2,xacc
- EXTERNAL funcd
- integer :: j
- !real(kind=CUSTOM_REAL) :: df,dx,dxold,f,fh,fl,temp,xh,xl
- double precision :: df,dx,dxold,f,fh,fl,temp,xh,xl,rtsafe
- real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
- integer :: statelaw
-
- call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- if( (fl>0 .and. fh>0) .or. (fl<0 .and. fh<0) ) stop 'root must be bracketed in rtsafe'
- if(fl==0.) then
- rtsafe=x2
- return
- elseif(fh==0.) then
- rtsafe=x2
- return
- elseif(fl<0) then
- xl=x1
- xh=x2
- else
- xh=x1
- xl=x2
- endif
-
- rtsafe=0.5d0*(x1+x2)
- dxold=abs(x2-x1)
- dx=dxold
- call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- do j=1,MAXIT
- if( ((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0 .or. abs(2.*f)>abs(dxold*df) ) then
- dxold=dx
- dx=0.5d0*(xh-xl)
- rtsafe=xl+dx
- if(xl==rtsafe) return
- else
- dxold=dx
- dx=f/df
- temp=rtsafe
- rtsafe=rtsafe-dx
- if(temp==rtsafe) return
- endif
- if(abs(dx)<xacc) return
- call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- if(f<0.) then
- xl=rtsafe
- else
- xh=rtsafe
- endif
- enddo
- stop 'rtsafe exceeding maximum iterations'
- return
-
-end function rtsafe
-
-!---------------------------------------------------------------------
-subroutine rsf_update_state(V,dt,f)
-
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- type(rsf_type), intent(inout) :: f
- real(kind=CUSTOM_REAL), intent(in) :: dt
-
- real(kind=CUSTOM_REAL) :: vDtL(size(V))
- real(kind=CUSTOM_REAL) :: f_ss(size(V)),xi_ss(size(V)),fLV(size(V))
-
- vDtL = V * dt / f%L
-
- ! ageing law
- if (f%StateLaw == 1) then
- where(vDtL > 1.e-5_CUSTOM_REAL)
- f%theta = f%theta * exp(-vDtL) + f%L/V * (ONE - exp(-vDtL))
- elsewhere
- f%theta = f%theta * exp(-vDtL) + dt*( ONE - HALF*vDtL )
- endwhere
-
- ! slip law : by default use strong rate-weakening
- else
-! f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
- where(V /= 0._CUSTOM_REAL)
- fLV = f%f0 - (f%b - f%a)*log(V/f%V0)
- f_ss = f%fw + (fLV - f%fw)/(ONE + (V/f%Vw)**8)**0.125
- xi_ss = f%a * log( TWO*f%V0/V * sinh(f_ss/f%a) )
- f%theta = xi_ss + (f%theta - xi_ss) * exp(-vDtL)
- elsewhere
- f%theta = f%theta
- endwhere
- endif
-
-end subroutine rsf_update_state
-
-
-!===============================================================
-! OUTPUTS
-subroutine init_dataT(DataT,coord,nglob,NT,iflt)
-
- use specfem_par, only : NPROC,myrank
- ! NT = total number of time steps
-
- integer, intent(in) :: nglob,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
- type (dataT_type), intent(out) :: DataT
-
- real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
- integer :: i, iglob , IIN, ier, jflt, np, k
- character(len=70) :: tmpname
-
- integer :: ipoin, ipoin_local, iproc
-
- ! 1. read fault output coordinates from user file,
- ! 2. define iglob: the fault global index of the node nearest to user
- ! requested coordinate
-
- IIN = 251
- open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
- read(IIN,*) np
- DataT%npoin =0
- do i=1,np
- read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
- if (jflt==iflt) DataT%npoin = DataT%npoin +1
- enddo
- close(IIN)
-
- if (DataT%npoin == 0) return
-
- allocate(DataT%iglob(DataT%npoin))
- allocate(DataT%name(DataT%npoin))
- allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
-
- open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
- if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
- read(IIN,*) np
- k = 0
- do i=1,np
- read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
- if (jflt/=iflt) cycle
- k = k+1
- DataT%name(k) = tmpname
- !search nearest node
- distkeep = huge(distkeep)
-
- do iglob=1,nglob
- dist = sqrt((coord(1,iglob)-xtarget)**2 &
- + (coord(2,iglob)-ytarget)**2 &
- + (coord(3,iglob)-ztarget)**2)
- if (dist < distkeep) then
- distkeep = dist
- DataT%iglob(k) = iglob
- endif
- enddo
- DataT%dist(k) = distkeep !Surendra : for parallel fault
- enddo
-
- !Surendra : for parallel fault
- if (PARALLEL_FAULT) then
- allocate(DataT%islice(DataT%npoin))
- allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
- allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
- call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
- call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
- if(myrank==0) then
- do ipoin = 1,DataT%npoin
- iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
- DataT%islice(ipoin) = iproc
- DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
- enddo
- endif
-
- call bcast_all_i(DataT%islice,DataT%npoin)
- call bcast_all_i(DataT%iglob,DataT%npoin)
-
- DataT%npoin_local = 0
- do ipoin = 1,DataT%npoin
- if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
- enddo
- allocate(DataT%glob_indx(DataT%npoin_local))
- do ipoin = 1,DataT%npoin
- if(myrank == DataT%islice(ipoin)) then
- ipoin_local = ipoin_local + 1
- DataT%glob_indx(ipoin_local) = ipoin
- endif
- enddo
- else
- DataT%npoin_local = DataT%npoin
- endif !Parallel_fault
-
-
- ! 3. allocate arrays and set to zero
- allocate(DataT%d1(NT,DataT%npoin))
- allocate(DataT%v1(NT,DataT%npoin))
- allocate(DataT%t1(NT,DataT%npoin))
- allocate(DataT%d2(NT,DataT%npoin))
- allocate(DataT%v2(NT,DataT%npoin))
- allocate(DataT%t2(NT,DataT%npoin))
- allocate(DataT%t3(NT,DataT%npoin))
- allocate(DataT%theta(NT,DataT%npoin))
- DataT%d1 = 0e0_CUSTOM_REAL
- DataT%v1 = 0e0_CUSTOM_REAL
- DataT%t1 = 0e0_CUSTOM_REAL
- DataT%d2 = 0e0_CUSTOM_REAL
- DataT%v2 = 0e0_CUSTOM_REAL
- DataT%t2 = 0e0_CUSTOM_REAL
- DataT%t3 = 0e0_CUSTOM_REAL
- DataT%theta = 0e0_CUSTOM_REAL
-
- close(IIN)
-
-end subroutine init_dataT
-
-!---------------------------------------------------------------
-subroutine store_dataT(dataT,d,v,t,theta,itime)
-
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
- integer, intent(in) :: itime
-
- integer :: i,k,ipoin
-
- do ipoin=1,dataT%npoin_local
- if (PARALLEL_FAULT) then
- i = DataT%glob_indx(ipoin)
- else
- i = ipoin
- endif
- k = dataT%iglob(i)
- dataT%d1(itime,i) = d(1,k)
- dataT%d2(itime,i) = d(2,k)
- dataT%v1(itime,i) = v(1,k)
- dataT%v2(itime,i) = v(2,k)
- dataT%t1(itime,i) = t(1,k)
- dataT%t2(itime,i) = t(2,k)
- dataT%t3(itime,i) = t(3,k)
- dataT%theta(itime,i) = theta(k)
- enddo
-
-end subroutine store_dataT
-
-!-----------------------------------------------------------------
-subroutine write_dataT_all(nt,statelaw)
-
- integer, intent(in) :: nt
- integer, intent(in) :: statelaw
-
- integer :: i
-
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
- enddo
-
-end subroutine write_dataT_all
-
-!------------------------------------------------------------------------
-subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
-
- type(dataT_type), intent(in) :: dataT
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT
- integer, intent(in) :: statelaw
-
- integer :: i,k,IOUT,ipoin
- character :: NTchar*5
- integer :: today(3), now(3)
-
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
-
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- write(NTchar,1) NT
- NTchar = adjustl(NTchar)
-
-1 format(I5)
- do ipoin=1,dataT%npoin_local
- if (PARALLEL_FAULT) then
- i = DataT%glob_indx(ipoin)
- else
- i = ipoin
- endif
-
- open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV104"
- write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) today(2), today(1), today(3), now
- write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
- write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
- write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# location=",trim(dataT%name(i))
- write(IOUT,*) "# Column #1 = Time (s)"
- write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
- write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
- write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
- write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
- write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
- write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
- write(IOUT,*) "# Column #8 = normal stress (MPa)"
- if (RATE_AND_STATE) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
- write(IOUT,*) "#"
- write(IOUT,*) "# The line below lists the names of the data fields:"
- if (.not. RATE_AND_STATE) then
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
- -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
- dataT%t3(k,i)/1.0e6_CUSTOM_REAL
- enddo
- else
- if(statelaw == 1) then
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
- -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
- dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
- enddo
- else
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress psi"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
- -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
- dataT%t3(k,i)/1.0e6_CUSTOM_REAL, dataT%theta(k,i)
- enddo
- endif
- endif
- close(IOUT)
-
-
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
-
-
- enddo
-
-end subroutine SCEC_write_dataT
-
-!-------------------------------------------------------------------------------------------------
-
-subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
-
- type(dataXZ_type), intent(in) :: dataXZ
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT,iflt
-
- integer :: i,IOUT
- character(len=70) :: filename
- integer*4 today(3), now(3)
-
- call idate(today) ! today(1)=day, (2)=month, (3)=year
- call itime(now) ! now(1)=hour, (2)=minute, (3)=second
-
-
- write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
-
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- open(IOUT,file=trim(filename),status='replace')
- write(IOUT,*) "# problem=TPV104"
- write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) today(2), today(1), today(3), now
- write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
- write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
- write(IOUT,*) "# Column #1 = horizontal coordinate, distance along strike (m)"
- write(IOUT,*) "# Column #2 = vertical coordinate, distance down-dip (m)"
- write(IOUT,*) "# Column #3 = rupture time (s)"
- write(IOUT,*) "# "
- write(IOUT,*) "j k t"
- do i = 1,size(dataXZ%tRUP)
- write(IOUT,'(3(E15.7))') dataXZ%xcoord(i), -dataXZ%zcoord(i), dataXZ%tRUP(i)
- end do
-
- close(IOUT)
-
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
-
-end subroutine SCEC_Write_RuptureTime
-
-!-------------------------------------------------------------------------------------------------
-subroutine init_dataXZ(DataXZ,bc)
-
- use specfem_par, only : NPROC,myrank
-
- type(dataXZ_type), intent(inout) :: DataXZ
- type(bc_dynflt_type) :: bc
-
- integer :: npoin_all,iproc
-
- DataXZ%npoin = bc%nglob
-
- if(bc%nglob > 0) then
-
- allocate(DataXZ%stg(bc%nglob))
- if(.not. RATE_AND_STATE) then
- DataXZ%sta => bc%swf%theta
- else
- DataXZ%sta => bc%rsf%theta
- endif
- DataXZ%d1 => bc%d(1,:)
- DataXZ%d2 => bc%d(2,:)
- DataXZ%v1 => bc%v(1,:)
- DataXZ%v2 => bc%v(2,:)
- DataXZ%t1 => bc%t(1,:)
- DataXZ%t2 => bc%t(2,:)
- DataXZ%t3 => bc%t(3,:)
- DataXZ%xcoord => bc%coord(1,:)
- DataXZ%ycoord => bc%coord(2,:)
- DataXZ%zcoord => bc%coord(3,:)
- allocate(DataXZ%tRUP(bc%nglob))
- allocate(DataXZ%tPZ(bc%nglob))
-
- !Percy , setting up initial rupture time null for all faults.
- DataXZ%tRUP = 0e0_CUSTOM_REAL
- DataXZ%tPZ = 0e0_CUSTOM_REAL
-
- endif
-
- !Surendra : for parallel fault
- if (PARALLEL_FAULT) then
- npoin_all = 0
- call sum_all_i(bc%nglob,npoin_all)
- if (myrank==0 .and. npoin_all>0) then
- bc%DataXZ_all%npoin = npoin_all
- allocate(bc%DataXZ_all%xcoord(npoin_all))
- allocate(bc%DataXZ_all%ycoord(npoin_all))
- allocate(bc%DataXZ_all%zcoord(npoin_all))
- allocate(bc%DataXZ_all%t1(npoin_all))
- allocate(bc%DataXZ_all%t2(npoin_all))
- allocate(bc%DataXZ_all%t3(npoin_all))
- allocate(bc%DataXZ_all%d1(npoin_all))
- allocate(bc%DataXZ_all%d2(npoin_all))
- allocate(bc%DataXZ_all%v1(npoin_all))
- allocate(bc%DataXZ_all%v2(npoin_all))
- allocate(bc%DataXZ_all%tRUP(npoin_all))
- allocate(bc%DataXZ_all%tPZ(npoin_all))
- allocate(bc%DataXZ_all%stg(npoin_all))
- allocate(bc%DataXZ_all%sta(npoin_all))
- endif
-
- allocate(bc%npoin_perproc(NPROC))
- bc%npoin_perproc=0
- call gather_all_i(DataXZ%npoin,1,bc%npoin_perproc,1,NPROC)
-
- allocate(bc%poin_offset(NPROC))
- bc%poin_offset(1)=0
- do iproc=2,NPROC
- bc%poin_offset(iproc) = sum(bc%npoin_perproc(1:iproc-1))
- enddo
-
- call gatherv_all_cr(DataXZ%xcoord,DataXZ%npoin,bc%DataXZ_all%xcoord,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(DataXZ%ycoord,DataXZ%npoin,bc%DataXZ_all%ycoord,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(DataXZ%zcoord,DataXZ%npoin,bc%DataXZ_all%zcoord,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- endif
-
-end subroutine init_dataXZ
-!---------------------------------------------------------------
-
-subroutine gather_dataXZ(bc)
-
- use specfem_par, only : NPROC
-
- type(bc_dynflt_type), intent(inout) :: bc
-
- call gatherv_all_cr(bc%DataXZ%t1,bc%DataXZ%npoin,bc%DataXZ_all%t1,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%t2,bc%DataXZ%npoin,bc%DataXZ_all%t2,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%t3,bc%DataXZ%npoin,bc%DataXZ_all%t3,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%d1,bc%DataXZ%npoin,bc%DataXZ_all%d1,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%d2,bc%DataXZ%npoin,bc%DataXZ_all%d2,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%v1,bc%DataXZ%npoin,bc%DataXZ_all%v1,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%v2,bc%DataXZ%npoin,bc%DataXZ_all%v2,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%tRUP,bc%DataXZ%npoin,bc%DataXZ_all%tRUP,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%tPZ,bc%DataXZ%npoin,bc%DataXZ_all%tPZ,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%stg,bc%DataXZ%npoin,bc%DataXZ_all%stg,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
- call gatherv_all_cr(bc%DataXZ%sta,bc%DataXZ%npoin,bc%DataXZ_all%sta,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
-
-end subroutine gather_dataXZ
-
-!---------------------------------------------------------------
-subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
-
- type(dataXZ_type), intent(inout) :: dataXZ
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: stg,dold,dnew,dc,vold,vnew
- real(kind=CUSTOM_REAL), intent(in) :: time,dt
-
- integer :: i
-
- ! "stg" : strength .
-
- dataXZ%stg = stg
-
- do i = 1,size(stg)
- ! process zone time = first time when slip = dc (break down process).
- ! with linear time interpolation
- if (dataXZ%tPZ(i)==0e0_CUSTOM_REAL) then
- if (dold(i)<=dc(i) .and. dnew(i) >= dc(i)) then
- dataXZ%tPZ(i) = time-dt*(dnew(i)-dc(i))/(dnew(i)-dold(i))
- endif
- endif
- ! rupture time = first time when slip velocity = vc
- ! with linear time interpolation
- ! vc should be pre-defined as input data .
-
- if (dataXZ%tRUP(i)==0e0_CUSTOM_REAL) then
- if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
- endif
- enddo
-
- ! To do : add stress criteria (firs time strength is reached).
-
- ! note: the other arrays in dataXZ are pointers to arrays in bc
- ! they do not need to be updated here
-
-end subroutine store_dataXZ
-
-!---------------------------------------------------------------
-subroutine write_dataXZ(dataXZ,itime,iflt)
-
- type(dataXZ_type), intent(in) :: dataXZ
- integer, intent(in) :: itime,iflt
-
- character(len=70) :: filename
-
- write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
- ! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
- ! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
- ! compilers.
-
- ! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
-
- open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
-
- write(IOUT) dataXZ%xcoord
- write(IOUT) dataXZ%ycoord
- write(IOUT) dataXZ%zcoord
- write(IOUT) dataXZ%d1
- write(IOUT) dataXZ%d2
- write(IOUT) dataXZ%v1
- write(IOUT) dataXZ%v2
- write(IOUT) dataXZ%t1
- write(IOUT) dataXZ%t2
- write(IOUT) dataXZ%t3
- write(IOUT) dataXZ%sta
- write(IOUT) dataXZ%stg
- write(IOUT) dataXZ%tRUP
- write(IOUT) dataXZ%tPZ
- close(IOUT)
-
-end subroutine write_dataXZ
-
-
-end module fault_solver
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90 (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90 2012-10-27 03:15:37 UTC (rev 20950)
@@ -0,0 +1,1455 @@
+!=====================================================================
+!
+! s p e c f e m 3 d v e r s i o n 1 . 4
+! ---------------------------------------
+!
+! dimitri komatitsch and jeroen tromp
+! seismological laboratory - california institute of technology
+! (c) california institute of technology september 2006
+!
+! 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 2 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.
+!
+!===============================================================================================================
+
+! This module was written by:
+! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
+! Surendra Nadh Somala : Added Heterogenous initial stress capabilities (based on TPV16)
+! Surendra Nadh Somala : Added Rate and State Friction
+! Somala/Ampuero : fault parallelization
+
+module fault_solver_dynamic
+
+ use fault_solver_common
+
+ implicit none
+
+ include 'constants.h'
+
+ private
+
+ ! outputs on selected fault nodes at every time step:
+ ! slip, slip velocity, fault stresses
+ type dataT_type
+ integer :: npoin=0,npoin_local=0
+ integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
+ integer, dimension(:,:), pointer :: iglob_all=>null()
+ integer, dimension(:), pointer :: islice=>null(),glob_indx=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: dist=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: dist_all=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1=>null(),v1=>null(),t1=>null(), &
+ d2=>null(),v2=>null(),t2=>null(), &
+ t3=>null(),theta=>null()
+ character(len=70), dimension(:), pointer :: name=>null()
+ end type dataT_type
+
+
+ ! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
+ ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
+ ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
+ ! process zone time = first time when slip = Dc
+ type dataXZ_type
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: stg=>null(), sta=>null(), d1=>null(), d2=>null(), v1=>null(), v2=>null(), &
+ t1=>null(), t2=>null(), t3=>null(), tRUP=>null(), tPZ=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord=>null(), ycoord=>null(), zcoord=>null()
+ integer :: npoin=0
+ end type dataXZ_type
+
+ type swf_type
+ private
+ integer :: kind
+ logical :: healing = .false.
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), &
+ theta=>null(), T=>null(), C=>null()
+ end type swf_type
+
+ type rsf_type
+ private
+ integer :: StateLaw = 1 ! 1=ageing law, 2=slip law
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), &
+ V_init=>null(), &
+ a=>null(), b=>null(), theta=>null(), &
+ T=>null(), C=>null(), &
+ fw=>null(), Vw=>null()
+ end type rsf_type
+
+ type, extends (fault_type) :: bc_dynflt_type
+ private
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null(), Fload=>null()
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ,dataXZ_all
+ type(swf_type), pointer :: swf => null()
+ type(rsf_type), pointer :: rsf => null()
+ logical :: allow_opening = .false. ! default : do not allow opening
+ end type bc_dynflt_type
+
+ type(bc_dynflt_type), allocatable, save :: faults(:)
+
+ !slip velocity threshold for healing
+ !WARNING: not very robust
+ real(kind=CUSTOM_REAL), save :: V_HEALING
+
+ !slip velocity threshold for definition of rupture front
+ real(kind=CUSTOM_REAL), save :: V_RUPT
+
+ !Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
+
+ logical, save :: SIMULATION_TYPE_DYN = .false.
+
+ logical, save :: TPV16 = .false.
+
+ logical, save :: RATE_AND_STATE = .false.
+
+ real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
+
+ public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
+ SIMULATION_TYPE_DYN,RATE_AND_STATE
+
+
+contains
+
+!=====================================================================
+! BC_DYNFLT_init initializes dynamic faults
+!
+! prname fault database is read from file prname_fault_db.bin
+! Minv inverse mass matrix
+! dt global time step
+!
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,vel)
+
+ character(len=256), intent(in) :: prname ! 'proc***'
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ double precision, intent(in) :: DTglobal
+ integer, intent(in) :: nt
+ real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
+
+ real(kind=CUSTOM_REAL) :: dt
+ integer :: iflt,ier,dummy_idfault
+ integer :: nbfaults
+ integer :: size_Kelvin_Voigt
+ integer :: SIMULATION_TYPE
+ character(len=256) :: filename
+ integer, parameter :: IIN_PAR =151
+ integer, parameter :: IIN_BIN =170
+
+ NAMELIST / BEGIN_FAULT / dummy_idfault
+
+ dummy_idfault = 0
+
+ open(unit=IIN_PAR,file='DATA/FAULT/Par_file_faults',status='old',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File Par_file_faults not found: assume no faults'
+ close(IIN_PAR)
+ return
+ endif
+
+ read(IIN_PAR,*) nbfaults
+ if (nbfaults==0) return
+
+ filename = prname(1:len_trim(prname))//'fault_db.bin'
+ open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
+ ! WARNING TO DO: should be an MPI abort
+
+ ! Reading etas of each fault
+ do iflt = 1,nbfaults
+ read(IIN_PAR,*) ! etas
+ enddo
+ read(IIN_PAR,*) SIMULATION_TYPE
+ if ( SIMULATION_TYPE /= 1 ) then
+ close(IIN_BIN)
+ close(IIN_PAR)
+ return
+ endif
+ SIMULATION_TYPE_DYN = .true.
+ read(IIN_PAR,*) NTOUT
+ read(IIN_PAR,*) NSNAP
+ read(IIN_PAR,*) V_HEALING
+ read(IIN_PAR,*) V_RUPT
+
+ read(IIN_BIN) nbfaults ! should be the same as in IIN_PAR
+ allocate( faults(nbfaults) )
+ dt = real(DTglobal)
+ do iflt=1,nbfaults
+ read(IIN_PAR,nml=BEGIN_FAULT,end=100)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,vel,iflt)
+ enddo
+ close(IIN_BIN)
+ close(IIN_PAR)
+
+ filename = prname(1:len_trim(prname))//'Kelvin_voigt_eta.bin'
+ open(unit=IIN_BIN,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(6,*) 'File ',trim(filename),' not found. Abort'
+ stop
+ endif
+ read(IIN_BIN) size_Kelvin_Voigt
+ if (size_Kelvin_Voigt > 0) then
+ allocate(Kelvin_Voigt_eta(size_Kelvin_Voigt))
+ read(IIN_BIN) Kelvin_Voigt_eta
+ endif
+ close(IIN_BIN)
+ return
+100 stop 'Did not find BEGIN_FAULT block #'
+ ! WARNING TO DO: should be an MPI abort
+
+end subroutine BC_DYNFLT_init
+
+!---------------------------------------------------------------------
+
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,vel,iflt)
+
+ real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
+
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
+ integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+
+ real(kind=CUSTOM_REAL) :: S1,S2,S3
+ integer :: n1,n2,n3
+
+ NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
+
+ call initialize_fault(bc,IIN_BIN,Minv,dt)
+
+ if (bc%nspec>0) then
+
+ allocate(bc%T(3,bc%nglob))
+ allocate(bc%D(3,bc%nglob))
+ allocate(bc%V(3,bc%nglob))
+ bc%T = 0e0_CUSTOM_REAL
+ bc%D = 0e0_CUSTOM_REAL
+ bc%V = 0e0_CUSTOM_REAL
+
+ ! Set initial fault stresses
+ allocate(bc%T0(3,bc%nglob))
+ S1 = 0e0_CUSTOM_REAL
+ S2 = 0e0_CUSTOM_REAL
+ S3 = 0e0_CUSTOM_REAL
+ n1=0
+ n2=0
+ n3=0
+ read(IIN_PAR, nml=INIT_STRESS)
+ bc%T0(1,:) = S1
+ bc%T0(2,:) = S2
+ bc%T0(3,:) = S3
+ call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
+ call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
+ call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
+
+ !WARNING : Quick and dirty free surface condition at z=0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
+ ! end do
+
+ ! Set friction parameters and initialize friction variables
+ allocate(bc%MU(bc%nglob))
+ if (RATE_AND_STATE) then
+ allocate(bc%rsf)
+ call rsf_init(bc%rsf,bc%T0,bc%Fload,bc%coord,IIN_PAR)
+ ! WARNING: the line below is only valid for pure strike-slip faulting
+ bc%V(1,:) = bc%rsf%V_init
+ else
+ allocate(bc%swf)
+ call swf_init(bc%swf,bc%MU,bc%coord,IIN_PAR)
+ if (TPV16) call TPV16_init() !WARNING: ad hoc, initializes T0 and swf
+ endif
+
+ endif
+
+! call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc)
+
+!--------------------------------------------------------
+contains
+
+subroutine TPV16_init
+
+ integer :: ier, ipar
+ integer, parameter :: IIN_NUC =270 ! WARNING: not safe, should look for an available unit
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
+ dyn_fc,swcd,cohes,tim_forcedRup
+ integer, dimension(bc%nglob) :: inp_nx,inp_nz
+ real(kind=CUSTOM_REAL) :: minX, siz_str,siz_dip, hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ integer :: relz_num,sub_relz_num, num_cell_str,num_cell_dip, hypo_cell_str,hypo_cell_dip
+ integer :: i
+
+ open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
+ read(IIN_NUC,*) relz_num,sub_relz_num
+ read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
+ read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ do ipar=1,bc%nglob
+ read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
+ Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
+ enddo
+ close(IIN_NUC)
+
+ minX = minval(bc%coord(1,:))
+
+ do i=1,bc%nglob
+
+ ! WARNING: nearest neighbor interpolation
+ ipar = minloc( (minX+loc_str(:)-bc%coord(1,i))**2 + (-loc_dip(:)-bc%coord(3,i))**2 , 1)
+ !loc_dip is negative of Z-coord
+
+ bc%T0(3,i) = -sigma0(ipar)
+ bc%T0(1,i) = tau0_str(ipar)
+ bc%T0(2,i) = tau0_dip(ipar)
+
+ bc%swf%mus(i) = static_fc(ipar)
+ bc%swf%mud(i) = dyn_fc(ipar)
+ bc%swf%Dc(i) = swcd(ipar)
+ bc%swf%C(i) = cohes(ipar)
+ bc%swf%T(i) = tim_forcedRup(ipar)
+
+ enddo
+
+end subroutine TPV16_init
+
+end subroutine init_one_fault
+
+!---------------------------------------------------------------------
+! REPLACES the value of a fault parameter inside an area with prescribed shape
+subroutine init_2d_distribution(a,coord,iin,n)
+!JPA refactor: background value should be an argument
+
+ real(kind=CUSTOM_REAL), intent(inout) :: a(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: iin,n
+
+ real(kind=CUSTOM_REAL) :: b(size(a))
+ character(len=20) :: shape
+ real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
+ real(kind=CUSTOM_REAL) :: r1(size(a))
+ integer :: i
+
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+
+ if (n==0) return
+
+ do i=1,n
+ shape = ''
+ xc = 0e0_CUSTOM_REAL
+ yc = 0e0_CUSTOM_REAL
+ zc = 0e0_CUSTOM_REAL
+ r = 0e0_CUSTOM_REAL
+ l = 0e0_CUSTOM_REAL
+ lx = 0e0_CUSTOM_REAL
+ ly = 0e0_CUSTOM_REAL
+ lz = 0e0_CUSTOM_REAL
+ valh = 0e0_CUSTOM_REAL
+
+ read(iin,DIST2D)
+ select case(shape)
+ case ('circle')
+ b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
+ case ('circle-exp')
+ r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
+ where(r1<r)
+ b =exp(r1**2/(r1**2 - r**2) )
+ elsewhere
+ b =0._CUSTOM_REAL
+ endwhere
+ case ('ellipse')
+ b = heaviside( 1e0_CUSTOM_REAL - sqrt( (coord(1,:)-xc)**2/lx**2 + (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2 ) )
+ case ('square')
+ b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('cylinder')
+ b = heaviside(r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2)) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle-taper')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
+ case default
+ stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
+ end select
+
+ ! REPLACE the value inside the prescribed area
+ where (b /= 0e0_CUSTOM_REAL) a = b
+ enddo
+
+end subroutine init_2d_distribution
+
+!---------------------------------------------------------------------
+elemental function heaviside(x)
+
+ real(kind=CUSTOM_REAL), intent(in) :: x
+ real(kind=CUSTOM_REAL) :: heaviside
+
+ if (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
+ else
+ heaviside = 0e0_CUSTOM_REAL
+ endif
+
+end function heaviside
+
+!=====================================================================
+! adds boundary term Bt into Force array for each fault.
+!
+! NOTE: On non-split nodes at fault edges, dD=dV=dA=0
+! and the net contribution of B*T is =0
+!
+subroutine bc_dynflt_set3d_all(F,V,D)
+
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V,D
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
+
+ integer :: i
+
+ if (.not. allocated(faults)) return
+ do i=1,size(faults)
+ call BC_DYNFLT_set3d(faults(i),F,V,D,i)
+ enddo
+
+end subroutine bc_dynflt_set3d_all
+
+!---------------------------------------------------------------------
+subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
+
+ use specfem_par, only: it,NSTEP,myrank
+
+ real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
+ type(bc_dynflt_type), intent(inout) :: bc
+ real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
+ integer, intent(in) :: iflt
+
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T,dD,dV,dA
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength,tStick,tnew, &
+ theta_old, theta_new, dc, &
+ ta,Vf_old,Vf_new,TxExt
+ real(kind=CUSTOM_REAL) :: half_dt,TLoad,DTau0,GLoad,time
+ integer :: i
+
+ if (bc%nspec > 0) then !Surendra : for parallel faults
+
+ half_dt = 0.5e0_CUSTOM_REAL*bc%dt
+ Vf_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+
+ ! get predicted values
+ dD = get_jump(bc,D) ! dD_predictor
+ dV = get_jump(bc,V) ! dV_predictor
+ dA = get_weighted_jump(bc,MxA) ! dA_free
+
+ ! rotate to fault frame (tangent,normal)
+ ! component 3 is normal to the fault
+ dD = rotate(bc,dD,1)
+ dV = rotate(bc,dV,1)
+ dA = rotate(bc,dA,1)
+
+ ! T_stick
+ T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
+ T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
+ T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
+
+ !Warning : dirty particular free surface condition z = 0.
+ ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
+ ! end do
+
+ ! add initial stress
+ T = T + bc%T0
+
+ ! Solve for normal stress (negative is compressive)
+ ! Opening implies free stress
+ if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
+
+ ! smooth loading within nucleation patch
+ !WARNING : ad hoc for SCEC benchmark TPV10x
+ if (RATE_AND_STATE) then
+ TxExt = 0._CUSTOM_REAL
+ TLoad = 1.0_CUSTOM_REAL
+ DTau0 = 25e6_CUSTOM_REAL
+ time = it*bc%dt !time will never be zero. it starts from 1
+ if (time <= TLoad) then
+ GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
+ else
+ GLoad = 1.0_CUSTOM_REAL
+ endif
+ TxExt = DTau0 * bc%Fload * GLoad
+ T(1,:) = T(1,:) + TxExt
+ endif
+
+ tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+
+ if (.not. RATE_AND_STATE) then ! Update slip weakening friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ theta_old = bc%swf%theta
+ call swf_update_state(bc%D,dD,bc%V,bc%swf)
+
+ ! Update friction coeficient
+ bc%MU = swf_mu(bc%swf)
+
+ ! combined with time-weakening for nucleation
+ ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ if (TPV16) then
+ where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
+ endif
+
+ ! Update strength
+ strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
+
+ ! Solve for shear stress
+ tnew = min(tStick,strength)
+
+ else ! Update rate and state friction:
+ !JPA the solver below can be refactored into a loop with two passes
+
+ ! first pass
+ theta_old = bc%rsf%theta
+ call rsf_update_state(Vf_old,bc%dt,bc%rsf)
+ do i=1,bc%nglob
+ Vf_new(i)=rtsafe(funcd,0.0,Vf_old(i)+5.0,1e-5,tStick(i),-T(3,i),bc%Z(i),bc%rsf%f0(i), &
+ bc%rsf%V0(i),bc%rsf%a(i),bc%rsf%b(i),bc%rsf%L(i),bc%rsf%theta(i),bc%rsf%StateLaw)
+ enddo
+
+ ! second pass
+ bc%rsf%theta = theta_old
+ call rsf_update_state(0.5e0_CUSTOM_REAL*(Vf_old + Vf_new),bc%dt,bc%rsf)
+ do i=1,bc%nglob
+ Vf_new(i)=rtsafe(funcd,0.0,Vf_old(i)+5.0,1e-5,tStick(i),-T(3,i),bc%Z(i),bc%rsf%f0(i), &
+ bc%rsf%V0(i),bc%rsf%a(i),bc%rsf%b(i),bc%rsf%L(i),bc%rsf%theta(i),bc%rsf%StateLaw)
+ enddo
+
+ tnew = tStick - bc%Z*Vf_new
+
+ endif
+
+ tStick = max(tStick,1e0_CUSTOM_REAL) ! to avoid division by zero
+ T(1,:) = tnew * T(1,:)/tStick
+ T(2,:) = tnew * T(2,:)/tStick
+
+ ! Save total tractions
+ bc%T = T
+
+ ! Subtract initial stress
+ T = T - bc%T0
+
+ if (RATE_AND_STATE) T(1,:) = T(1,:) - TxExt
+ !JPA: this eliminates the effect of TxExt on the equations of motion. Why is it needed?
+
+ ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
+ dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
+ dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
+ dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
+
+ ! Update slip and slip rate, in fault frame
+ bc%D = dD
+ bc%V = dV + half_dt*dA
+
+ ! Rotate tractions back to (x,y,z) frame
+ T = rotate(bc,T,-1)
+
+ ! Add boundary term B*T to M*a
+ call add_BT(bc,MxA,T)
+
+ !-- intermediate storage of outputs --
+ Vf_new = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ if(.not. RATE_AND_STATE) then
+ theta_new = bc%swf%theta
+ dc = bc%swf%dc
+ else
+ theta_new = bc%rsf%theta
+ dc = bc%rsf%L
+ endif
+
+ call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
+ Vf_old, Vf_new, it*bc%dt,bc%dt)
+! call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
+
+ !-- outputs --
+ ! write dataT every NTOUT time step or at the end of simulation
+! if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it,bc%rsf%StateLaw)
+
+ endif
+
+ ! write dataXZ every NSNAP time step
+ if ( mod(it,NSNAP) == 0) then
+ if (.NOT. PARALLEL_FAULT) then
+ if (bc%nspec > 0) call write_dataXZ(bc%dataXZ,it,iflt)
+ else
+ call gather_dataXZ(bc)
+ if (myrank==0) call write_dataXZ(bc%dataXZ_all,it,iflt)
+ endif
+ endif
+
+ if ( it == NSTEP) then
+ if (.NOT. PARALLEL_FAULT) then
+ call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
+ else
+ if (myrank==0) call SCEC_Write_RuptureTime(bc%dataXZ_all,bc%dt,NSTEP,iflt)
+ endif
+ endif
+
+end subroutine BC_DYNFLT_set3d
+
+!===============================================================
+
+subroutine swf_init(f,mu,coord,IIN_PAR)
+
+ type(swf_type), intent(out) :: f
+ real(kind=CUSTOM_REAL), intent(out) :: mu(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: IIN_PAR
+
+ integer :: nglob
+ real(kind=CUSTOM_REAL) :: mus,mud,dc,C,T
+ integer :: nmus,nmud,ndc,nC,nForcedRup
+
+ NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
+
+ nglob = size(mu)
+ allocate( f%mus(nglob) )
+ allocate( f%mud(nglob) )
+ allocate( f%Dc(nglob) )
+ allocate( f%theta(nglob) )
+ allocate( f%C(nglob) )
+ allocate( f%T(nglob) )
+
+ ! WARNING: if V_HEALING is negative we turn off healing
+ f%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+
+ mus = 0.6e0_CUSTOM_REAL
+ mud = 0.1e0_CUSTOM_REAL
+ dc = 1e0_CUSTOM_REAL
+ C = 0._CUSTOM_REAL
+ T = HUGEVAL
+
+ nmus = 0
+ nmud = 0
+ ndc = 0
+ nC = 0
+ nForcedRup = 0
+
+ read(IIN_PAR, nml=SWF)
+
+ f%mus = mus
+ f%mud = mud
+ f%Dc = dc
+ f%C = C
+ f%T = T
+ call init_2d_distribution(f%mus,coord,IIN_PAR,nmus)
+ call init_2d_distribution(f%mud,coord,IIN_PAR,nmud)
+ call init_2d_distribution(f%Dc ,coord,IIN_PAR,ndc)
+ call init_2d_distribution(f%C ,coord,IIN_PAR,nC)
+ call init_2d_distribution(f%T ,coord,IIN_PAR,nForcedRup)
+
+ f%theta = 0e0_CUSTOM_REAL
+
+ mu = swf_mu(f)
+
+end subroutine swf_init
+
+!---------------------------------------------------------------------
+subroutine swf_update_state(dold,dnew,vold,f)
+
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
+ type(swf_type), intent(inout) :: f
+
+ real(kind=CUSTOM_REAL) :: vnorm
+ integer :: k,npoin
+
+ f%theta = f%theta + sqrt( (dold(1,:)-dnew(1,:))**2 + (dold(2,:)-dnew(2,:))**2 )
+
+ if (f%healing) then
+ npoin = size(vold,2)
+ do k=1,npoin
+ vnorm = sqrt(vold(1,k)**2 + vold(2,k)**2)
+ if (vnorm<V_HEALING) f%theta(k) = 0e0_CUSTOM_REAL
+ enddo
+ endif
+end subroutine swf_update_state
+
+!---------------------------------------------------------------------
+function swf_mu(f) result(mu)
+
+ type(swf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+
+ mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
+ mu = max( mu, f%mud)
+
+end function swf_mu
+
+!=====================================================================
+
+subroutine rsf_init(f,T0,nucFload,coord,IIN_PAR)
+
+ type(rsf_type), intent(out) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: T0(:,:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ real(kind=CUSTOM_REAL), pointer :: nucFload(:)
+ integer, intent(in) :: IIN_PAR
+
+ real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init,fw,Vw, C,T
+ integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init,nfw,nVw, nC,nForcedRup
+ real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
+ real(kind=CUSTOM_REAL) :: x,z
+ logical :: c1,c2,c3,c4
+ real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
+ integer :: i !,nglob_bulk
+ real(kind=CUSTOM_REAL) :: Fload
+ integer :: nFload
+! real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: init_vel
+ integer :: nglob
+
+ NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup,Vw,fw,nVw,nfw
+ NAMELIST / ASP / Fload,nFload
+
+ nglob = size(coord,2)
+
+ allocate( f%V0(nglob) )
+ allocate( f%f0(nglob) )
+ allocate( f%a(nglob) )
+ allocate( f%b(nglob) )
+ allocate( f%L(nglob) )
+ allocate( f%V_init(nglob) )
+ allocate( f%theta(nglob) )
+ allocate( f%C(nglob) )
+ allocate( f%T(nglob) )
+ allocate( f%fw(nglob) )
+ allocate( f%Vw(nglob) )
+
+ V0 =1.e-6_CUSTOM_REAL
+ f0 =0.6_CUSTOM_REAL
+ a =0.0080_CUSTOM_REAL !0.0080_CUSTOM_REAL
+ b =0.0040_CUSTOM_REAL !0.0120_CUSTOM_REAL
+ L =0.0135_CUSTOM_REAL
+ V_init =1.e-12_CUSTOM_REAL
+ theta_init =1.084207680000000e+09_CUSTOM_REAL
+ C = 0._CUSTOM_REAL
+ T = HUGEVAL
+ fw = 0.2_CUSTOM_REAL
+ Vw = 0.1_CUSTOM_REAL
+
+ nV0 =0
+ nf0 =0
+ na =0
+ nb =0
+ nL =0
+ nV_init =0
+ ntheta_init =0
+ nC = 0
+ nForcedRup = 0
+ nfw =0
+ nVw =0
+
+ read(IIN_PAR, nml=RSF)
+
+ f%V0 = V0
+ f%f0 = f0
+ f%a = a
+ f%b = b
+ f%L = L
+ f%V_init = V_init
+ f%theta = theta_init
+ f%C = C
+ f%T = T
+ f%fw = fw
+ f%Vw = Vw
+
+ call init_2d_distribution(f%V0,coord,IIN_PAR,nV0)
+ call init_2d_distribution(f%f0,coord,IIN_PAR,nf0)
+ call init_2d_distribution(f%a,coord,IIN_PAR,na)
+ call init_2d_distribution(f%b,coord,IIN_PAR,nb)
+ call init_2d_distribution(f%L,coord,IIN_PAR,nL)
+ call init_2d_distribution(f%V_init,coord,IIN_PAR,nV_init)
+ call init_2d_distribution(f%theta,coord,IIN_PAR,ntheta_init)
+ call init_2d_distribution(f%C,coord,IIN_PAR,nC)
+ call init_2d_distribution(f%T,coord,IIN_PAR,nForcedRup)
+ call init_2d_distribution(f%fw,coord,IIN_PAR,nfw)
+ call init_2d_distribution(f%Vw,coord,IIN_PAR,nVw)
+
+!!$ ! WARNING : Not general enough
+!!$ nglob_bulk = size(vel,2)
+!!$ allocate(init_vel(3,nglob_bulk))
+!!$ init_vel = 0._CUSTOM_REAL
+!!$ init_vel(1,bc%ibulk1) = -f%V_init/2._CUSTOM_REAL
+!!$ init_vel(1,bc%ibulk2) = f%V_init/2._CUSTOM_REAL
+!!$ where(ystore > 0) init_vel(1,:) = -V_init/2._CUSTOM_REAL
+!!$ where(ystore < 0) init_vel(1,:) = V_init/2._CUSTOM_REAL
+!!$ !init_vel = rotate(bc,init_vel,-1) ! directly assigned in global coordinates here as it is a simplified case
+!!$ vel = vel + init_vel
+
+ ! WARNING: below is an ad-hoc setting of a(x,z) for some SCEC benchmark
+ ! This should be instead an option in init_2d_distribution
+ W1=15000._CUSTOM_REAL
+ W2=7500._CUSTOM_REAL
+ w=3000._CUSTOM_REAL
+ hypo_z = -7500._CUSTOM_REAL
+ do i=1,nglob
+ x=coord(1,i)
+ z=coord(3,i)
+ c1=abs(x)<W1+w
+ c2=abs(x)>W1
+ c3=abs(z-hypo_z)<W2+w
+ c4=abs(z-hypo_z)>W2
+ if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
+
+ if (c1 .and. c2) then
+ b11 = w/(abs(x)-W1-w)
+ b12 = w/(abs(x)-W1)
+ B1 = HALF * (ONE + tanh(b11 + b12))
+ elseif(abs(x)<=W1) then
+ B1 = 1._CUSTOM_REAL
+ else
+ B1 = 0._CUSTOM_REAL
+ endif
+
+ if (c3 .and. c4) then
+ b21 = w/(abs(z-hypo_z)-W2-w)
+ b22 = w/(abs(z-hypo_z)-W2)
+ B2 = HALF * (ONE + tanh(b21 + b22))
+ elseif(abs(z-hypo_z)<=W2) then
+ B2 = 1._CUSTOM_REAL
+ else
+ B2 = 0._CUSTOM_REAL
+ endif
+
+ f%a(i) = 0.008 + 0.008 * (ONE - B1*B2)
+ f%Vw(i) = 0.1 + 0.9 * (ONE - B1*B2)
+
+ elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
+ f%a(i) = 0.008
+ f%Vw(i) = 0.1_CUSTOM_REAL
+ else
+ f%a(i) = 0.016
+ f%Vw(i) = 1.0_CUSTOM_REAL
+ endif
+
+ enddo
+
+ ! WARNING: The line below scratches an earlier initialization of theta through theta_init
+ ! We should implement it as an option for the user
+ if(f%stateLaw == 1) then
+ f%theta = f%L/f%V0 &
+ * exp( ( f%a * log(TWO*sinh(-sqrt(T0(1,:)**2+T0(2,:)**2)/T0(3,:)/f%a)) &
+ - f%f0 - f%a*log(f%V_init/f%V0) ) &
+ / f%b )
+ else
+ f%theta = f%a * log(TWO*V0/V_init * sinh(-sqrt(T0(1,:)**2+T0(2,:)**2)/T0(3,:)/f%a))
+ endif
+
+ ! WARNING : ad hoc for SCEC benchmark TPV10x
+ allocate( nucFload(nglob) )
+ Fload = 0.e0_CUSTOM_REAL
+ nFload = 0
+ read(IIN_PAR, nml=ASP)
+ nucFload = Fload
+ call init_2d_distribution(nucFload,coord,IIN_PAR,nFload)
+
+end subroutine rsf_init
+
+!---------------------------------------------------------------------
+! Rate and state friction coefficient
+function rsf_mu(f,V) result(mu)
+
+ type(rsf_type), intent(in) :: f
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ real(kind=CUSTOM_REAL) :: mu(size(V))
+ mu = f%a * asinh( V/2.0/f%V0 * exp((f%f0 + f%b*log(f%theta*f%V0/f%L))/f%a ) ) ! Regularized
+
+end function rsf_mu
+
+!---------------------------------------------------------------------
+subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
+
+ real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+ double precision :: arg,fn,df,x
+ integer :: statelaw
+
+ if(statelaw == 1) then
+ arg = exp((f0+dble(b)*log(V0*theta/L))/a)/TWO/V0
+ else
+ arg = exp(theta/a)/TWO/V0
+ endif
+ fn = tStick - Z*x - a*Seff*asinh(x*arg)
+ df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
+end subroutine funcd
+
+!---------------------------------------------------------------------
+function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
+
+ integer, parameter :: MAXIT=200
+ real(kind=CUSTOM_REAL) :: x1,x2,xacc
+ EXTERNAL funcd
+ integer :: j
+ !real(kind=CUSTOM_REAL) :: df,dx,dxold,f,fh,fl,temp,xh,xl
+ double precision :: df,dx,dxold,f,fh,fl,temp,xh,xl,rtsafe
+ real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+ integer :: statelaw
+
+ call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ if( (fl>0 .and. fh>0) .or. (fl<0 .and. fh<0) ) stop 'root must be bracketed in rtsafe'
+ if(fl==0.) then
+ rtsafe=x2
+ return
+ elseif(fh==0.) then
+ rtsafe=x2
+ return
+ elseif(fl<0) then
+ xl=x1
+ xh=x2
+ else
+ xh=x1
+ xl=x2
+ endif
+
+ rtsafe=0.5d0*(x1+x2)
+ dxold=abs(x2-x1)
+ dx=dxold
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ do j=1,MAXIT
+ if( ((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0 .or. abs(2.*f)>abs(dxold*df) ) then
+ dxold=dx
+ dx=0.5d0*(xh-xl)
+ rtsafe=xl+dx
+ if(xl==rtsafe) return
+ else
+ dxold=dx
+ dx=f/df
+ temp=rtsafe
+ rtsafe=rtsafe-dx
+ if(temp==rtsafe) return
+ endif
+ if(abs(dx)<xacc) return
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ if(f<0.) then
+ xl=rtsafe
+ else
+ xh=rtsafe
+ endif
+ enddo
+ stop 'rtsafe exceeding maximum iterations'
+ return
+
+end function rtsafe
+
+!---------------------------------------------------------------------
+subroutine rsf_update_state(V,dt,f)
+
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ type(rsf_type), intent(inout) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+
+ real(kind=CUSTOM_REAL) :: vDtL(size(V))
+ real(kind=CUSTOM_REAL) :: f_ss(size(V)),xi_ss(size(V)),fLV(size(V))
+
+ vDtL = V * dt / f%L
+
+ ! ageing law
+ if (f%StateLaw == 1) then
+ where(vDtL > 1.e-5_CUSTOM_REAL)
+ f%theta = f%theta * exp(-vDtL) + f%L/V * (ONE - exp(-vDtL))
+ elsewhere
+ f%theta = f%theta * exp(-vDtL) + dt*( ONE - HALF*vDtL )
+ endwhere
+
+ ! slip law : by default use strong rate-weakening
+ else
+! f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+ where(V /= 0._CUSTOM_REAL)
+ fLV = f%f0 - (f%b - f%a)*log(V/f%V0)
+ f_ss = f%fw + (fLV - f%fw)/(ONE + (V/f%Vw)**8)**0.125
+ xi_ss = f%a * log( TWO*f%V0/V * sinh(f_ss/f%a) )
+ f%theta = xi_ss + (f%theta - xi_ss) * exp(-vDtL)
+ elsewhere
+ f%theta = f%theta
+ endwhere
+ endif
+
+end subroutine rsf_update_state
+
+
+!===============================================================
+! OUTPUTS
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+
+ use specfem_par, only : NPROC,myrank
+ ! NT = total number of time steps
+
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
+
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
+
+ integer :: ipoin, ipoin_local, iproc
+
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
+
+ IIN = 251
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ read(IIN,*) np
+ DataT%npoin =0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt==iflt) DataT%npoin = DataT%npoin +1
+ enddo
+ close(IIN)
+
+ if (DataT%npoin == 0) return
+
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
+ allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
+
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
+ read(IIN,*) np
+ k = 0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt/=iflt) cycle
+ k = k+1
+ DataT%name(k) = tmpname
+ !search nearest node
+ distkeep = huge(distkeep)
+
+ do iglob=1,nglob
+ dist = sqrt((coord(1,iglob)-xtarget)**2 &
+ + (coord(2,iglob)-ytarget)**2 &
+ + (coord(3,iglob)-ztarget)**2)
+ if (dist < distkeep) then
+ distkeep = dist
+ DataT%iglob(k) = iglob
+ endif
+ enddo
+ DataT%dist(k) = distkeep !Surendra : for parallel fault
+ enddo
+
+ !Surendra : for parallel fault
+ if (PARALLEL_FAULT) then
+ allocate(DataT%islice(DataT%npoin))
+ allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
+ allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
+ call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
+ call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
+ if(myrank==0) then
+ do ipoin = 1,DataT%npoin
+ iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
+ DataT%islice(ipoin) = iproc
+ DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
+ enddo
+ endif
+
+ call bcast_all_i(DataT%islice,DataT%npoin)
+ call bcast_all_i(DataT%iglob,DataT%npoin)
+
+ DataT%npoin_local = 0
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
+ enddo
+ allocate(DataT%glob_indx(DataT%npoin_local))
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) then
+ ipoin_local = ipoin_local + 1
+ DataT%glob_indx(ipoin_local) = ipoin
+ endif
+ enddo
+ else
+ DataT%npoin_local = DataT%npoin
+ endif !Parallel_fault
+
+
+ ! 3. allocate arrays and set to zero
+ allocate(DataT%d1(NT,DataT%npoin))
+ allocate(DataT%v1(NT,DataT%npoin))
+ allocate(DataT%t1(NT,DataT%npoin))
+ allocate(DataT%d2(NT,DataT%npoin))
+ allocate(DataT%v2(NT,DataT%npoin))
+ allocate(DataT%t2(NT,DataT%npoin))
+ allocate(DataT%t3(NT,DataT%npoin))
+ allocate(DataT%theta(NT,DataT%npoin))
+ DataT%d1 = 0e0_CUSTOM_REAL
+ DataT%v1 = 0e0_CUSTOM_REAL
+ DataT%t1 = 0e0_CUSTOM_REAL
+ DataT%d2 = 0e0_CUSTOM_REAL
+ DataT%v2 = 0e0_CUSTOM_REAL
+ DataT%t2 = 0e0_CUSTOM_REAL
+ DataT%t3 = 0e0_CUSTOM_REAL
+ DataT%theta = 0e0_CUSTOM_REAL
+
+ close(IIN)
+
+end subroutine init_dataT
+
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,theta,itime)
+
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
+ integer, intent(in) :: itime
+
+ integer :: i,k,ipoin
+
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+ k = dataT%iglob(i)
+ dataT%d1(itime,i) = d(1,k)
+ dataT%d2(itime,i) = d(2,k)
+ dataT%v1(itime,i) = v(1,k)
+ dataT%v2(itime,i) = v(2,k)
+ dataT%t1(itime,i) = t(1,k)
+ dataT%t2(itime,i) = t(2,k)
+ dataT%t3(itime,i) = t(3,k)
+ dataT%theta(itime,i) = theta(k)
+ enddo
+
+end subroutine store_dataT
+
+!-----------------------------------------------------------------
+subroutine write_dataT_all(nt,statelaw)
+
+ integer, intent(in) :: nt
+ integer, intent(in) :: statelaw
+
+ integer :: i
+
+ if (.not.allocated(faults)) return
+ do i = 1,size(faults)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
+ enddo
+
+end subroutine write_dataT_all
+
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
+
+ type(dataT_type), intent(in) :: dataT
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT
+ integer, intent(in) :: statelaw
+
+ integer :: i,k,IOUT,ipoin
+ character :: NTchar*5
+ integer :: today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
+
+1 format(I5)
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+
+ open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ write(IOUT,*) "# problem=TPV104"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
+ write(IOUT,*) "# code_version=1.1"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# location=",trim(dataT%name(i))
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ write(IOUT,*) "# Column #8 = normal stress (MPa)"
+ if (RATE_AND_STATE) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
+ write(IOUT,*) "#"
+ write(IOUT,*) "# The line below lists the names of the data fields:"
+ if (.not. RATE_AND_STATE) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL
+ enddo
+ else
+ if(statelaw == 1) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
+ enddo
+ else
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress psi"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, dataT%theta(k,i)
+ enddo
+ endif
+ endif
+ close(IOUT)
+
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
+
+ enddo
+
+end subroutine SCEC_write_dataT
+
+!-------------------------------------------------------------------------------------------------
+
+subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
+
+ type(dataXZ_type), intent(in) :: dataXZ
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT,iflt
+
+ integer :: i,IOUT
+ character(len=70) :: filename
+ integer*4 today(3), now(3)
+
+ call idate(today) ! today(1)=day, (2)=month, (3)=year
+ call itime(now) ! now(1)=hour, (2)=minute, (3)=second
+
+
+ write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ open(IOUT,file=trim(filename),status='replace')
+ write(IOUT,*) "# problem=TPV104"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ write(IOUT,1000) today(2), today(1), today(3), now
+ write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
+ write(IOUT,*) "# code_version=1.1"
+ write(IOUT,*) "# element_size=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# Column #1 = horizontal coordinate, distance along strike (m)"
+ write(IOUT,*) "# Column #2 = vertical coordinate, distance down-dip (m)"
+ write(IOUT,*) "# Column #3 = rupture time (s)"
+ write(IOUT,*) "# "
+ write(IOUT,*) "j k t"
+ do i = 1,size(dataXZ%tRUP)
+ write(IOUT,'(3(E15.7))') dataXZ%xcoord(i), -dataXZ%zcoord(i), dataXZ%tRUP(i)
+ end do
+
+ close(IOUT)
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
+end subroutine SCEC_Write_RuptureTime
+
+!-------------------------------------------------------------------------------------------------
+subroutine init_dataXZ(DataXZ,bc)
+
+ use specfem_par, only : NPROC,myrank
+
+ type(dataXZ_type), intent(inout) :: DataXZ
+ type(bc_dynflt_type) :: bc
+
+ integer :: npoin_all,iproc
+
+ DataXZ%npoin = bc%nglob
+
+ if(bc%nglob > 0) then
+
+ allocate(DataXZ%stg(bc%nglob))
+ if(.not. RATE_AND_STATE) then
+ DataXZ%sta => bc%swf%theta
+ else
+ DataXZ%sta => bc%rsf%theta
+ endif
+ DataXZ%d1 => bc%d(1,:)
+ DataXZ%d2 => bc%d(2,:)
+ DataXZ%v1 => bc%v(1,:)
+ DataXZ%v2 => bc%v(2,:)
+ DataXZ%t1 => bc%t(1,:)
+ DataXZ%t2 => bc%t(2,:)
+ DataXZ%t3 => bc%t(3,:)
+ DataXZ%xcoord => bc%coord(1,:)
+ DataXZ%ycoord => bc%coord(2,:)
+ DataXZ%zcoord => bc%coord(3,:)
+ allocate(DataXZ%tRUP(bc%nglob))
+ allocate(DataXZ%tPZ(bc%nglob))
+
+ !Percy , setting up initial rupture time null for all faults.
+ DataXZ%tRUP = 0e0_CUSTOM_REAL
+ DataXZ%tPZ = 0e0_CUSTOM_REAL
+
+ endif
+
+ !Surendra : for parallel fault
+ if (PARALLEL_FAULT) then
+ npoin_all = 0
+ call sum_all_i(bc%nglob,npoin_all)
+ if (myrank==0 .and. npoin_all>0) then
+ bc%DataXZ_all%npoin = npoin_all
+ allocate(bc%DataXZ_all%xcoord(npoin_all))
+ allocate(bc%DataXZ_all%ycoord(npoin_all))
+ allocate(bc%DataXZ_all%zcoord(npoin_all))
+ allocate(bc%DataXZ_all%t1(npoin_all))
+ allocate(bc%DataXZ_all%t2(npoin_all))
+ allocate(bc%DataXZ_all%t3(npoin_all))
+ allocate(bc%DataXZ_all%d1(npoin_all))
+ allocate(bc%DataXZ_all%d2(npoin_all))
+ allocate(bc%DataXZ_all%v1(npoin_all))
+ allocate(bc%DataXZ_all%v2(npoin_all))
+ allocate(bc%DataXZ_all%tRUP(npoin_all))
+ allocate(bc%DataXZ_all%tPZ(npoin_all))
+ allocate(bc%DataXZ_all%stg(npoin_all))
+ allocate(bc%DataXZ_all%sta(npoin_all))
+ endif
+
+ allocate(bc%npoin_perproc(NPROC))
+ bc%npoin_perproc=0
+ call gather_all_i(DataXZ%npoin,1,bc%npoin_perproc,1,NPROC)
+
+ allocate(bc%poin_offset(NPROC))
+ bc%poin_offset(1)=0
+ do iproc=2,NPROC
+ bc%poin_offset(iproc) = sum(bc%npoin_perproc(1:iproc-1))
+ enddo
+
+ call gatherv_all_cr(DataXZ%xcoord,DataXZ%npoin,bc%DataXZ_all%xcoord,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(DataXZ%ycoord,DataXZ%npoin,bc%DataXZ_all%ycoord,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(DataXZ%zcoord,DataXZ%npoin,bc%DataXZ_all%zcoord,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ endif
+
+end subroutine init_dataXZ
+!---------------------------------------------------------------
+
+subroutine gather_dataXZ(bc)
+
+ use specfem_par, only : NPROC
+
+ type(bc_dynflt_type), intent(inout) :: bc
+
+ call gatherv_all_cr(bc%DataXZ%t1,bc%DataXZ%npoin,bc%DataXZ_all%t1,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%t2,bc%DataXZ%npoin,bc%DataXZ_all%t2,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%t3,bc%DataXZ%npoin,bc%DataXZ_all%t3,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%d1,bc%DataXZ%npoin,bc%DataXZ_all%d1,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%d2,bc%DataXZ%npoin,bc%DataXZ_all%d2,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%v1,bc%DataXZ%npoin,bc%DataXZ_all%v1,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%v2,bc%DataXZ%npoin,bc%DataXZ_all%v2,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%tRUP,bc%DataXZ%npoin,bc%DataXZ_all%tRUP,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%tPZ,bc%DataXZ%npoin,bc%DataXZ_all%tPZ,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%stg,bc%DataXZ%npoin,bc%DataXZ_all%stg,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+ call gatherv_all_cr(bc%DataXZ%sta,bc%DataXZ%npoin,bc%DataXZ_all%sta,bc%npoin_perproc,bc%poin_offset,bc%DataXZ_all%npoin,NPROC)
+
+end subroutine gather_dataXZ
+
+!---------------------------------------------------------------
+subroutine store_dataXZ(dataXZ,stg,dold,dnew,dc,vold,vnew,time,dt)
+
+ type(dataXZ_type), intent(inout) :: dataXZ
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: stg,dold,dnew,dc,vold,vnew
+ real(kind=CUSTOM_REAL), intent(in) :: time,dt
+
+ integer :: i
+
+ ! "stg" : strength .
+
+ dataXZ%stg = stg
+
+ do i = 1,size(stg)
+ ! process zone time = first time when slip = dc (break down process).
+ ! with linear time interpolation
+ if (dataXZ%tPZ(i)==0e0_CUSTOM_REAL) then
+ if (dold(i)<=dc(i) .and. dnew(i) >= dc(i)) then
+ dataXZ%tPZ(i) = time-dt*(dnew(i)-dc(i))/(dnew(i)-dold(i))
+ endif
+ endif
+ ! rupture time = first time when slip velocity = vc
+ ! with linear time interpolation
+ ! vc should be pre-defined as input data .
+
+ if (dataXZ%tRUP(i)==0e0_CUSTOM_REAL) then
+ if (vold(i)<=V_RUPT .and. vnew(i)>=V_RUPT) dataXZ%tRUP(i)= time-dt*(vnew(i)-V_RUPT)/(vnew(i)-vold(i))
+ endif
+ enddo
+
+ ! To do : add stress criteria (firs time strength is reached).
+
+ ! note: the other arrays in dataXZ are pointers to arrays in bc
+ ! they do not need to be updated here
+
+end subroutine store_dataXZ
+
+!---------------------------------------------------------------
+subroutine write_dataXZ(dataXZ,itime,iflt)
+
+ type(dataXZ_type), intent(in) :: dataXZ
+ integer, intent(in) :: itime,iflt
+
+ character(len=70) :: filename
+
+ write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
+ ! open(unit=IOUT, file= trim(filename), status='replace', form='formatted',action='write')
+ ! NOTE : It had to be adopted formatted output to avoid conflicts readings with different
+ ! compilers.
+
+ ! write(IOUT,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
+
+ open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
+
+ write(IOUT) dataXZ%xcoord
+ write(IOUT) dataXZ%ycoord
+ write(IOUT) dataXZ%zcoord
+ write(IOUT) dataXZ%d1
+ write(IOUT) dataXZ%d2
+ write(IOUT) dataXZ%v1
+ write(IOUT) dataXZ%v2
+ write(IOUT) dataXZ%t1
+ write(IOUT) dataXZ%t2
+ write(IOUT) dataXZ%t3
+ write(IOUT) dataXZ%sta
+ write(IOUT) dataXZ%stg
+ write(IOUT) dataXZ%tRUP
+ write(IOUT) dataXZ%tPZ
+ close(IOUT)
+
+end subroutine write_dataXZ
+
+
+end module fault_solver_dynamic
More information about the CIG-COMMITS
mailing list