[cig-commits] r21622 - in seismo/3D/SPECFEM3D/trunk: examples/CPML/purely_elastic/HOMO8_FREE_SURFACE examples/CPML/purely_elastic/HOMO8_NOFREESURFACE examples/homogeneous_halfspace_HEX8 examples/meshfem3D_examples/simple_model src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Sun Mar 24 11:58:32 PDT 2013
Author: joseph.charles
Date: 2013-03-24 11:58:31 -0700 (Sun, 24 Mar 2013)
New Revision: 21622
Added:
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/process.sh
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh
Removed:
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt
Modified:
seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace_HEX8/process.sh
seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/process.sh
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
Log:
fixes a bug for runs without C-PML conditions; creates bash scripts for both CPML/purely_elastic/ examples; modifies loops using the CPML_mask_ibool array; comments out boundary conditions for acoustic pml elements
Deleted: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt 2013-03-24 18:58:31 UTC (rev 21622)
@@ -1,11 +0,0 @@
-
-To run this example on 8 processor cores for instance, type this (replace "8" with another value in the three lines below if you want to use more processor cores):
-
-cd bin/
-
-./xdecompose_mesh 8 ../DATA/MESH ../OUTPUT_FILES/DATABASES_MPI
-
-mpirun -np 8 ./xgenerate_databases
-
-mpirun -np 8 ./xspecfem3D
-
Added: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/process.sh (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_FREE_SURFACE/process.sh 2013-03-24 18:58:31 UTC (rev 21622)
@@ -0,0 +1,91 @@
+#!/bin/bash
+#
+# script runs decomposition,database generation and solver
+# using this example setup
+#
+# prior to running this script, you must create the mesh files
+# in directory MESH/
+# (see section 3.1 "Meshing with CUBIT" in user guide)
+#
+
+##################################################
+
+# number of processes
+NPROC=8
+
+##################################################
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take about 5 minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo " setting up example..."
+echo
+
+mkdir -p bin
+mkdir -p OUTPUT_FILES/DATABASES_MPI
+
+rm -f OUTPUT_FILES/*
+rm -rf OUTPUT_FILES/DATABASES_MPI/*
+
+# compiles executables in root directory
+cd ../../../..
+
+rm -fr DATA/*
+cd $currentdir
+cp -fr DATA/* ../../../../DATA/.
+
+cd ../../../..
+
+make clean
+./configure
+make all > $currentdir/tmp.log
+cd $currentdir
+
+# links executables
+cd bin/
+rm -f *
+ln -s ../../../../../bin/xdecompose_mesh .
+ln -s ../../../../../bin/xgenerate_databases .
+ln -s ../../../../../bin/xspecfem3D .
+cd ../
+
+if [ ! -e bin/xspecfem3D ]; then echo "compilation failed, please check..."; exit 1; fi
+
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/CMTSOLUTION OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+
+# decomposes mesh
+echo
+echo " decomposing mesh..."
+echo
+cd bin/
+./xdecompose_mesh $NPROC ../DATA/MESH/ ../OUTPUT_FILES/DATABASES_MPI/
+
+# runs database generation
+echo
+echo " running database generation..."
+echo
+mpirun -n $NPROC ./xgenerate_databases
+
+# runs simulation
+echo
+echo " running solver..."
+echo
+mpirun -n $NPROC ./xspecfem3D
+cd ../
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
+
+
Deleted: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt 2013-03-24 18:58:31 UTC (rev 21622)
@@ -1,11 +0,0 @@
-
-To run this example on 8 processor cores for instance, type this (replace "8" with another value in the three lines below if you want to use more processor cores):
-
-cd bin/
-
-./xdecompose_mesh 8 ../DATA/MESH ../OUTPUT_FILES/DATABASES_MPI
-
-mpirun -np 8 ./xgenerate_databases
-
-mpirun -np 8 ./xspecfem3D
-
Added: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh 2013-03-24 18:58:31 UTC (rev 21622)
@@ -0,0 +1,91 @@
+#!/bin/bash
+#
+# script runs decomposition,database generation and solver
+# using this example setup
+#
+# prior to running this script, you must create the mesh files
+# in directory MESH/
+# (see section 3.1 "Meshing with CUBIT" in user guide)
+#
+
+##################################################
+
+# number of processes
+NPROC=8
+
+##################################################
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take about 5 minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo " setting up example..."
+echo
+
+mkdir -p bin
+mkdir -p OUTPUT_FILES/DATABASES_MPI
+
+rm -f OUTPUT_FILES/*
+rm -rf OUTPUT_FILES/DATABASES_MPI/*
+
+# compiles executables in root directory
+cd ../../../..
+
+rm -fr DATA/*
+cd $currentdir
+cp -fr DATA/* ../../../../DATA/.
+
+cd ../../../..
+
+make clean
+./configure
+make all > $currentdir/tmp.log
+cd $currentdir
+
+# links executables
+cd bin/
+rm -f bin/*
+ln -s ../../../../../bin/xdecompose_mesh .
+ln -s ../../../../../bin/xgenerate_databases .
+ln -s ../../../../../bin/xspecfem3D .
+cd ../
+
+if [ ! -e bin/xspecfem3D ]; then echo "compilation failed, please check..."; exit 1; fi
+
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/CMTSOLUTION OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+
+# decomposes mesh
+echo
+echo " decomposing mesh..."
+echo
+cd bin/
+./xdecompose_mesh $NPROC ../DATA/MESH/ ../OUTPUT_FILES/DATABASES_MPI/
+
+# runs database generation
+echo
+echo " running database generation..."
+echo
+mpirun -n $NPROC ./xgenerate_databases
+
+# runs simulation
+echo
+echo " running solver..."
+echo
+mpirun -n $NPROC ./xspecfem3D
+cd ../
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
+
+
Modified: seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace_HEX8/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace_HEX8/process.sh 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/examples/homogeneous_halfspace_HEX8/process.sh 2013-03-24 18:58:31 UTC (rev 21622)
@@ -5,6 +5,7 @@
#
# prior to running this script, you must create the mesh files
# in directory MESH/
+# (see section 3.1 "Meshing with CUBIT" in user guide)
#
###################################################
@@ -27,49 +28,54 @@
echo
mkdir -p bin
-mkdir -p OUTPUT_FILES
mkdir -p OUTPUT_FILES/DATABASES_MPI
-rm -rf OUTPUT_FILES/*
+rm -f OUTPUT_FILES/*
rm -rf OUTPUT_FILES/DATABASES_MPI/*
# compiles executables in root directory
cd ../../
-make > tmp.log
+
+rm -fr DATA/*
cd $currentdir
+cp -fr DATA/* ../../DATA/.
+cd ../../
+
+make clean
+./configure
+make all > $currentdir/tmp.log
+cd $currentdir
+
# links executables
cd bin/
-rm -f ./x*
-cp ../../../bin/xdecompose_mesh ./
-cp ../../../bin/xgenerate_databases ./
-cp ../../../bin/xspecfem3D ./
-cd ../
+rm -f *
+cp ../../../bin/xdecompose_mesh .
+cp ../../../bin/xgenerate_databases .
+cp ../../../bin/xspecfem3D .
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/CMTSOLUTION OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+
# decomposes mesh
echo
echo " decomposing mesh..."
echo
-./bin/xdecompose_mesh $NPROC DATA/MESH-default/ OUTPUT_FILES/DATABASES_MPI/
+cd bin/
+./xdecompose_mesh $NPROC ../DATA/MESH-default/ ../OUTPUT_FILES/DATABASES_MPI/
-# stores setup
-cp DATA/Par_file OUTPUT_FILES/
-cp DATA/CMTSOLUTION OUTPUT_FILES/
-cp DATA/STATIONS OUTPUT_FILES/
-
# runs database generation
echo
echo " running database generation..."
echo
-cd bin/
mpirun -np $NPROC ./xgenerate_databases
-cd ../
# runs simulation
echo
echo " running solver..."
echo
-cd bin/
mpirun -np $NPROC ./xspecfem3D
cd ../
Modified: seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/process.sh 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/examples/meshfem3D_examples/simple_model/process.sh 2013-03-24 18:58:31 UTC (rev 21622)
@@ -61,13 +61,11 @@
echo
cd bin/
mpirun -np $NPROC ./xgenerate_databases
-cd ../
# runs simulation
echo
echo " running solver..."
echo
-cd bin/
mpirun -np $NPROC ./xspecfem3D
cd ../
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_noDev.f90 2013-03-24 18:58:31 UTC (rev 21622)
@@ -147,28 +147,32 @@
temp3l = temp3l + chi_elem(i,j,l)*hprime_zz(k,l)
enddo
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- temp1l_new = temp1l
- temp2l_new = temp2l
- temp3l_new = temp3l
-
- do l=1,NGLLX
- hp1 = hprime_xx(l,i)
- iglob = ibool(l,j,k,ispec)
- temp1l_new = temp1l_new + deltat*potential_dot_acoustic(iglob)*hp1
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ temp1l_new = temp1l
+ temp2l_new = temp2l
+ temp3l_new = temp3l
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(l,i)
+ iglob = ibool(l,j,k,ispec)
+ temp1l_new = temp1l_new + deltat*potential_dot_acoustic(iglob)*hp1
!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
+
!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- hp2 = hprime_yy(l,j)
- iglob = ibool(i,l,k,ispec)
- temp2l_new = temp2l_new + deltat*potential_dot_acoustic(iglob)*hp2
+ hp2 = hprime_yy(l,j)
+ iglob = ibool(i,l,k,ispec)
+ temp2l_new = temp2l_new + deltat*potential_dot_acoustic(iglob)*hp2
!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
+
!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- hp3 = hprime_zz(l,k)
- iglob = ibool(i,j,l,ispec)
- temp3l_new = temp3l_new + deltat*potential_dot_acoustic(iglob)*hp3
- enddo
+ hp3 = hprime_zz(l,k)
+ iglob = ibool(i,j,l,ispec)
+ temp3l_new = temp3l_new + deltat*potential_dot_acoustic(iglob)*hp3
+ enddo
+ endif
endif
! get derivatives of potential with respect to x, y and z
@@ -189,20 +193,24 @@
dpotentialdzl = xizl*temp1l + etazl*temp2l + gammazl*temp3l
! stores derivatives of ux, uy and uz with respect to x, y and z
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- ispec_CPML = spec_to_CPML(ispec)
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ ispec_CPML = spec_to_CPML(ispec)
- PML_dpotential_dxl(i,j,k,ispec_CPML) = dpotentialdxl
- PML_dpotential_dyl(i,j,k,ispec_CPML) = dpotentialdyl
- PML_dpotential_dzl(i,j,k,ispec_CPML) = dpotentialdzl
+ PML_dpotential_dxl(i,j,k,ispec_CPML) = dpotentialdxl
+ PML_dpotential_dyl(i,j,k,ispec_CPML) = dpotentialdyl
+ PML_dpotential_dzl(i,j,k,ispec_CPML) = dpotentialdzl
- dpotentialdxl_new = xixl*temp1l_new + etaxl*temp2l_new + gammaxl*temp3l_new
- dpotentialdyl_new = xiyl*temp1l_new + etayl*temp2l_new + gammayl*temp3l_new
- dpotentialdzl_new = xizl*temp1l_new + etazl*temp2l_new + gammazl*temp3l_new
+ dpotentialdxl_new = xixl*temp1l_new + etaxl*temp2l_new + gammaxl*temp3l_new
+ dpotentialdyl_new = xiyl*temp1l_new + etayl*temp2l_new + gammayl*temp3l_new
+ dpotentialdzl_new = xizl*temp1l_new + etazl*temp2l_new + gammazl*temp3l_new
- PML_dpotential_dxl_new(i,j,k,ispec_CPML) = dpotentialdxl_new
- PML_dpotential_dyl_new(i,j,k,ispec_CPML) = dpotentialdyl_new
- PML_dpotential_dzl_new(i,j,k,ispec_CPML) = dpotentialdzl_new
+ PML_dpotential_dxl_new(i,j,k,ispec_CPML) = dpotentialdxl_new
+ PML_dpotential_dyl_new(i,j,k,ispec_CPML) = dpotentialdyl_new
+ PML_dpotential_dzl_new(i,j,k,ispec_CPML) = dpotentialdzl_new
+ endif
endif
! density (reciproc)
@@ -220,15 +228,19 @@
enddo
enddo
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- ! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
- call pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
- sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,xixl,xiyl,xizl, &
- etaxl,etayl,etazl,gammaxl,gammayl,gammazl)
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ ! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
+ call pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
+ tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
+ sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,xixl,xiyl,xizl, &
+ etaxl,etayl,etazl,gammaxl,gammayl,gammazl)
- ! calculates contribution from each C-PML element to update acceleration
- call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+ ! calculates contribution from each C-PML element to update acceleration
+ call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+ endif
endif
! second double-loop over GLL to compute all the terms
@@ -255,9 +267,13 @@
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - ( temp1l + temp2l + temp3l )
! updates potential_dot_dot_acoustic with contribution from each C-PML element
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
- potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML)
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - &
+ potential_dot_dot_acoustic_CPML(i,j,k,ispec_CPML)
+ endif
endif
enddo
@@ -285,111 +301,111 @@
!! DK DK thus test something like: if (is_CPML(ispec) .and. acoustic(ispec)) then
!! DK DK or something like that
- ! C-PML boundary
- if( PML_CONDITIONS ) then
- ! xmin
- do ispec2D=1,nspec2D_xmin
- ispec = ibelm_xmin(ispec2D)
+!!$ ! C-PML boundary
+!!$ if( PML_CONDITIONS ) then
+!!$ ! xmin
+!!$ do ispec2D=1,nspec2D_xmin
+!!$ ispec = ibelm_xmin(ispec2D)
+!!$
+!!$ i = 1
+!!$
+!!$ do k=1,NGLLZ
+!!$ do j=1,NGLLY
+!!$ iglob = ibool(i,j,k,ispec)
+!!$
+!!$ potential_dot_dot_acoustic(iglob) = 0.d0
+!!$ potential_dot_acoustic(iglob) = 0.d0
+!!$ potential_acoustic(iglob) = 0.d0
+!!$ enddo
+!!$ enddo
+!!$ enddo
+!!$
+!!$ ! xmax
+!!$ do ispec2D=1,nspec2D_xmax
+!!$ ispec = ibelm_xmax(ispec2D)
+!!$
+!!$ i = NGLLX
+!!$
+!!$ do k=1,NGLLZ
+!!$ do j=1,NGLLY
+!!$ iglob = ibool(i,j,k,ispec)
+!!$
+!!$ potential_dot_dot_acoustic(iglob) = 0.d0
+!!$ potential_dot_acoustic(iglob) = 0.d0
+!!$ potential_acoustic(iglob) = 0.d0
+!!$ enddo
+!!$ enddo
+!!$ enddo
+!!$
+!!$ ! ymin
+!!$ do ispec2D=1,nspec2D_ymin
+!!$ ispec = ibelm_ymin(ispec2D)
+!!$
+!!$ j = 1
+!!$
+!!$ do k=1,NGLLZ
+!!$ do i=1,NGLLX
+!!$ iglob = ibool(i,j,k,ispec)
+!!$
+!!$ potential_dot_dot_acoustic(iglob) = 0.d0
+!!$ potential_dot_acoustic(iglob) = 0.d0
+!!$ potential_acoustic(iglob) = 0.d0
+!!$ enddo
+!!$ enddo
+!!$ enddo
+!!$
+!!$ ! ymax
+!!$ do ispec2D=1,nspec2D_ymax
+!!$ ispec = ibelm_ymax(ispec2D)
+!!$
+!!$ j = NGLLY
+!!$
+!!$ do k=1,NGLLZ
+!!$ do i=1,NGLLX
+!!$ iglob = ibool(i,j,k,ispec)
+!!$
+!!$ potential_dot_dot_acoustic(iglob) = 0.d0
+!!$ potential_dot_acoustic(iglob) = 0.d0
+!!$ potential_acoustic(iglob) = 0.d0
+!!$ enddo
+!!$ enddo
+!!$ enddo
+!!$
+!!$ ! bottom (zmin)
+!!$ do ispec2D=1,NSPEC2D_BOTTOM
+!!$ ispec = ibelm_bottom(ispec2D)
+!!$
+!!$ k = 1
+!!$
+!!$ do j=1,NGLLY
+!!$ do i=1,NGLLX
+!!$ iglob = ibool(i,j,k,ispec)
+!!$
+!!$ potential_dot_dot_acoustic(iglob) = 0.d0
+!!$ potential_dot_acoustic(iglob) = 0.d0
+!!$ potential_acoustic(iglob) = 0.d0
+!!$ enddo
+!!$ enddo
+!!$ enddo
+!!$
+!!$ ! top (zmax)
+!!$ do ispec2D=1,NSPEC2D_BOTTOM
+!!$ ispec = ibelm_top(ispec2D)
+!!$
+!!$ k = NGLLZ
+!!$
+!!$ do j=1,NGLLY
+!!$ do i=1,NGLLX
+!!$ iglob = ibool(i,j,k,ispec)
+!!$
+!!$ potential_dot_dot_acoustic(iglob) = 0.d0
+!!$ potential_dot_acoustic(iglob) = 0.d0
+!!$ potential_acoustic(iglob) = 0.d0
+!!$ enddo
+!!$ enddo
+!!$ enddo
+!!$
+!!$ endif ! if( PML_CONDITIONS )
- i = 1
-
- do k=1,NGLLZ
- do j=1,NGLLY
- iglob = ibool(i,j,k,ispec)
-
- potential_dot_dot_acoustic(iglob) = 0.d0
- potential_dot_acoustic(iglob) = 0.d0
- potential_acoustic(iglob) = 0.d0
- enddo
- enddo
- enddo
-
- ! xmax
- do ispec2D=1,nspec2D_xmax
- ispec = ibelm_xmax(ispec2D)
-
- i = NGLLX
-
- do k=1,NGLLZ
- do j=1,NGLLY
- iglob = ibool(i,j,k,ispec)
-
- potential_dot_dot_acoustic(iglob) = 0.d0
- potential_dot_acoustic(iglob) = 0.d0
- potential_acoustic(iglob) = 0.d0
- enddo
- enddo
- enddo
-
- ! ymin
- do ispec2D=1,nspec2D_ymin
- ispec = ibelm_ymin(ispec2D)
-
- j = 1
-
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- potential_dot_dot_acoustic(iglob) = 0.d0
- potential_dot_acoustic(iglob) = 0.d0
- potential_acoustic(iglob) = 0.d0
- enddo
- enddo
- enddo
-
- ! ymax
- do ispec2D=1,nspec2D_ymax
- ispec = ibelm_ymax(ispec2D)
-
- j = NGLLY
-
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- potential_dot_dot_acoustic(iglob) = 0.d0
- potential_dot_acoustic(iglob) = 0.d0
- potential_acoustic(iglob) = 0.d0
- enddo
- enddo
- enddo
-
- ! bottom (zmin)
- do ispec2D=1,NSPEC2D_BOTTOM
- ispec = ibelm_bottom(ispec2D)
-
- k = 1
-
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- potential_dot_dot_acoustic(iglob) = 0.d0
- potential_dot_acoustic(iglob) = 0.d0
- potential_acoustic(iglob) = 0.d0
- enddo
- enddo
- enddo
-
- ! top (zmax)
- do ispec2D=1,NSPEC2D_BOTTOM
- ispec = ibelm_top(ispec2D)
-
- k = NGLLZ
-
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
-
- potential_dot_dot_acoustic(iglob) = 0.d0
- potential_dot_acoustic(iglob) = 0.d0
- potential_acoustic(iglob) = 0.d0
- enddo
- enddo
- enddo
-
- endif ! if( PML_CONDITIONS )
-
end subroutine compute_forces_acoustic_noDev
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 2013-03-24 18:58:31 UTC (rev 21622)
@@ -90,7 +90,7 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic)
+ phase_ispec_inner_elastic,ispec_is_elastic)
! adjoint simulations: backward/reconstructed wavefield
if( SIMULATION_TYPE == 3 ) &
@@ -122,7 +122,7 @@
b_dsdx_top,b_dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic)
+ phase_ispec_inner_elastic,ispec_is_elastic)
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_viscoelastic_noDev.f90 2013-03-24 18:58:31 UTC (rev 21622)
@@ -52,7 +52,7 @@
dsdx_top,dsdx_bot, &
ispec2D_moho_top,ispec2D_moho_bot, &
num_phase_ispec_elastic,nspec_inner_elastic,nspec_outer_elastic, &
- phase_ispec_inner_elastic)
+ phase_ispec_inner_elastic,ispec_is_elastic)
use constants, only: NGLLX,NGLLY,NGLLZ,NDIM,N_SLS,SAVE_MOHO_MESH,ONE_THIRD,FOUR_THIRDS,FULL_ATTENUATION_SOLID,IOUT !ZN
use pml_par
@@ -135,6 +135,8 @@
integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
integer, dimension(NSPEC2D_TOP) :: ibelm_top
+ logical, dimension(NSPEC_AB) :: ispec_is_elastic
+
! local parameters
integer :: i_SLS,imodulo_N_SLS
integer :: ispec,ispec2D,iglob,ispec_p,num_elements
@@ -181,9 +183,6 @@
tempx1_att,tempx2_att,tempx3_att,tempy1_att,tempy2_att,tempy3_att,tempz1_att,tempz2_att,tempz3_att
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc_att,dummyy_loc_att,dummyz_loc_att
- real(kind=CUSTOM_REAL) :: duxdxl_new,duxdyl_new,duxdzl_new,duydxl_new
- real(kind=CUSTOM_REAL) :: duydyl_new,duydzl_new,duzdxl_new,duzdyl_new,duzdzl_new;
-
real(kind=CUSTOM_REAL) :: eta
! local C-PML absorbing boundary conditions parameters
@@ -254,6 +253,21 @@
enddo
enddo
enddo
+ else if(PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc_att(i,j,k) = deltat*veloc(1,iglob)
+ dummyy_loc_att(i,j,k) = deltat*veloc(2,iglob)
+ dummyz_loc_att(i,j,k) = deltat*veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+ endif
endif
do k=1,NGLLZ
@@ -291,8 +305,7 @@
tempz3(i,j,k) = tempz3(i,j,k) + dummyz_loc(i,j,l)*hp3
enddo
- if( (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) .or. &
- (PML_CONDITIONS .and. CPML_mask_ibool(ispec)) ) then
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
tempx1_att(i,j,k) = tempx1(i,j,k)
tempx2_att(i,j,k) = tempx2(i,j,k)
tempx3_att(i,j,k) = tempx3(i,j,k)
@@ -308,24 +321,61 @@
! use first order Taylor expansion of displacement for local storage of stresses
! at this current time step, to fix attenuation in a consistent way
do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- tempx1_att(i,j,k) = tempx1_att(i,j,k) + dummyx_loc_att(l,j,k)*hp1
- tempy1_att(i,j,k) = tempy1_att(i,j,k) + dummyy_loc_att(l,j,k)*hp1
- tempz1_att(i,j,k) = tempz1_att(i,j,k) + dummyz_loc_att(l,j,k)*hp1
+ hp1 = hprime_xx(i,l)
+ tempx1_att(i,j,k) = tempx1_att(i,j,k) + dummyx_loc_att(l,j,k)*hp1
+ tempy1_att(i,j,k) = tempy1_att(i,j,k) + dummyy_loc_att(l,j,k)*hp1
+ tempz1_att(i,j,k) = tempz1_att(i,j,k) + dummyz_loc_att(l,j,k)*hp1
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ
- hp2 = hprime_yy(j,l)
- tempx2_att(i,j,k) = tempx2_att(i,j,k) + dummyx_loc_att(i,l,k)*hp2
- tempy2_att(i,j,k) = tempy2_att(i,j,k) + dummyy_loc_att(i,l,k)*hp2
- tempz2_att(i,j,k) = tempz2_att(i,j,k) + dummyz_loc_att(i,l,k)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp2 = hprime_yy(j,l)
+ tempx2_att(i,j,k) = tempx2_att(i,j,k) + dummyx_loc_att(i,l,k)*hp2
+ tempy2_att(i,j,k) = tempy2_att(i,j,k) + dummyy_loc_att(i,l,k)*hp2
+ tempz2_att(i,j,k) = tempz2_att(i,j,k) + dummyz_loc_att(i,l,k)*hp2
- !!! can merge these loops because NGLLX = NGLLY = NGLLZ
- hp3 = hprime_zz(k,l)
- tempx3_att(i,j,k) = tempx3_att(i,j,k) + dummyx_loc_att(i,j,l)*hp3
- tempy3_att(i,j,k) = tempy3_att(i,j,k) + dummyy_loc_att(i,j,l)*hp3
- tempz3_att(i,j,k) = tempz3_att(i,j,k) + dummyz_loc_att(i,j,l)*hp3
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp3 = hprime_zz(k,l)
+ tempx3_att(i,j,k) = tempx3_att(i,j,k) + dummyx_loc_att(i,j,l)*hp3
+ tempy3_att(i,j,k) = tempy3_att(i,j,k) + dummyy_loc_att(i,j,l)*hp3
+ tempz3_att(i,j,k) = tempz3_att(i,j,k) + dummyz_loc_att(i,j,l)*hp3
enddo
+ else if(PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ tempx1_att(i,j,k) = tempx1(i,j,k)
+ tempx2_att(i,j,k) = tempx2(i,j,k)
+ tempx3_att(i,j,k) = tempx3(i,j,k)
+
+ tempy1_att(i,j,k) = tempy1(i,j,k)
+ tempy2_att(i,j,k) = tempy2(i,j,k)
+ tempy3_att(i,j,k) = tempy3(i,j,k)
+
+ tempz1_att(i,j,k) = tempz1(i,j,k)
+ tempz2_att(i,j,k) = tempz2(i,j,k)
+ tempz3_att(i,j,k) = tempz3(i,j,k)
+
+ ! use first order Taylor expansion of displacement for local storage of stresses
+ ! at this current time step, to fix attenuation in a consistent way
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ tempx1_att(i,j,k) = tempx1_att(i,j,k) + dummyx_loc_att(l,j,k)*hp1
+ tempy1_att(i,j,k) = tempy1_att(i,j,k) + dummyy_loc_att(l,j,k)*hp1
+ tempz1_att(i,j,k) = tempz1_att(i,j,k) + dummyz_loc_att(l,j,k)*hp1
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp2 = hprime_yy(j,l)
+ tempx2_att(i,j,k) = tempx2_att(i,j,k) + dummyx_loc_att(i,l,k)*hp2
+ tempy2_att(i,j,k) = tempy2_att(i,j,k) + dummyy_loc_att(i,l,k)*hp2
+ tempz2_att(i,j,k) = tempz2_att(i,j,k) + dummyz_loc_att(i,l,k)*hp2
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ
+ hp3 = hprime_zz(k,l)
+ tempx3_att(i,j,k) = tempx3_att(i,j,k) + dummyx_loc_att(i,j,l)*hp3
+ tempy3_att(i,j,k) = tempy3_att(i,j,k) + dummyy_loc_att(i,j,l)*hp3
+ tempz3_att(i,j,k) = tempz3_att(i,j,k) + dummyz_loc_att(i,j,l)*hp3
+ enddo
+ endif
endif
! get derivatives of ux, uy and uz with respect to x, y and z
@@ -352,22 +402,26 @@
duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
- ! stores derivatives of ux, uy and uz with respect to x, y and z
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- ispec_CPML = spec_to_CPML(ispec)
+ ! stores derivatives of ux, uy and uz with respect to x, y and z
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ ispec_CPML = spec_to_CPML(ispec)
- PML_dux_dxl(i,j,k,ispec_CPML) = duxdxl
- PML_dux_dyl(i,j,k,ispec_CPML) = duxdyl
- PML_dux_dzl(i,j,k,ispec_CPML) = duxdzl
+ PML_dux_dxl(i,j,k,ispec_CPML) = duxdxl
+ PML_dux_dyl(i,j,k,ispec_CPML) = duxdyl
+ PML_dux_dzl(i,j,k,ispec_CPML) = duxdzl
- PML_duy_dxl(i,j,k,ispec_CPML) = duydxl
- PML_duy_dyl(i,j,k,ispec_CPML) = duydyl
- PML_duy_dzl(i,j,k,ispec_CPML) = duydzl
+ PML_duy_dxl(i,j,k,ispec_CPML) = duydxl
+ PML_duy_dyl(i,j,k,ispec_CPML) = duydyl
+ PML_duy_dzl(i,j,k,ispec_CPML) = duydzl
- PML_duz_dxl(i,j,k,ispec_CPML) = duzdxl
- PML_duz_dyl(i,j,k,ispec_CPML) = duzdyl
- PML_duz_dzl(i,j,k,ispec_CPML) = duzdzl
- endif
+ PML_duz_dxl(i,j,k,ispec_CPML) = duzdxl
+ PML_duz_dyl(i,j,k,ispec_CPML) = duzdyl
+ PML_duz_dzl(i,j,k,ispec_CPML) = duzdzl
+ endif
+ endif
! save strain on the Moho boundary
if (SAVE_MOHO_MESH ) then
@@ -402,9 +456,7 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- if( (ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) .or. &
- (PML_CONDITIONS .and. CPML_mask_ibool(ispec)) ) then
-
+ if(ATTENUATION .and. COMPUTE_AND_STORE_STRAIN) then
! temporary variables used for fixing attenuation in a consistent way
duxdxl_att = xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
duxdyl_att = xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
@@ -432,37 +484,33 @@
epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl_att
epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl_att
epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl_att
+
+ else if(PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ PML_dux_dxl_new(i,j,k,ispec_CPML) = &
+ xixl*tempx1_att(i,j,k) + etaxl*tempx2_att(i,j,k) + gammaxl*tempx3_att(i,j,k)
+ PML_dux_dyl_new(i,j,k,ispec_CPML) = &
+ xiyl*tempx1_att(i,j,k) + etayl*tempx2_att(i,j,k) + gammayl*tempx3_att(i,j,k)
+ PML_dux_dzl_new(i,j,k,ispec_CPML) = &
+ xizl*tempx1_att(i,j,k) + etazl*tempx2_att(i,j,k) + gammazl*tempx3_att(i,j,k)
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
-!! DK DK to Jo and Zhinan to debug CPML: duxdxl_new is undefined
-!! DK DK to Jo and Zhinan to debug CPML: maybe replace it with "duxdxl_att"??
-!! DK DK (this is a variable name that has changed in this routine since Jo left, thus it cannot explain the previous CPML bug)
-! PML_dux_dxl_new(i,j,k,ispec_CPML) = duxdxl_new
-! PML_dux_dyl_new(i,j,k,ispec_CPML) = duxdyl_new
-! PML_dux_dzl_new(i,j,k,ispec_CPML) = duxdzl_new
+ PML_duy_dxl_new(i,j,k,ispec_CPML) = &
+ xixl*tempy1_att(i,j,k) + etaxl*tempy2_att(i,j,k) + gammaxl*tempy3_att(i,j,k)
+ PML_duy_dyl_new(i,j,k,ispec_CPML) = &
+ xiyl*tempy1_att(i,j,k) + etayl*tempy2_att(i,j,k) + gammayl*tempy3_att(i,j,k)
+ PML_duy_dzl_new(i,j,k,ispec_CPML) = &
+ xizl*tempy1_att(i,j,k) + etazl*tempy2_att(i,j,k) + gammazl*tempy3_att(i,j,k)
-! PML_duy_dxl_new(i,j,k,ispec_CPML) = duydxl_new
-! PML_duy_dyl_new(i,j,k,ispec_CPML) = duydyl_new
-! PML_duy_dzl_new(i,j,k,ispec_CPML) = duydzl_new
+ PML_duz_dxl_new(i,j,k,ispec_CPML) = &
+ xixl*tempz1_att(i,j,k) + etaxl*tempz2_att(i,j,k) + gammaxl*tempz3_att(i,j,k)
+ PML_duz_dyl_new(i,j,k,ispec_CPML) = &
+ xiyl*tempz1_att(i,j,k) + etayl*tempz2_att(i,j,k) + gammayl*tempz3_att(i,j,k)
+ PML_duz_dzl_new(i,j,k,ispec_CPML) = &
+ xizl*tempz1_att(i,j,k) + etazl*tempz2_att(i,j,k) + gammazl*tempz3_att(i,j,k)
+ endif
-! PML_duz_dxl_new(i,j,k,ispec_CPML) = duzdxl_new
-! PML_duz_dyl_new(i,j,k,ispec_CPML) = duzdyl_new
-! PML_duz_dzl_new(i,j,k,ispec_CPML) = duzdzl_new
-
-!! DK DK 23 March 2013: I make the change from "duxdxl_new" to "duxdxl_att" here, see my comment above
- PML_dux_dxl_new(i,j,k,ispec_CPML) = duxdxl_att
- PML_dux_dyl_new(i,j,k,ispec_CPML) = duxdyl_att
- PML_dux_dzl_new(i,j,k,ispec_CPML) = duxdzl_att
-
- PML_duy_dxl_new(i,j,k,ispec_CPML) = duydxl_att
- PML_duy_dyl_new(i,j,k,ispec_CPML) = duydyl_att
- PML_duy_dzl_new(i,j,k,ispec_CPML) = duydzl_att
-
- PML_duz_dxl_new(i,j,k,ispec_CPML) = duzdxl_att
- PML_duz_dyl_new(i,j,k,ispec_CPML) = duzdyl_att
- PML_duz_dzl_new(i,j,k,ispec_CPML) = duzdzl_att
- endif
-
else
! computes deviatoric strain attenuation and/or for kernel calculations
if (COMPUTE_AND_STORE_STRAIN) then
@@ -631,41 +679,64 @@
!! DK DK that when PML_CONDITIONS is on then you do not compute the tempx, tempy, tempz arrays
!! DK DK (even in non-PML elements!!), even though such arrays are needed below;
!! DK DK shouldn't there be at least a "if (is_CPML(ispec))" test as well here, or something like that?
- if( .not. PML_CONDITIONS ) then
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(.not.CPML_mask_ibool(ispec)) then
+ ! define symmetric components of sigma
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
- ! define symmetric components of sigma
- sigma_yx = sigma_xy
- sigma_zx = sigma_xz
- sigma_zy = sigma_yz
+ ! form dot product with test vector, non-symmetric form (which is useful in the case of PML)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- ! form dot product with test vector, non-symmetric form (which is useful in the case of PML)
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ endif
+ else
+ ! define symmetric components of sigma
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ ! form dot product with test vector, non-symmetric form (which is useful in the case of PML)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- endif
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
- enddo
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ endif
+ enddo
enddo
enddo
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- ! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
- call pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
- tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
- sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,xixl,xiyl,xizl, &
- etaxl,etayl,etazl,gammaxl,gammayl,gammazl)
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ ! sets C-PML elastic memory variables to compute stress sigma and form dot product with test vector
+ call pml_compute_memory_variables(ispec,ispec_CPML,deltat,jacobianl,tempx1,tempy1,tempz1,tempx2,tempy2,tempz2, &
+ tempx3,tempy3,tempz3,sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz, &
+ sigma_yx,sigma_zx,sigma_zy,lambdal,mul,lambdalplus2mul,xixl,xiyl,xizl, &
+ etaxl,etayl,etazl,gammaxl,gammayl,gammazl)
- ! calculates contribution from each C-PML element to update acceleration
- call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+ ! calculates contribution from each C-PML element to update acceleration
+ call pml_compute_accel_contribution(ispec,ispec_CPML,deltat,jacobianl,accel_elastic_CPML)
+ endif
endif
! second double-loop over GLL to compute all the terms
@@ -718,64 +789,68 @@
fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
! updates acceleration with contribution from each C-PML element
- if( PML_CONDITIONS .and. CPML_mask_ibool(ispec) ) then
- accel(1,iglob) = accel(1,iglob) - accel_elastic_CPML(1,i,j,k,ispec_CPML)
- accel(2,iglob) = accel(2,iglob) - accel_elastic_CPML(2,i,j,k,ispec_CPML)
- accel(3,iglob) = accel(3,iglob) - accel_elastic_CPML(3,i,j,k,ispec_CPML)
+ if (PML_CONDITIONS) then
+ ! do not merge this second line with the first using an ".and." statement
+ ! because array CPML_mask_ibool() is unallocated when PML_CONDITIONS is false
+ if(CPML_mask_ibool(ispec)) then
+ accel(1,iglob) = accel(1,iglob) - accel_elastic_CPML(1,i,j,k,ispec_CPML)
+ accel(2,iglob) = accel(2,iglob) - accel_elastic_CPML(2,i,j,k,ispec_CPML)
+ accel(3,iglob) = accel(3,iglob) - accel_elastic_CPML(3,i,j,k,ispec_CPML)
+ endif
endif
- ! update memory variables based upon the Runge-Kutta scheme
- if(ATTENUATION) then
+ ! update memory variables based upon the Runge-Kutta scheme
+ if(ATTENUATION) then
- ! use Runge-Kutta scheme to march in time
- do i_sls = 1,N_SLS
+ ! use Runge-Kutta scheme to march in time
+ do i_sls = 1,N_SLS
- alphaval_loc = alphaval(i_sls)
- betaval_loc = betaval(i_sls)
- gammaval_loc = gammaval(i_sls)
+ alphaval_loc = alphaval(i_sls)
+ betaval_loc = betaval(i_sls)
+ gammaval_loc = gammaval(i_sls)
- if(FULL_ATTENUATION_SOLID) then
- ! term in trace !ZN
- factor_loc = kappastore(i,j,k,ispec) * factor_common_kappa(i_sls,i,j,k,ispec) !ZN
+ if(FULL_ATTENUATION_SOLID) then
+ ! term in trace !ZN
+ factor_loc = kappastore(i,j,k,ispec) * factor_common_kappa(i_sls,i,j,k,ispec) !ZN
- Sn = factor_loc * epsilondev_trace(i,j,k,ispec) !ZN
- Snp1 = factor_loc * epsilondev_trace_loc(i,j,k) !ZN
- R_trace(i,j,k,ispec,i_sls) = alphaval_loc * R_trace(i,j,k,ispec,i_sls) + & !ZN
- betaval_loc * Sn + gammaval_loc * Snp1 !ZN
- endif
+ Sn = factor_loc * epsilondev_trace(i,j,k,ispec) !ZN
+ Snp1 = factor_loc * epsilondev_trace_loc(i,j,k) !ZN
+ R_trace(i,j,k,ispec,i_sls) = alphaval_loc * R_trace(i,j,k,ispec,i_sls) + & !ZN
+ betaval_loc * Sn + gammaval_loc * Snp1 !ZN
+ endif
- ! term in xx yy zz xy xz yz
- factor_loc = mustore(i,j,k,ispec) * factor_common(i_sls,i,j,k,ispec)
+ ! term in xx yy zz xy xz yz
+ factor_loc = mustore(i,j,k,ispec) * factor_common(i_sls,i,j,k,ispec)
- ! term in xx
- Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
- R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
- ! term in yy
- Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
- R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
- ! term in zz not computed since zero trace
- ! term in xy
- Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
- R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
- ! term in xz
- Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
- R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
- ! term in yz
- Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
- R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + &
- betaval_loc * Sn + gammaval_loc * Snp1
- enddo ! end of loop on memory variables
+ ! term in xx
+ Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
+ R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+ ! term in yy
+ Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
+ R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+ ! term in zz not computed since zero trace
+ ! term in xy
+ Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
+ R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+ ! term in xz
+ Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
+ R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+ ! term in yz
+ Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
+ R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + &
+ betaval_loc * Sn + gammaval_loc * Snp1
+ enddo ! end of loop on memory variables
- endif ! end of if attenuation
+ endif ! end of if attenuation
enddo
enddo
@@ -807,150 +882,162 @@
do ispec2D=1,nspec2D_xmin
ispec = ibelm_xmin(ispec2D)
- i = 1
+ if(CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec)) then
+ i = 1
- do k=1,NGLLZ
- do j=1,NGLLY
- iglob = ibool(i,j,k,ispec)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ accel(1,iglob) = 0._CUSTOM_REAL
+ accel(2,iglob) = 0._CUSTOM_REAL
+ accel(3,iglob) = 0._CUSTOM_REAL
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ veloc(1,iglob) = 0._CUSTOM_REAL
+ veloc(2,iglob) = 0._CUSTOM_REAL
+ veloc(3,iglob) = 0._CUSTOM_REAL
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
+ displ(1,iglob) = 0._CUSTOM_REAL
+ displ(2,iglob) = 0._CUSTOM_REAL
+ displ(3,iglob) = 0._CUSTOM_REAL
+ enddo
enddo
- enddo
+ endif
enddo
! xmax
do ispec2D=1,nspec2D_xmax
ispec = ibelm_xmax(ispec2D)
- i = NGLLX
+ if(CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec)) then
+ i = NGLLX
- do k=1,NGLLZ
- do j=1,NGLLY
- iglob = ibool(i,j,k,ispec)
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ accel(1,iglob) = 0._CUSTOM_REAL
+ accel(2,iglob) = 0._CUSTOM_REAL
+ accel(3,iglob) = 0._CUSTOM_REAL
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ veloc(1,iglob) = 0._CUSTOM_REAL
+ veloc(2,iglob) = 0._CUSTOM_REAL
+ veloc(3,iglob) = 0._CUSTOM_REAL
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
+ displ(1,iglob) = 0._CUSTOM_REAL
+ displ(2,iglob) = 0._CUSTOM_REAL
+ displ(3,iglob) = 0._CUSTOM_REAL
+ enddo
enddo
- enddo
+ endif
enddo
! ymin
do ispec2D=1,nspec2D_ymin
ispec = ibelm_ymin(ispec2D)
- j = 1
+ if(CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec)) then
+ j = 1
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
+ do k=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ accel(1,iglob) = 0._CUSTOM_REAL
+ accel(2,iglob) = 0._CUSTOM_REAL
+ accel(3,iglob) = 0._CUSTOM_REAL
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ veloc(1,iglob) = 0._CUSTOM_REAL
+ veloc(2,iglob) = 0._CUSTOM_REAL
+ veloc(3,iglob) = 0._CUSTOM_REAL
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
+ displ(1,iglob) = 0._CUSTOM_REAL
+ displ(2,iglob) = 0._CUSTOM_REAL
+ displ(3,iglob) = 0._CUSTOM_REAL
+ enddo
enddo
- enddo
+ endif
enddo
! ymax
do ispec2D=1,nspec2D_ymax
ispec = ibelm_ymax(ispec2D)
- j = NGLLY
+ if(CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec)) then
+ j = NGLLY
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
+ do k=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ accel(1,iglob) = 0._CUSTOM_REAL
+ accel(2,iglob) = 0._CUSTOM_REAL
+ accel(3,iglob) = 0._CUSTOM_REAL
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ veloc(1,iglob) = 0._CUSTOM_REAL
+ veloc(2,iglob) = 0._CUSTOM_REAL
+ veloc(3,iglob) = 0._CUSTOM_REAL
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
+ displ(1,iglob) = 0._CUSTOM_REAL
+ displ(2,iglob) = 0._CUSTOM_REAL
+ displ(3,iglob) = 0._CUSTOM_REAL
+ enddo
enddo
- enddo
+ endif
enddo
! bottom (zmin)
do ispec2D=1,NSPEC2D_BOTTOM
ispec = ibelm_bottom(ispec2D)
- k = 1
+ if(CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec)) then
+ k = 1
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ accel(1,iglob) = 0._CUSTOM_REAL
+ accel(2,iglob) = 0._CUSTOM_REAL
+ accel(3,iglob) = 0._CUSTOM_REAL
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ veloc(1,iglob) = 0._CUSTOM_REAL
+ veloc(2,iglob) = 0._CUSTOM_REAL
+ veloc(3,iglob) = 0._CUSTOM_REAL
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
+ displ(1,iglob) = 0._CUSTOM_REAL
+ displ(2,iglob) = 0._CUSTOM_REAL
+ displ(3,iglob) = 0._CUSTOM_REAL
+ enddo
enddo
- enddo
+ endif
enddo
! top (zmax)
do ispec2D=1,NSPEC2D_BOTTOM
ispec = ibelm_top(ispec2D)
- k = NGLLZ
+ if(CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec)) then
+ k = NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec)
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
- accel(1,iglob) = 0._CUSTOM_REAL
- accel(2,iglob) = 0._CUSTOM_REAL
- accel(3,iglob) = 0._CUSTOM_REAL
+ accel(1,iglob) = 0._CUSTOM_REAL
+ accel(2,iglob) = 0._CUSTOM_REAL
+ accel(3,iglob) = 0._CUSTOM_REAL
- veloc(1,iglob) = 0._CUSTOM_REAL
- veloc(2,iglob) = 0._CUSTOM_REAL
- veloc(3,iglob) = 0._CUSTOM_REAL
+ veloc(1,iglob) = 0._CUSTOM_REAL
+ veloc(2,iglob) = 0._CUSTOM_REAL
+ veloc(3,iglob) = 0._CUSTOM_REAL
- displ(1,iglob) = 0._CUSTOM_REAL
- displ(2,iglob) = 0._CUSTOM_REAL
- displ(3,iglob) = 0._CUSTOM_REAL
+ displ(1,iglob) = 0._CUSTOM_REAL
+ displ(2,iglob) = 0._CUSTOM_REAL
+ displ(3,iglob) = 0._CUSTOM_REAL
+ enddo
enddo
- enddo
+ endif
enddo
endif ! if( PML_CONDITIONS )
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-03-24 00:54:00 UTC (rev 21621)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_allocate_arrays.f90 2013-03-24 18:58:31 UTC (rev 21622)
@@ -158,7 +158,7 @@
! stores C-PML memory variables needed for potential
allocate(rmemory_potential_acoustic(NGLLX,NGLLY,NGLLZ,NSPEC_CPML,3),stat=ier)
- if(ier /= 0) stop 'error allocating rmemory_displ_elastic array'
+ if(ier /= 0) stop 'error allocating rmemory_potential_acoustic array'
! stores C-PML contribution to update acceleration to the global mesh
allocate(accel_elastic_CPML(NDIM,NGLLX,NGLLY,NGLLZ,NSPEC_CPML),stat=ier)
More information about the CIG-COMMITS
mailing list