[cig-commits] [commit] devel,master: adds configuration option --with-cem; adds file src/meshfem3D/model_cem.f90 for comprehensive earth model (8f81cfe)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:29:31 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 8f81cfeaeac461f4becc68841f66a95658c5a5d8
Author: daniel peter <peterda at ethz.ch>
Date: Thu Aug 28 14:40:24 2014 +0200
adds configuration option --with-cem; adds file src/meshfem3D/model_cem.f90 for comprehensive earth model
>---------------------------------------------------------------
8f81cfeaeac461f4becc68841f66a95658c5a5d8
Makefile.in | 14 ++
configure | 76 +++++++++-
configure.ac | 45 ++++++
src/meshfem3D/model_cem.f90 | 341 ++++++++++++++++++++++++++++++++++++++++++++
src/meshfem3D/rules.mk | 12 ++
5 files changed, 487 insertions(+), 1 deletion(-)
diff --git a/Makefile.in b/Makefile.in
index 7114287..a471fa4 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -169,6 +169,20 @@ MPILIBS += @ADIOS_LIBS@
#######################################
####
+#### CEM
+#### with configure: ./configure --with-cem CEM_LIBS=.. CEM_FCFLAGS=..
+####
+#######################################
+
+ at COND_CEM_TRUE@CEM = yes
+ at COND_CEM_FALSE@CEM = no
+
+FCFLAGS += @CEM_FCFLAGS@
+MPILIBS += @CEM_LIBS@
+
+
+#######################################
+####
#### directories
####
#######################################
diff --git a/configure b/configure
index 0486526..1056b4e 100755
--- a/configure
+++ b/configure
@@ -623,6 +623,8 @@ ac_includes_default="\
ac_subst_vars='LTLIBOBJS
LIBOBJS
+CEM_LIBS
+CEM_FCFLAGS
MPI_INCLUDES
MPI_INC
VTK_LIBS
@@ -678,6 +680,8 @@ ac_ct_FC
LDFLAGS
FCFLAGS
FC
+COND_CEM_FALSE
+COND_CEM_TRUE
COND_VECTORIZATION_FALSE
COND_VECTORIZATION_TRUE
COND_ADIOS_FALSE
@@ -747,6 +751,7 @@ with_opencl
enable_vtk
with_adios
enable_vectorization
+with_cem
with_vtk
with_vtk_suffix
with_vtk_version
@@ -778,7 +783,9 @@ OCL_CPU_FLAGS
OCL_GPU_FLAGS
OCL_INC
OCL_LIB
-MPI_INC'
+MPI_INC
+CEM_FCFLAGS
+CEM_LIBS'
# Initialize some variables set by options.
@@ -1404,6 +1411,7 @@ Optional Packages:
--with-cuda build cuda GPU enabled version [default=no]
--with-opencl build OpenCL GPU enabled version [default=no]
--with-adios build ADIOS enabled version [default=no]
+ --with-cem build CEM enabled version [default=no]
--with-vtk The prefix where VTK is installed [default=/usr]
--with-vtk-suffix Suffix to append to VTK's include directory, e.g.,
for vtk-5.2/ the suffix is "-5.2"
@@ -1443,6 +1451,8 @@ Some influential environment variables:
OCL_LIB Location of OpenCL library libOpenCL
MPI_INC Location of MPI include mpi.h, needed by CUDA/VTK (default:
auto-detect)
+ CEM_FCFLAGS CEM Fortran compiler flags
+ CEM_LIBS CEM extra libraries for linking programs
Use these variables to override the choices made by `configure' or to help
it to find libraries and programs with nonstandard names/locations.
@@ -2833,6 +2843,27 @@ else
fi
+###
+### CEM
+###
+
+
+# Check whether --with-cem was given.
+if test "${with_cem+set}" = set; then :
+ withval=$with_cem; want_cem="$withval"
+else
+ want_cem=no
+fi
+
+ if test x"$want_cem" != xno; then
+ COND_CEM_TRUE=
+ COND_CEM_FALSE='#'
+else
+ COND_CEM_TRUE='#'
+ COND_CEM_FALSE=
+fi
+
+
############################################################
# Checks for programs.
@@ -7456,6 +7487,45 @@ ac_compiler_gnu=$ac_cv_fc_compiler_gnu
fi
+
+###
+### CEM
+###
+
+if test x"$want_cem" != xno; then :
+
+ $as_echo "## ------------------------------- ##
+## Comprehensive Earth Model (CEM) ##
+## ------------------------------- ##"
+ { $as_echo "$as_me:${as_lineno-$LINENO}: CEM is enabled" >&5
+$as_echo "$as_me: CEM is enabled" >&6;}
+
+
+
+
+ # adds compiler flag
+ FLAGS_CHECK="${FLAGS_CHECK} ${FC_DEFINE}CEM"
+
+ # daniel: todo add NetCDF checking for Fortran, something like:
+ # netcdf checking
+ #AC_LANG_PUSH(Fortran)
+ #AC_MSG_CHECKING([for NetCDF modules])
+ #FCFLAGS="$CEM_FCFLAGS $FCFLAGS"
+ #AC_COMPILE_IFELSE([
+ # AC_LANG_PROGRAM([], [[
+ # use netcdf
+ # ]])
+ #], [
+ # AC_MSG_RESULT(yes)
+ #], [
+ # AC_MSG_RESULT(no)
+ # AC_MSG_ERROR([NetCDF modules not found; is NetCDF built with Fortran support for this compiler?])
+ #])
+ #AC_LANG_POP([Fortran])
+
+fi
+
+
############################################################
$as_echo "## ----------------------------------- ##
@@ -7609,6 +7679,10 @@ if test -z "${COND_VECTORIZATION_TRUE}" && test -z "${COND_VECTORIZATION_FALSE}"
as_fn_error $? "conditional \"COND_VECTORIZATION\" was never defined.
Usually this means the macro was only invoked conditionally." "$LINENO" 5
fi
+if test -z "${COND_CEM_TRUE}" && test -z "${COND_CEM_FALSE}"; then
+ as_fn_error $? "conditional \"COND_CEM\" was never defined.
+Usually this means the macro was only invoked conditionally." "$LINENO" 5
+fi
case $FC_MODINC in #(
*\ ) FC_MODINC=$FC_MODINC'${ac_empty}' ;;
esac
diff --git a/configure.ac b/configure.ac
index eb809dd..12d4334 100644
--- a/configure.ac
+++ b/configure.ac
@@ -94,6 +94,17 @@ AC_ARG_ENABLE([vectorization],
[want_vec=yes])
AM_CONDITIONAL([COND_VECTORIZATION], [test x"$want_vec" != xno])
+###
+### CEM
+###
+
+AC_ARG_WITH([cem],
+ [AC_HELP_STRING([--with-cem],
+ [build CEM enabled version @<:@default=no@:>@])],
+ [want_cem="$withval"],
+ [want_cem=no])
+AM_CONDITIONAL([COND_CEM], [test x"$want_cem" != xno])
+
############################################################
# Checks for programs.
@@ -237,6 +248,40 @@ AS_IF([test x"$want_cuda" != xno -o x"$want_vtk" != xno], [
CIT_MPI_INCDIR([$MPIFC])
])
+
+###
+### CEM
+###
+
+AS_IF([test x"$want_cem" != xno], [
+ AS_BOX([Comprehensive Earth Model (CEM)])
+ AC_MSG_NOTICE([CEM is enabled])
+
+ AC_ARG_VAR(CEM_FCFLAGS, [CEM Fortran compiler flags])
+ AC_ARG_VAR(CEM_LIBS, [CEM extra libraries for linking programs])
+
+ # adds compiler flag
+ FLAGS_CHECK="${FLAGS_CHECK} ${FC_DEFINE}CEM"
+
+ # daniel: todo add NetCDF checking for Fortran, something like:
+ # netcdf checking
+ #AC_LANG_PUSH(Fortran)
+ #AC_MSG_CHECKING([for NetCDF modules])
+ #FCFLAGS="$CEM_FCFLAGS $FCFLAGS"
+ #AC_COMPILE_IFELSE([
+ # AC_LANG_PROGRAM([], [[
+ # use netcdf
+ # ]])
+ #], [
+ # AC_MSG_RESULT(yes)
+ #], [
+ # AC_MSG_RESULT(no)
+ # AC_MSG_ERROR([NetCDF modules not found; is NetCDF built with Fortran support for this compiler?])
+ #])
+ #AC_LANG_POP([Fortran])
+])
+
+
############################################################
AS_BOX([setting up default simulation setup])
diff --git a/src/meshfem3D/model_cem.f90 b/src/meshfem3D/model_cem.f90
new file mode 100644
index 0000000..c4be778
--- /dev/null
+++ b/src/meshfem3D/model_cem.f90
@@ -0,0 +1,341 @@
+module cem_par
+
+ use :: constants, only: PI, GRAV, RHOAV, R_EARTH
+
+ implicit none
+
+ double precision, parameter :: scaleval=dsqrt(PI*GRAV*RHOAV)
+ double precision, parameter :: scale_GPa = (RHOAV / 1000.d0) * &
+ ((R_EARTH * scaleval / 1000.d0) ** 2)
+
+
+ integer, dimension (:), allocatable :: regCode
+ integer, parameter :: shuOn=1
+ integer, parameter :: comLvl=9
+ integer, parameter :: comOn=1
+ integer, parameter :: NDIMS=3
+ integer :: nnodes_cem
+ integer :: nelem_cem
+ integer :: rank
+
+ real, dimension (:,:), allocatable :: xyzOut
+
+ type par
+
+ double precision, dimension (:), allocatable :: vsv, vsh, vpv, vph, rho
+ double precision, dimension (:), allocatable :: c11, c12, c13, c14
+ double precision, dimension (:), allocatable :: c15, c16, c26, c33
+ double precision, dimension (:), allocatable :: c22, c23, c24, c25
+ double precision, dimension (:), allocatable :: c34, c35, c36, c44
+ double precision, dimension (:), allocatable :: c55, c56, c66, c45
+ double precision, dimension (:), allocatable :: c46
+
+ end type par
+
+ type (par) :: reg1Bc, reg2Bc, reg3Bc
+
+end module cem_par
+
+subroutine model_cem_broadcast ( myrank )
+
+ use :: cem_par
+ use :: netcdf
+ use :: meshfem3D_models_par, only: CEM_ACCEPT
+
+ integer, intent (in) :: myrank
+ integer :: wSize
+
+ rank = myrank
+ call world_size (wSize)
+
+ if ( CEM_ACCEPT ) then
+
+ call return_populated_arrays (reg1Bc, "vsv", 1)
+ call return_populated_arrays (reg2Bc, "vsv", 2)
+ call return_populated_arrays (reg3Bc, "vsv", 3)
+
+ call return_populated_arrays (reg1Bc, "rho", 1)
+ call return_populated_arrays (reg2Bc, "rho", 2)
+ call return_populated_arrays (reg3Bc, "rho", 3)
+
+ call return_populated_arrays (reg1Bc, "vsh", 1)
+ call return_populated_arrays (reg2Bc, "vsh", 2)
+ call return_populated_arrays (reg3Bc, "vsh", 3)
+
+ call return_populated_arrays (reg1Bc, "vpp", 1)
+ call return_populated_arrays (reg2Bc, "vpp", 2)
+ call return_populated_arrays (reg3Bc, "vpp", 3)
+
+ call synchronize_all ()
+ if (myrank == 0 ) print *, "Finished reading in netcdf."
+
+ end if
+
+end subroutine model_cem_broadcast
+
+subroutine request_cem ( vsh, vsv, vph, vpv, rho, iregion_code, ispec, i, j, k )
+
+ use :: cem_par
+ use :: constants
+
+ use :: meshfem3D_par, only: ibool
+
+ implicit none
+
+ double precision, intent (out) :: vsh, vsv, vph, vpv, rho
+
+ integer, intent (in) :: iregion_code, ispec, i, j, k
+ integer :: iglob
+
+ if ( iregion_code == IREGION_CRUST_MANTLE ) then
+
+ iglob = ibool(i,j,k,ispec)
+
+ rho = reg1Bc%rho(iglob) * 1000.0d0 / (RHOAV)
+ vpv = reg1Bc%vpv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+ vsv = reg1Bc%vsv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+ vsh = reg1Bc%vsh(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+
+ else if ( iregion_code == IREGION_OUTER_CORE ) then
+
+ iglob = ibool(i,j,k,ispec)
+
+ rho = reg2Bc%rho(iglob) * 1000.0d0 / (RHOAV)
+ vpv = reg2Bc%vpv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+ vsv = reg2Bc%vsv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+ vsh = reg2Bc%vsh(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+
+ else if ( iregion_code == IREGION_INNER_CORE ) then
+
+ iglob = ibool(i,j,k,ispec)
+
+ rho = reg3Bc%rho(iglob) * 1000.0d0 / (RHOAV)
+ vpv = reg3Bc%vpv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+ vsv = reg3Bc%vsv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+ vsh = reg3Bc%vsh(iglob) * 1000.0d0 / (R_EARTH * scaleval)
+
+ end if
+
+ vph = vpv
+
+end subroutine request_cem
+
+subroutine return_populated_arrays ( structure, param, reg )
+
+ use :: cem_par
+ use :: netcdf
+
+ implicit none
+
+ integer, intent (in) :: reg
+ integer :: ncid, dimid, nPar
+ integer :: status, varid
+
+
+ character(3), intent(in) :: param
+
+ type (par) :: structure
+
+ character (len=1024) :: fileName
+ character (len=1024) :: fileNameTrim
+ character (len=1024) :: formatString
+
+
+ formatString = "(A,I0.2,A,I0.4,A1,A3,A3)"
+ write (fileName,formatString) "./DATA/cemRequest/xyz_reg", reg, &
+ "_proc", rank, ".", param, ".nc"
+
+ fileNameTrim = trim(fileName)
+
+ status = nf90_open (fileNameTrim, NF90_NOWRITE, ncid)
+ status = nf90_inq_dimid (ncid, "param", dimid)
+ status = nf90_inquire_dimension (ncid, dimid, len = nPar)
+
+ if ( param == "vsv" ) allocate(structure%vsv(nPar))
+ if ( param == "rho" ) allocate(structure%rho(nPar))
+ if ( param == "vpp" ) allocate(structure%vpv(nPar))
+ if ( param == "vsh" ) allocate(structure%vsh(nPar))
+
+ status = nf90_inq_varid (ncid, "data", varid)
+
+ if ( param == "vsv" ) status = nf90_get_var (ncid, varid, structure%vsv)
+ if ( param == "rho" ) status = nf90_get_var (ncid, varid, structure%rho)
+ if ( param == "vpp" ) status = nf90_get_var (ncid, varid, structure%vpv)
+ if ( param == "vsh" ) status = nf90_get_var (ncid, varid, structure%vsh)
+
+
+end subroutine return_populated_arrays
+
+subroutine write_cem_request ( iregion_code )
+
+ use :: netcdf
+ use :: constants
+ use :: CEM_par
+
+ implicit none
+
+ integer :: ncid, status, x_dimind, y_dimind, z_dimind, varidX, varidY, varidZ
+ integer :: varidR, r_dimind, iregion_code
+ character (len = 1024) :: fileName, fileNameTrim, formatString
+
+! This line does not seem portable.
+! call execute_command_line ('mkdir -p cemRequest/')
+
+ formatString = "(A,I0.2,A,I0.4)"
+ write (fileName,formatString) "DATA/cemRequest/xyz_reg", iregion_code, &
+ "_proc", rank
+ fileNameTrim = trim(fileName)
+
+ status = nf90_create (path = fileNameTrim, cmode = NF90_CLOBBER, ncid = ncid)
+
+ status = nf90_def_dim (ncid, "x", size (xyzOut(:,1)), x_dimind )
+ status = nf90_def_dim (ncid, "y", size (xyzOut(:,2)), y_dimind )
+ status = nf90_def_dim (ncid, "z", size (xyzOut(:,3)), z_dimind )
+ status = nf90_def_dim (ncid, "r", size (regCode), r_dimind )
+
+ status = nf90_def_var (ncid, "dataX", NF90_FLOAT, x_dimind, varidX)
+ status = nf90_def_var (ncid, "dataY", NF90_FLOAT, y_dimind, varidY)
+ status = nf90_def_var (ncid, "dataZ", NF90_FLOAT, z_dimind, varidZ)
+ status = nf90_def_var (ncid, "regC_", NF90_SHORT, r_dimind, varidR)
+
+ status = nf90_def_var_deflate (ncid, varidX, shuOn, comOn, comLvl)
+ status = nf90_def_var_deflate (ncid, varidY, shuOn, comOn, comLvl)
+ status = nf90_def_var_deflate (ncid, varidZ, shuOn, comOn, comLvl)
+ status = nf90_def_var_deflate (ncid, varidR, shuOn, comOn, comLvl)
+
+ status = nf90_enddef (ncid)
+
+ status = nf90_put_var (ncid, varidX, xyzOut(:,1))
+ status = nf90_put_var (ncid, varidY, xyzOut(:,2))
+ status = nf90_put_var (ncid, varidZ, xyzOut(:,3))
+ status = nf90_put_var (ncid, varidR, regCode)
+ status = nf90_close (ncid)
+
+ deallocate(xyzOut)
+ deallocate(regCode)
+
+end subroutine write_cem_request
+
+subroutine build_global_coordinates ( nspec, nglob, iregion_code )
+
+ use :: constants
+ use :: cem_par
+
+ use :: meshfem3D_par, only: &
+ ibool,idoubling,xstore,ystore,zstore
+
+ implicit none
+
+ integer, intent (in) :: nspec, nglob, iregion_code
+ integer :: i, j, k, iglob, ispec, region
+
+ double precision, parameter :: R_020_KM=6351.0d0, R_052_KM=6319.0d0
+ double precision, parameter :: R_100_KM=6271.0d0, R_400_KM=5971.0d0
+ double precision, parameter :: R_670_KM=5701.0d0, R_OCR_KM=3480.0d0
+ double precision, parameter :: R_ICR_KM=1221.0d0, R_THO_KM=5371.0d0
+ double precision :: x, y, z, rad
+
+ allocate(xyzOut(nglob,NDIMS))
+ allocate(regCode(nglob))
+
+ ! Global to local x co-ordinate numbering.
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ x = sngl(xstore(i,j,k,ispec)) * R_EARTH_KM
+ iglob = ibool(i,j,k,ispec)
+ xyzOut(iglob,1) = x
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! Global to local y co-ordinate numbering.
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ y = sngl(ystore(i,j,k,ispec)) * R_EARTH_KM
+ iglob = ibool(i,j,k,ispec)
+ xyzOut(iglob,2) = y
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ ! Global to local z co-ordinate numbering.
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ z = sngl(zstore(i,j,k,ispec)) * R_EARTH_KM
+ iglob = ibool(i,j,k,ispec)
+ xyzOut(iglob,3) = z
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ if (iregion_code == 1) then
+
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ x = sngl(xstore(i,j,k,ispec)) * R_EARTH_KM
+ iglob = ibool(i,j,k,ispec)
+ xyzOut(iglob,1) = x
+
+ y = sngl(ystore(i,j,k,ispec)) * R_EARTH_KM
+ iglob = ibool(i,j,k,ispec)
+ xyzOut(iglob,2) = y
+
+ z = sngl(zstore(i,j,k,ispec)) * R_EARTH_KM
+ iglob = ibool(i,j,k,ispec)
+ xyzOut(iglob,3) = z
+
+ rad = dsqrt( x * x + y * y + z * z )
+
+ if ( rad < R_THO_KM .and. rad >= R_ICR_KM ) then
+ region = 7
+ else if ( rad < R_670_KM .and. rad >= R_THO_KM ) then
+ region = 6
+ else if ( rad < R_400_KM .and. rad >= R_670_KM ) then
+ region = 5
+ else if ( rad < R_100_KM .and. rad >= R_400_KM ) then
+ region = 4
+ else if ( rad < R_052_KM .and. rad >= R_100_KM ) then
+ region = 3
+ else if ( rad < R_020_KM .and. rad >= R_052_KM ) then
+ region = 2
+ else if ( rad >= R_020_KM ) then
+ region = 1
+ end if
+
+ iglob = ibool(i,j,k,ispec)
+ regCode(iglob) = region
+
+ end do
+ end do
+ end do
+ end do
+
+ else if (iregion_code == 2) then
+
+ regCode(:) = 8
+
+ else if (iregion_code == 3) then
+
+ regCode(:) = 9
+
+ end if
+
+end subroutine build_global_coordinates
diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk
index 9b18936..c0e4a59 100644
--- a/src/meshfem3D/rules.mk
+++ b/src/meshfem3D/rules.mk
@@ -227,6 +227,14 @@ meshfem3D_OBJECTS += $(adios_meshfem3D_STUBS)
meshfem3D_SHARED_OBJECTS += $(adios_meshfem3D_SHARED_STUBS)
endif
+# conditional CEM model
+ifeq ($(CEM),yes)
+meshfem3D_OBJECTS += $O/model_cem.checknetcdf.o
+NETCDF_INCLUDE = -I$(netcdf_specfem_dir)/include -I$(hdf5_specfem_dir)/include -I$(zlib_specfem_dir)/include
+LDFLAGS += -L$(netcdf_specfem_dir)/lib -lnetcdff -lnetcdf -L$(hdf5_specfem_dir)/lib -lhdf5_hl -lhdf5 -L$(zlib_specfem_dir)/lib -lz
+endif
+
+
#######################################
####
@@ -292,3 +300,7 @@ $O/%.check_adios.o: $S/%.f90 $O/shared_par.shared_module.o $O/meshfem3D_par.chec
$O/%.check_adios.o: $S/%.F90 $O/shared_par.shared_module.o $O/meshfem3D_par.check_module.o $O/adios_helpers.shared_adios.o
${FCCOMPILE_CHECK} ${FCFLAGS_f90} -c -o $@ $<
+## CEM
+
+$O/%.checknetcdf.o: $S/%.f90 $O/shared_par.shared_module.o $O/meshfem3D_par.check_module.o
+ ${FCCOMPILE_CHECK} ${FCFLAGS_f90} $(NETCDF_INCLUDE) -c -o $@ $<
More information about the CIG-COMMITS
mailing list