[cig-commits] r21625 - in seismo/3D/SPECFEM3D/trunk: examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE examples/CPML/purely_elastic/HOMO8_NOFREESURFACE src/generate_databases src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Sun Mar 24 15:29:45 PDT 2013
Author: joseph.charles
Date: 2013-03-24 15:29:45 -0700 (Sun, 24 Mar 2013)
New Revision: 21625
Added:
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/process.sh
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/process.sh
Removed:
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt
Modified:
seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
Log:
fixes bugs in acoustic C-PML simulations; create bash scripts for both C-PML examples
Deleted: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt 2013-03-24 22:10:15 UTC (rev 21624)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/HOWTO_run_this_example.txt 2013-03-24 22:29:45 UTC (rev 21625)
@@ -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_acoustic/HOMO8_FREE_SURFACE/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/process.sh (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/process.sh 2013-03-24 22:29:45 UTC (rev 21625)
@@ -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`
+
+
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_FREE_SURFACE/process.sh
___________________________________________________________________
Name: svn:executable
+ *
Deleted: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt 2013-03-24 22:10:15 UTC (rev 21624)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/HOWTO_run_this_example.txt 2013-03-24 22:29:45 UTC (rev 21625)
@@ -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_acoustic/HOMO8_NOFREESURFACE/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/process.sh (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/process.sh 2013-03-24 22:29:45 UTC (rev 21625)
@@ -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`
+
+
Property changes on: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_acoustic/HOMO8_NOFREESURFACE/process.sh
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh 2013-03-24 22:10:15 UTC (rev 21624)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/purely_elastic/HOMO8_NOFREESURFACE/process.sh 2013-03-24 22:29:45 UTC (rev 21625)
@@ -49,7 +49,7 @@
# links executables
cd bin/
-rm -f bin/*
+rm -f *
ln -s ../../../../../bin/xdecompose_mesh .
ln -s ../../../../../bin/xgenerate_databases .
ln -s ../../../../../bin/xspecfem3D .
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-03-24 22:10:15 UTC (rev 21624)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2013-03-24 22:29:45 UTC (rev 21625)
@@ -450,16 +450,12 @@
! loops over physical mesh elements
do ispec=1,nspec
- if( .not. CPML_mask_ibool(ispec) ) then
+ if( .not. CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec) ) then
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -476,6 +472,28 @@
enddo
enddo
enddo
+ else if( .not. CPML_mask_ibool(ispec) .and. ispec_is_acoustic(ispec) ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef
+ endif
+ enddo
+ enddo
+ enddo
endif
enddo
@@ -483,18 +501,14 @@
do ispec_CPML=1,nspec_cpml
ispec = CPML_to_spec(ispec_CPML)
- if( CPML_mask_ibool(ispec) ) then
+ if( CPML_mask_ibool(ispec) .and. ispec_is_elastic(ispec) ) then
! X_surface C-PML
if( CPML_regions(ispec_CPML) == 1 ) then
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -520,11 +534,7 @@
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -550,11 +560,7 @@
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -580,11 +586,7 @@
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -614,11 +616,7 @@
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -648,11 +646,7 @@
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -682,11 +676,7 @@
do j=1,NGLLY
do i=1,NGLLX
! defines the material coefficient associated to the domain
- if( ispec_is_elastic(ispec) ) then
- mat_coef = rhostore(i,j,k,ispec)
- elseif( ispec_is_acoustic(ispec) ) then
- mat_coef = 1.d0 / kappastore(i,j,k,ispec)
- endif
+ mat_coef = rhostore(i,j,k,ispec)
iglob = ibool(i,j,k,ispec)
@@ -719,4 +709,216 @@
endif
enddo ! do ispec_CPML=1,nspec_cpml
+ ! loops over C-PML elements
+ do ispec_CPML=1,nspec_cpml
+ ispec = CPML_to_spec(ispec_CPML)
+
+ if( CPML_mask_ibool(ispec) .and. ispec_is_acoustic(ispec) ) then
+ ! X_surface C-PML
+ if( CPML_regions(ispec_CPML) == 1 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_x(i,j,k,ispec_CPML)) + dble(d_store_x(i,j,k,ispec_CPML)) * deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_x(i,j,k,ispec_CPML) + d_store_x(i,j,k,ispec_CPML) * deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! Y_surface C-PML
+ elseif( CPML_regions(ispec_CPML) == 2 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_y(i,j,k,ispec_CPML)) + dble(d_store_y(i,j,k,ispec_CPML)) * deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_y(i,j,k,ispec_CPML) + d_store_y(i,j,k,ispec_CPML) * deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! Z_surface C-PML
+ elseif( CPML_regions(ispec_CPML) == 3 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_z(i,j,k,ispec_CPML)) + dble(d_store_z(i,j,k,ispec_CPML)) * deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_z(i,j,k,ispec_CPML) + d_store_z(i,j,k,ispec_CPML) * deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! XY_edge C-PML
+ elseif( CPML_regions(ispec_CPML) == 4 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_x(i,j,k,ispec_CPML)) * dble(K_store_y(i,j,k,ispec_CPML)) + &
+ (dble(d_store_x(i,j,k,ispec_CPML)) * dble(K_store_y(i,j,k,ispec_CPML)) + &
+ dble(d_store_y(i,j,k,ispec_CPML)) * dble(K_store_x(i,j,k,ispec_CPML))) * deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_x(i,j,k,ispec_CPML) * K_store_y(i,j,k,ispec_CPML) + &
+ (d_store_x(i,j,k,ispec_CPML) * K_store_y(i,j,k,ispec_CPML) + &
+ d_store_y(i,j,k,ispec_CPML) * K_store_x(i,j,k,ispec_CPML)) * deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! XZ_edge C-PML
+ elseif( CPML_regions(ispec_CPML) == 5 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_x(i,j,k,ispec_CPML)) * dble(K_store_z(i,j,k,ispec_CPML)) + &
+ (dble(d_store_x(i,j,k,ispec_CPML)) * dble(K_store_z(i,j,k,ispec_CPML)) + &
+ dble(d_store_z(i,j,k,ispec_CPML)) * dble(K_store_x(i,j,k,ispec_CPML))) * deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_x(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ (d_store_x(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ d_store_z(i,j,k,ispec_CPML) * K_store_x(i,j,k,ispec_CPML)) * deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! YZ_edge C-PML
+ elseif( CPML_regions(ispec_CPML) == 6 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_y(i,j,k,ispec_CPML)) * dble(K_store_z(i,j,k,ispec_CPML)) + &
+ (dble(d_store_y(i,j,k,ispec_CPML)) * dble(K_store_z(i,j,k,ispec_CPML)) + &
+ dble(d_store_z(i,j,k,ispec_CPML)) * dble(K_store_y(i,j,k,ispec_CPML))) * deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_y(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ (d_store_y(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ d_store_z(i,j,k,ispec_CPML) * K_store_y(i,j,k,ispec_CPML)) * deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+
+ ! XYZ_corner C-PML
+ elseif( CPML_regions(ispec_CPML) == 7 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! defines the material coefficient associated to the domain
+ mat_coef = 1.d0 / kappastore(i,j,k,ispec)
+
+ iglob = ibool(i,j,k,ispec)
+
+ weight = wxgll(i)*wygll(j)*wzgll(k)
+ jacobianl = jacobianstore(i,j,k,ispec)
+
+ if(CUSTOM_REAL == SIZE_REAL) then
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ sngl( dble(jacobianl) * weight * dble(mat_coef) * &
+ (dble(K_store_x(i,j,k,ispec_CPML)) * dble(K_store_y(i,j,k,ispec_CPML)) * &
+ dble(K_store_z(i,j,k,ispec_CPML)) + (dble(d_store_x(i,j,k,ispec_CPML)) * &
+ dble(K_store_y(i,j,k,ispec_CPML)) * dble(K_store_z(i,j,k,ispec_CPML)) + &
+ dble(d_store_y(i,j,k,ispec_CPML)) * dble(K_store_x(i,j,k,ispec_CPML)) * &
+ dble(K_store_z(i,j,k,ispec_CPML)) + dble(d_store_z(i,j,k,ispec_CPML)) * &
+ dble(K_store_y(i,j,k,ispec_CPML)) * dble(K_store_z(i,j,k,ispec_CPML))) * &
+ deltat / 2.d0) )
+ else
+ rmass_acoustic(iglob) = rmass_acoustic(iglob) + &
+ jacobianl * weight * mat_coef * &
+ (K_store_x(i,j,k,ispec_CPML) * K_store_y(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ (d_store_x(i,j,k,ispec_CPML) * K_store_y(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ d_store_y(i,j,k,ispec_CPML) * K_store_x(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML) + &
+ d_store_z(i,j,k,ispec_CPML) * K_store_y(i,j,k,ispec_CPML) * K_store_z(i,j,k,ispec_CPML)) * &
+ deltat / 2.d0)
+ endif
+ enddo
+ enddo
+ enddo
+ endif
+ endif
+ enddo ! do ispec_CPML=1,nspec_cpml
+
end subroutine create_mass_matrices_pml
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-03-24 22:10:15 UTC (rev 21624)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_accel_contribution.f90 2013-03-24 22:29:45 UTC (rev 21625)
@@ -61,8 +61,11 @@
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
- rhol = rho_vp(i,j,k,ispec)
- kappal = kappastore(i,j,k,ispec)
+ if( ispec_is_elastic(ispec) ) then
+ rhol = rho_vp(i,j,k,ispec)
+ else if( ispec_is_acoustic(ispec) ) then
+ kappal = kappastore(i,j,k,ispec)
+ endif
iglob = ibool(i,j,k,ispec)
More information about the CIG-COMMITS
mailing list