[cig-commits] r12456 - seismo/2D/SPECFEM2D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Jul 21 07:57:24 PDT 2008
Author: dkomati1
Date: 2008-07-21 07:57:24 -0700 (Mon, 21 Jul 2008)
New Revision: 12456
Removed:
seismo/2D/SPECFEM2D/trunk/constants_unstruct.h
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/attenuation_compute_param.c
seismo/2D/SPECFEM2D/trunk/checkgrid.F90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/get_perm_cuthill_mckee.f90
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt
Log:
cleaned the source code, following the commit of version 3 developed in Barcelona;
done cleaning everything except specfem2D.F90
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2008-07-21 14:57:24 UTC (rev 12456)
@@ -71,7 +71,7 @@
#F90 = /opt/openmpi-1.2.1/gfortran64/bin/mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
CC = gcc
#FLAGS_NOCHECK = -O3 -march=opteron -m64 -mfpmath=sse,387
-FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O2 -Wunused-labels -Waliasing -Wampersand -Wsurprising -Wline-truncation -Wunderflow -fbounds-check
+FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -mcmodel=medium
FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
# IBM
@@ -231,7 +231,7 @@
$O/write_seismograms.o: write_seismograms.F90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/write_seismograms.o write_seismograms.F90
-$O/part_unstruct.o: part_unstruct.F90 constants_unstruct.h
+$O/part_unstruct.o: part_unstruct.F90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/part_unstruct.o part_unstruct.F90
$O/construct_acoustic_surface.o: construct_acoustic_surface.f90 constants.h
Modified: seismo/2D/SPECFEM2D/trunk/attenuation_compute_param.c
===================================================================
--- seismo/2D/SPECFEM2D/trunk/attenuation_compute_param.c 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/attenuation_compute_param.c 2008-07-21 14:57:24 UTC (rev 12456)
@@ -17,12 +17,12 @@
#define PI2 6.28318530717958
/* Underscores should or should not follow this function name, depending on the compiler and its options.
- It is called in "attenuation_model.f90".
+ It is called in "attenuation_model.f90".
*/
-int attenuation_compute_param_(int *nmech_in, double *Qp_in, double *Qs_in, double *f1_in, double *f2_in,
- double *tau_sigma_nu1, double *tau_sigma_nu2,
- double *tau_epsilon_nu1, double *tau_epsilon_nu2
- )
+int attenuation_compute_param_(int *nmech_in, double *Qp_in, double *Qs_in, double *f1_in, double *f2_in,
+ double *tau_sigma_nu1, double *tau_sigma_nu2,
+ double *tau_epsilon_nu1, double *tau_epsilon_nu2
+ )
{
int xmgr, n, i, j, plot, nu;
@@ -120,35 +120,35 @@
/* output in Fortran90 format */
for (i = 1; i <= n; i++) {
- /*
+ /*
printf(" tau_sigma_nu%d(%1d) = %30.20lfd0\n", nu, i, tau_s[i]);
*/
/* We put the results in tau_sigma_nu to get them in fortran. */
if ( nu == 1 ) {
- tau_sigma_nu1[i-1] = tau_s[i];
+ tau_sigma_nu1[i-1] = tau_s[i];
}
if ( nu == 2 ) {
- tau_sigma_nu2[i-1] = tau_s[i];
+ tau_sigma_nu2[i-1] = tau_s[i];
}
-
- }
+
+ }
//printf("\n");
-
+
for (i = 1; i <= n; i++) {
/*
- printf(" tau_epsilon_nu%d(%1d) = %30.20lfd0\n", nu, i, tau_e[i]);
+ printf(" tau_epsilon_nu%d(%1d) = %30.20lfd0\n", nu, i, tau_e[i]);
*/
/* We put the results in tau_epsilon_nu to get them in fortran. */
if ( nu == 1 ) {
- tau_epsilon_nu1[i-1] = tau_e[i];
+ tau_epsilon_nu1[i-1] = tau_e[i];
}
if ( nu == 2 ) {
- tau_epsilon_nu2[i-1] = tau_e[i];
+ tau_epsilon_nu2[i-1] = tau_e[i];
}
-
+
}
//printf("\n");
-
+
free_dvector(tau_s, 1, n);
free_dvector(tau_e, 1, n);
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2008-07-21 14:57:24 UTC (rev 12456)
@@ -41,9 +41,8 @@
!========================================================================
subroutine checkgrid(vpext,vsext,rhoext,density,elastcoef,ibool,kmato,coord,npoin,vpmin,vpmax, &
- assign_external_model,nspec,numat,deltat,f0,t0,initialfield,time_function_type, &
- coorg,xinterp,zinterp,shapeint,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic,myrank,nproc, &
- nspec_outer,nspec_inner)
+ assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat,f0,t0,initialfield,time_function_type, &
+ coorg,xinterp,zinterp,shapeint,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic,myrank,nproc)
! check the mesh, stability and number of points per wavelength
@@ -54,8 +53,9 @@
include 'mpif.h'
#endif
-!! DK DK to analyze Cuthill-McKee display
- integer :: UPPER_LIMIT_DISPLAY,nspec_outer,nspec_inner
+! option to display only part of the mesh and not the whole mesh,
+! for instance to analyze Cuthill-McKee mesh partitioning etc.
+ integer :: UPPER_LIMIT_DISPLAY
! color palette
integer, parameter :: NUM_COLORS = 236
@@ -127,12 +127,6 @@
! title of the plot
character(len=60) simulation_title
-!! DK DK to analyze Cuthill-McKee display
-! UPPER_LIMIT_DISPLAY = nspec
- UPPER_LIMIT_DISPLAY = nspec_outer
-! UPPER_LIMIT_DISPLAY = nspec_inner
-! UPPER_LIMIT_DISPLAY = 2300
-
if(UPPER_LIMIT_DISPLAY > nspec) stop 'cannot have UPPER_LIMIT_DISPLAY > nspec in checkgrid.F90'
#ifndef USE_MPI
@@ -147,7 +141,6 @@
deallocate(greyscale_recv)
#endif
-
! define percentage of smallest distance between GLL points for NGLLX points
! percentages were computed by calling the GLL points routine for each degree
percent_GLL(2) = 100.d0
@@ -2109,8 +2102,6 @@
coorg_send(2,(ispec-1)*5+5) = z2
end if
-
-
material = kmato(ispec)
mu = elastcoef(2,material)
@@ -2283,7 +2274,6 @@
end if
#endif
-
if ( myrank == 0 ) then
write(24,*) '%'
write(24,*) 'grestore'
@@ -2570,10 +2560,9 @@
end if
-
if ( myrank == 0 ) then
print *
- print *,'Creating PostScript file with partitioning'
+ print *,'Creating PostScript file with mesh partitioning'
!
!---- open PostScript file
!
@@ -2767,7 +2756,6 @@
coorg_send(2,(ispec-1)*5+5) = z2
end if
-
if ( myrank == 0 ) then
write(24,*) red(1), green(1), blue(1), 'RG GF 0 setgray ST'
end if
@@ -2821,8 +2809,6 @@
print *,'End of creation of PostScript file with partitioning'
end if
-
-
10 format('%!PS-Adobe-2.0',/,'%%',/,'%% Title: ',a50,/,'%% Created by: Specfem2D',/,'%% Author: Dimitri Komatitsch',/,'%%')
681 format(f6.2,1x,f6.2)
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2008-07-21 14:57:24 UTC (rev 12456)
@@ -84,7 +84,6 @@
! for overlapping MPI communications with computation
integer, intent(in) :: nspec_outer
-!!!!!!! integer, dimension(max(1,nspec_outer)), intent(in) :: ispec_inner_outer_to_glob
logical, intent(in) :: we_are_in_phase_outer
!---
@@ -105,11 +104,6 @@
! material properties of the elastic medium
real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,kappal,cpl,rhol
-! loop over spectral elements
-! do ispec_inner_outer = 1,nspec_outer
-
-! ispec = ispec_inner_outer_to_glob(ispec_inner_outer)
-
integer :: ifirstelem,ilastelem
if(we_are_in_phase_outer) then
@@ -120,6 +114,7 @@
ilastelem = nspec
endif
+! loop over spectral elements
do ispec = ifirstelem,ilastelem
!---
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-07-21 14:57:24 UTC (rev 12456)
@@ -101,7 +101,6 @@
! for overlapping MPI communications with computation
integer, intent(in) :: nspec_outer
-!!!!!!!!!!!!!!!!! integer, dimension(max(1,nspec_outer)), intent(in) :: ispec_inner_outer_to_glob
logical, intent(in) :: we_are_in_phase_outer
!---
@@ -150,12 +149,6 @@
dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
endif
-! loop over spectral elements
-! do ispec_inner_outer = 1,nspec_outer
-
-! get global numbering for inner or outer elements
-! ispec = ispec_inner_outer_to_glob(ispec_inner_outer)
-
if(we_are_in_phase_outer) then
ifirstelem = 1
ilastelem = nspec_outer
@@ -164,6 +157,7 @@
ilastelem = nspec
endif
+! loop over spectral elements
do ispec = ifirstelem,ilastelem
!---
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2008-07-21 14:57:24 UTC (rev 12456)
@@ -16,7 +16,11 @@
! the code does NOT work if NGLLZ /= NGLLX because it then cannot handle a non-structured mesh
integer, parameter :: NGLLZ = NGLLX
-! for Cuthill-McKee (1969) permutation
+! further reduce cache misses inner/outer in two passes in the case of an MPI simulation
+! this flag is ignored in the case of a serial simulation
+ logical, parameter :: FURTHER_REDUCE_CACHE_MISSES = .true.
+
+! for inverse Cuthill-McKee (1969) permutation
logical, parameter :: PERFORM_CUTHILL_MCKEE = .true.
logical, parameter :: INVERSE = .true.
logical, parameter :: FACE = .false.
@@ -26,6 +30,37 @@
! maximum size if multi-level Cuthill-McKee ordering
integer, parameter :: LIMIT_MULTI_CUTHILL = 50
+! implement Cuthill-McKee or replace with identity permutation
+ logical, parameter :: ACTUALLY_IMPLEMENT_PERM_OUT = .true.
+ logical, parameter :: ACTUALLY_IMPLEMENT_PERM_INN = .true.
+
+! add MPI barriers and suppress seismograms if we generate traces of the run for analysis with "ParaVer"
+ logical, parameter :: GENERATE_PARAVER_TRACES = .false.
+
+! option to display only part of the mesh and not the whole mesh,
+! for instance to analyze Cuthill-McKee mesh partitioning etc.
+! Possible values are:
+! 1: display all the elements (i.e., the whole mesh)
+! 2: display inner elements only
+! 3: display outer elements only
+! 4: display a fixed number of elements (in each partition) only
+ integer, parameter :: DISPLAY_SUBSET_OPTION = 1
+! number of spectral elements to display in each subset when a fixed subset size is used (option 4 above)
+ integer, parameter :: NSPEC_DISPLAY_SUBSET = 2300
+
+!--- beginning of Nicolas Le Goff's constants for an unstructured CUBIT/METIS/SCOTCH mesh
+
+! number of nodes per element
+ integer, parameter :: ESIZE = 4
+
+! maximum number of neighbors per element
+ integer, parameter :: max_neighbor = 30
+
+! maximum number of elements that can contain the same node
+ integer, parameter :: nsize = 20
+
+!--- end of Nicolas Le Goff's constants for an unstructured CUBIT/METIS/SCOTCH mesh
+
! compute and output acoustic and elastic energy (slows down the code significantly)
logical, parameter :: OUTPUT_ENERGY = .false.
@@ -35,10 +70,6 @@
! select fast (Paul Fischer) or slow (topology only) global numbering algorithm
logical, parameter :: FAST_NUMBERING = .true.
-! further reduce cache misses inner/outer in two passes in the case of an MPI simulation
-! this flag is ignored in the case of a serial simulation
- logical, parameter :: FURTHER_REDUCE_CACHE_MISSES = .true.
-
! mesh tolerance for fast global numbering
double precision, parameter :: SMALLVALTOL = 0.00001d0
Deleted: seismo/2D/SPECFEM2D/trunk/constants_unstruct.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants_unstruct.h 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/constants_unstruct.h 2008-07-21 14:57:24 UTC (rev 12456)
@@ -1,9 +0,0 @@
-
-! Number of nodes per elements.
-integer, parameter :: ESIZE = 4
-
-! Max number of neighbours per elements.
- integer :: max_neighbour=30
-
-! Max number of elements that can contain the same node.
- integer :: nsize=20
Modified: seismo/2D/SPECFEM2D/trunk/get_perm_cuthill_mckee.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/get_perm_cuthill_mckee.f90 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/get_perm_cuthill_mckee.f90 2008-07-21 14:57:24 UTC (rev 12456)
@@ -1,29 +1,44 @@
-!=====================================================================
+
+!========================================================================
!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
-! --------------------------------------------------
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory, California Institute of Technology, USA
-! and University of Pau / CNRS / INRIA, France
-! (c) California Institute of Technology and University of Pau / CNRS / INRIA
-! February 2008
+! Copyright Universite de Pau et des Pays de l'Adour, CNRS and INRIA, France.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
!
-! 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 software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic wave equation
+! using a spectral-element method (SEM).
!
-! 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.
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
!
-! 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.
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
!
-!=====================================================================
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
! implement reverse Cuthill-McKee (1969) ordering, introduced in
! E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices.
@@ -361,6 +376,12 @@
nelem = ne(jel)
if(maskel(nelem)) goto 120
if (FACE) then
+!! DK DK this below implemented by David Michea in 3D, but not true anymore in 2D: should be
+!! DK DK two corners instead of three. But does not matter because FACE is always .false.
+!! DK DK and therefore this part of the routine is currently never used.
+!! DK DK Let me add a stop statement just in case.
+ stop 'FACE = .true. not implemented, check the above comment in the source code'
+!! DK DK End of the stop statement added.
! if 2 elements share at least 3 corners, therefore they share a face
countel(nelem) = countel(nelem) + 1
if (countel(nelem)>=3) then
Modified: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2008-07-21 14:57:24 UTC (rev 12456)
@@ -49,8 +49,6 @@
implicit none
- include './constants_unstruct.h'
-
contains
!-----------------------------------------------
@@ -60,6 +58,8 @@
!-----------------------------------------------
subroutine read_mesh(filename, nelmnts, elmnts, nnodes, num_start)
+ include "constants.h"
+
character(len=256), intent(in) :: filename
integer, intent(out) :: nelmnts
integer, intent(out) :: nnodes
@@ -68,8 +68,6 @@
integer :: i
- print *, trim(filename)
-
open(unit=990, file=trim(filename), form='formatted' , status='old', action='read')
read(990,*) nelmnts
allocate(elmnts(0:ESIZE*nelmnts-1))
@@ -83,10 +81,8 @@
elmnts(:) = elmnts(:) - num_start
nnodes = maxval(elmnts) + 1
-
end subroutine read_mesh
-
!-----------------------------------------------
! Read the nodes coordinates and storing it in array 'nodes_coords'
!-----------------------------------------------
@@ -98,14 +94,11 @@
integer :: i
- print *, trim(filename)
-
open(unit=991, file=trim(filename), form='formatted' , status='old', action='read')
read(991,*) nnodes
allocate(nodes_coords(2,nnodes))
do i = 1, nnodes
read(991,*) nodes_coords(1,i), nodes_coords(2,i)
-
enddo
close(991)
@@ -123,12 +116,9 @@
integer :: i
- print *, trim(filename)
-
open(unit=992, file=trim(filename), form='formatted' , status='old', action='read')
do i = 1, nelmnts
read(992,*) num_material(i)
-
enddo
close(992)
@@ -144,7 +134,7 @@
subroutine read_acoustic_surface(filename, nelem_acoustic_surface, acoustic_surface, &
nelmnts, num_material, ANISOTROPIC_MATERIAL, nb_materials, icodemat, cs, num_start)
- include './constants.h'
+ include "constants.h"
character(len=256), intent(in) :: filename
integer, intent(out) :: nelem_acoustic_surface
@@ -210,7 +200,7 @@
!-----------------------------------------------
subroutine read_abs_surface(filename, nelemabs, abs_surface, num_start)
- include './constants.h'
+ include "constants.h"
character(len=256), intent(in) :: filename
integer, intent(out) :: nelemabs
@@ -246,6 +236,8 @@
!-----------------------------------------------
subroutine mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts, ncommonnodes)
+ include "constants.h"
+
integer, intent(in) :: nelmnts
integer, intent(in) :: nnodes
integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
@@ -261,10 +253,9 @@
integer :: elem_base, elem_target
integer :: connectivity
-
allocate(xadj(0:nelmnts))
xadj(:) = 0
- allocate(adjncy(0:max_neighbour*nelmnts-1))
+ allocate(adjncy(0:max_neighbor*nelmnts-1))
adjncy(:) = 0
allocate(nnodes_elmnts(0:nnodes-1))
nnodes_elmnts(:) = 0
@@ -280,8 +271,6 @@
enddo
- print *, 'nnodes_elmnts'
-
! checking which elements are neighbours ('ncommonnodes' criteria)
do j = 0, nnodes-1
do k = 0, nnodes_elmnts(j)-1
@@ -305,16 +294,16 @@
do m = 0, xadj(nodes_elmnts(k+j*nsize))
if ( .not.is_neighbour ) then
- if ( adjncy(nodes_elmnts(k+j*nsize)*max_neighbour+m) == nodes_elmnts(l+j*nsize) ) then
+ if ( adjncy(nodes_elmnts(k+j*nsize)*max_neighbor+m) == nodes_elmnts(l+j*nsize) ) then
is_neighbour = .true.
endif
endif
enddo
if ( .not.is_neighbour ) then
- adjncy(nodes_elmnts(k+j*nsize)*max_neighbour+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
+ adjncy(nodes_elmnts(k+j*nsize)*max_neighbor+xadj(nodes_elmnts(k+j*nsize))) = nodes_elmnts(l+j*nsize)
xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
- adjncy(nodes_elmnts(l+j*nsize)*max_neighbour+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
+ adjncy(nodes_elmnts(l+j*nsize)*max_neighbor+xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
endif
endif
@@ -327,7 +316,7 @@
k = xadj(i)
xadj(i) = nb_edges
do j = 0, k-1
- adjncy(nb_edges) = adjncy(i*max_neighbour+j)
+ adjncy(nb_edges) = adjncy(i*max_neighbor+j)
nb_edges = nb_edges + 1
enddo
enddo
@@ -392,6 +381,8 @@
subroutine Construct_glob2loc_nodes(nelmnts, nnodes, nnodes_elmnts, nodes_elmnts, part, nparts, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
+ include "constants.h"
+
integer, intent(in) :: nelmnts, nnodes, nparts
integer, dimension(0:nelmnts-1), intent(in) :: part
integer, dimension(0:nnodes-1), intent(in) :: nnodes_elmnts
@@ -475,13 +466,13 @@
subroutine Construct_interfaces(nelmnts, nparts, part, elmnts, xadj, adjncy, tab_interfaces, &
tab_size_interfaces, ninterfaces, nb_materials, cs_material, num_material)
- include 'constants.h'
+ include "constants.h"
integer, intent(in) :: nelmnts, nparts
integer, dimension(0:nelmnts-1), intent(in) :: part
integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
integer, dimension(0:nelmnts), intent(in) :: xadj
- integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(0:max_neighbor*nelmnts-1), intent(in) :: adjncy
integer, dimension(:),pointer :: tab_size_interfaces, tab_interfaces
integer, intent(out) :: ninterfaces
integer, dimension(1:nelmnts), intent(in) :: num_material
@@ -574,6 +565,7 @@
tab_interfaces(tab_size_interfaces(num_interface)*5+num_edge*5+2) = ncommon_nodes
else
print *, "Error while building interfaces!", ncommon_nodes
+ stop 'fatal error'
endif
num_edge = num_edge + 1
endif
@@ -586,7 +578,6 @@
enddo
enddo
-
end subroutine Construct_interfaces
@@ -880,7 +871,7 @@
elmnts, ngnod)
implicit none
- include 'constants.h'
+ include "constants.h"
integer, intent(inout) :: nelemabs
integer, intent(out) :: nelemabs_merge
@@ -1167,10 +1158,12 @@
!--------------------------------------------------
subroutine Part_metis(nelmnts, xadj, adjncy, vwgt, adjwgt, nparts, nb_edges, edgecut, part, metis_options)
+ include "constants.h"
+
integer, intent(in) :: nelmnts, nparts, nb_edges
integer, intent(inout) :: edgecut
integer, dimension(0:nelmnts), intent(in) :: xadj
- integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(0:max_neighbor*nelmnts-1), intent(in) :: adjncy
integer, dimension(0:nelmnts-1), intent(in) :: vwgt
integer, dimension(0:nb_edges-1), intent(in) :: adjwgt
integer, dimension(:), pointer :: part
@@ -1182,14 +1175,11 @@
num_start = 0
wgtflag = 0
- print *, 'avant', edgecut
call METIS_PartGraphRecursive(nelmnts, xadj(0), adjncy(0), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
metis_options, edgecut, part(0));
!call METIS_PartGraphVKway(nelmnts, xadj(0), adjncy(0), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
! options, edgecut, part(0));
- print *, 'apres', edgecut
-
end subroutine Part_metis
#endif
@@ -1200,12 +1190,14 @@
!--------------------------------------------------
subroutine Part_scotch(nelmnts, xadj, adjncy, vwgt, adjwgt, nparts, nedges, edgecut, part, scotch_strategy)
- include 'scotchf.h'
+ include "constants.h"
+ include "scotchf.h"
+
integer, intent(in) :: nelmnts, nparts, nedges
integer, intent(inout) :: edgecut
integer, dimension(0:nelmnts), intent(in) :: xadj
- integer, dimension(0:max_neighbour*nelmnts-1), intent(in) :: adjncy
+ integer, dimension(0:max_neighbor*nelmnts-1), intent(in) :: adjncy
integer, dimension(0:nelmnts-1), intent(in) :: vwgt
integer, dimension(:), pointer :: adjwgt
integer, dimension(:), pointer :: part
@@ -1284,7 +1276,7 @@
implicit none
- include 'constants.h'
+ include "constants.h"
integer, intent(in) :: nelmnts, nnodes, nproc, nb_materials
double precision, dimension(nb_materials), intent(in) :: cs_material
@@ -1327,8 +1319,6 @@
endif
enddo
- print *, 'nedges_coupled', nedges_coupled
-
allocate(edges_coupled(2,nedges_coupled))
nedges_coupled = 0
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-07-21 14:57:24 UTC (rev 12456)
@@ -144,10 +144,6 @@
implicit none
-! implement Cuthill-McKee or replace with identity permutation
- logical, parameter :: ACTUALLY_IMPLEMENT_PERM_OUT = .true.
- logical, parameter :: ACTUALLY_IMPLEMENT_PERM_INN = .true.
-
include "constants.h"
#ifdef USE_MPI
include "mpif.h"
@@ -236,7 +232,7 @@
double precision :: vpmin,vpmax
integer :: colors,numbers,subsamp,imagetype,NTSTEP_BETWEEN_OUTPUT_INFO,NTSTEP_BETWEEN_OUTPUT_SEISMO,seismotype
- integer :: numat,ngnod,nspec,pointsdisp,nelemabs,nelem_acoustic_surface,ispecabs
+ integer :: numat,ngnod,nspec,pointsdisp,nelemabs,nelem_acoustic_surface,ispecabs,UPPER_LIMIT_DISPLAY
logical interpol,meshvect,modelvect,boundvect,assign_external_model,initialfield, &
outputgrid,gnuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,output_postscript_snapshot,output_color_image, &
@@ -1595,9 +1591,6 @@
endif
-!! DK DK for ParaVer tests
- call MPI_BARRIER(MPI_COMM_WORLD,ier)
-
enddo ! end of further reduction of cache misses inner/outer in two passes
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
@@ -1608,12 +1601,21 @@
if(any_elastic) rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
if(any_acoustic) rmass_inverse_acoustic(:) = 1._CUSTOM_REAL / rmass_inverse_acoustic(:)
-
! check the mesh, stability and number of points per wavelength
+ if(DISPLAY_SUBSET_OPTION == 1) then
+ UPPER_LIMIT_DISPLAY = nspec
+ else if(DISPLAY_SUBSET_OPTION == 2) then
+ UPPER_LIMIT_DISPLAY = nspec_inner
+ else if(DISPLAY_SUBSET_OPTION == 3) then
+ UPPER_LIMIT_DISPLAY = nspec_outer
+ else if(DISPLAY_SUBSET_OPTION == 4) then
+ UPPER_LIMIT_DISPLAY = NSPEC_DISPLAY_SUBSET
+ else
+ stop 'incorrect value of DISPLAY_SUBSET_OPTION'
+ endif
call checkgrid(vpext,vsext,rhoext,density,elastcoef,ibool,kmato,coord,npoin,vpmin,vpmax, &
- assign_external_model,nspec,numat,deltat,f0,t0,initialfield,time_function_type, &
- coorg,xinterp,zinterp,shape2D_display,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic,myrank,nproc, &
- nspec_outer,nspec_inner)
+ assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat,f0,t0,initialfield,time_function_type, &
+ coorg,xinterp,zinterp,shape2D_display,knods,simulation_title,npgeo,pointsdisp,ngnod,any_elastic,myrank,nproc)
! convert receiver angle to radians
anglerec = anglerec * pi / 180.d0
@@ -1751,7 +1753,6 @@
enddo
enddo
-
! filling array iglob_image_color, containing info on which process owns which pixels.
#ifdef USE_MPI
allocate(nb_pixel_per_proc(nproc))
@@ -2441,8 +2442,10 @@
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
-!! DK DK for ParaVer tests
- call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#ifdef USE_MPI
+! add a barrier if we generate traces of the run for analysis with "ParaVer"
+ if(GENERATE_PARAVER_TRACES) call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#endif
do it = 1,NSTEP
@@ -3104,7 +3107,8 @@
endif
!---- save temporary or final seismograms
- call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
+! suppress seismograms if we generate traces of the run for analysis with "ParaVer", because time consuming
+ if(.not. GENERATE_PARAVER_TRACES) call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current)
@@ -3137,8 +3141,10 @@
endif
-!! DK DK for ParaVer tests
- call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#ifdef USE_MPI
+! add a barrier if we generate traces of the run for analysis with "ParaVer"
+ if(GENERATE_PARAVER_TRACES) call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#endif
enddo ! end of the main time loop
Modified: seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt 2008-07-21 14:07:04 UTC (rev 12455)
+++ seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt 2008-07-21 14:57:24 UTC (rev 12456)
@@ -1,6 +1,4 @@
-- splitting file part_unstruct.F90 in several files for clarity purpose.
-
- improving compiling with SCOTCH (issue with header file scotchf.h which is Fortran77 legal). Having our own scotchf.h file (without the comments) is not wise.
- comparing the different partitioning methods for METIS and SCOTCH, and finding a good default for SCOTCH.
@@ -9,22 +7,12 @@
- choosing a way to use assign_external_model?
-- adding comments.
-
- checking for points with different normals for absorbing conditions, when the absorbing edges are not in the same elements (similar to what is done for the corners).
-- scripts for translating GID/CUBIT meshes into files for xmeshfem2D.
+- scripts to translate GID/CUBIT meshes into files for xmeshfem2D.
-- modifying scripts for UPPA cluster (when FS sync issues are solved and remote commands are available).
+- user manual for unstructured meshes.
-- manual for unstructured meshes.
-
-- getting rid of constants_unstruct.h.
-
-- checking use of real or double precision. Gain in elapsed time is ok, but now we have to look for memory consumption.
-
-- checking output on stdout (for data that should be printed only once).
-
- Hi Jeroen, Perfect. I think talking to Jean-Paul Ampuero would be useful
as well because in Utrecht last year he had told us that
he had implemented some nice 4th-order symplectic schemes
More information about the cig-commits
mailing list