From 47332f75c69e6a7b8e0d5a89507d6ffd00da59f5 Mon Sep 17 00:00:00 2001 From: Dimitri Komatitsch Date: Mon, 5 May 2014 02:59:23 +0200 Subject: [PATCH] switched from a single observation point to a whole surface around the Earth for Roland_Sylvain integrals --- setup/constants.h.in | 49 ++--- src/meshfem3D/compute_coordinates_grid.f90 | 161 ++++++++++++++-- ...reas.f90 => compute_volumes_and_areas.F90} | 173 +++++++++--------- src/meshfem3D/create_meshes.f90 | 12 +- src/meshfem3D/create_regions_mesh.F90 | 36 ++-- src/meshfem3D/finalize_mesher.f90 | 160 +++++++++++++--- src/meshfem3D/get_ellipticity.f90 | 36 ++++ src/meshfem3D/meshfem3D_models.f90 | 6 +- src/meshfem3D/meshfem3D_par.f90 | 46 +++-- src/shared/model_topo_bathy.f90 | 13 +- src/shared/parallel.f90 | 41 +++-- src/shared/rthetaphi_xyz.f90 | 8 +- src/specfem3D/check_stability.f90 | 3 +- 13 files changed, 542 insertions(+), 202 deletions(-) rename src/meshfem3D/{compute_volumes_and_areas.f90 => compute_volumes_and_areas.F90} (68%) diff --git a/setup/constants.h.in b/setup/constants.h.in index 8a9d6a98f..49eacd8e1 100644 --- a/setup/constants.h.in +++ b/setup/constants.h.in @@ -91,6 +91,8 @@ ! average density in the full Earth to normalize equation double precision, parameter :: RHOAV = 5514.3d0 +! maximum number of chunks (full sphere) + integer, parameter :: NCHUNKS_MAX = 6 !!----------------------------------------------------------- !! @@ -402,30 +404,35 @@ !! DK DK for Roland_Sylvain logical, parameter :: ROLAND_SYLVAIN = .false. -! altitude of the observation point - double precision, parameter :: altitude_of_observation_point = 255.d0 ! in km +! altitude of the observation points in km + double precision, parameter :: altitude_of_observation_points = 255.d0 -! standard gravity at the surface of the Earth - double precision, parameter :: STANDARD_GRAVITY_EARTH = 9.80665d0 ! in m.s-2 +! elevation ratio of the points at which we observe + double precision, parameter :: observation_elevation_ratio = (R_EARTH_KM + altitude_of_observation_points) / R_EARTH_KM + +! compute the contribution of the crust only (otherwise by default compute the contribution of the whole Earth) + logical, parameter :: COMPUTE_CRUST_CONTRIB_ONLY = .false. - double precision, parameter :: SQUARE_ROOT_OF_THREE = 1.73205080756887729352d0 +! check for negative Jacobians in the calculation of integrals or not +! (can safely be done once to check that the mesh is OK at a given resolution, and then permanently +! turned off in future runs because the mesh does not change) + logical, parameter :: CHECK_FOR_NEGATIVE_JACOBIANS = .true. -! position of the observation point in non-dimensionalized value +! number of points in each horizontal direction of the observation grid of each cubed-sphere chunk +! at the altitude of the observation point + integer, parameter :: NX_OBSERVATION = 500 + integer, parameter :: NY_OBSERVATION = NX_OBSERVATION -!! DK DK along X only -! double precision, parameter :: x_observation = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM -! double precision, parameter :: y_observation = 0.d0 -! double precision, parameter :: z_observation = 0.d0 +! the code will display sample output values at this particular point as a check + integer, parameter :: ixr = NX_OBSERVATION / 3 + integer, parameter :: iyr = NY_OBSERVATION / 4 + integer, parameter :: ichunkr = 3 -!! DK DK along diagonal -! double precision, parameter :: x_observation = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE -! double precision, parameter :: y_observation = x_observation -! double precision, parameter :: z_observation = x_observation +! how often (every how many spectral elements computed) we print a timestamp to monitor the behavior of the code + integer, parameter :: NSPEC_DISPLAY_INTERVAL = 100 -!! DK DK point randomly chosen in space, not at the right elevation we want - double precision, parameter :: x_observation = (R_EARTH_KM + altitude_of_observation_point) / R_EARTH_KM / SQUARE_ROOT_OF_THREE - double precision, parameter :: y_observation = x_observation * 1.12d0 - double precision, parameter :: z_observation = x_observation * 1.38d0 +! for the FORCE_VECTORIZATION 1D version of some loops + integer, parameter :: NTOTAL_OBSERVATION = NX_OBSERVATION * NY_OBSERVATION * NCHUNKS_MAX !------------------------------------------------------ !----------- do not modify anything below ------------- @@ -471,6 +478,9 @@ ! double precision, parameter :: GRAV = 6.6723d-11 double precision, parameter :: GRAV = 6.67384d-11 +! standard gravity at the surface of the Earth + double precision, parameter :: STANDARD_GRAVITY_EARTH = 9.80665d0 ! in m.s-2 + ! a few useful constants double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0 double precision, parameter :: ONE_HALF = HALF, ONE_FOURTH = 0.25d0, ONE_EIGHTH = 0.125d0, ONE_SIXTEENTH = 0.0625d0 @@ -505,9 +515,6 @@ ! Mountains can create signals of several hundred Eotvos. double precision, parameter :: SI_UNITS_TO_EOTVOS = 1.d+9 -! maximum number of chunks (full sphere) - integer, parameter :: NCHUNKS_MAX = 6 - ! define block type based upon chunk number (between 1 and 6) ! do not change this numbering, chunk AB must be number 1 for central cube integer, parameter :: CHUNK_AB = 1 diff --git a/src/meshfem3D/compute_coordinates_grid.f90 b/src/meshfem3D/compute_coordinates_grid.f90 index 9193759e8..7ddd5ab57 100644 --- a/src/meshfem3D/compute_coordinates_grid.f90 +++ b/src/meshfem3D/compute_coordinates_grid.f90 @@ -57,8 +57,7 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, & double precision, dimension(NDIM) :: vector_ori,vector_rotated - double precision :: ratio_xi, ratio_eta, fact_xi, fact_eta, & - fact_xi_,fact_eta_ + double precision :: ratio_xi, ratio_eta, fact_xi, fact_eta, fact_xi_,fact_eta_ ! this to avoid compilation warnings x_=0 @@ -82,8 +81,8 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, & ! uncomment the corresponding lines in the else condition of this if statement too. ! note that the ratio bigger_edge_size/smaller_edge_size for the surface mesh is a bit higher (1.43 vs 1.41) -! fact_xi_= (3.d0*fact_xi+4.d0*fact_xi_)/7.d0 -! fact_eta_= (3.d0*fact_eta+4.d0*fact_eta_)/7.d0 +! fact_xi_= (3.d0*fact_xi+4.d0*fact_xi_)/7.d0 +! fact_eta_= (3.d0*fact_eta+4.d0*fact_eta_)/7.d0 xi = PI_OVER_TWO*fact_xi eta = PI_OVER_TWO*fact_eta @@ -134,15 +133,15 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, & case default stop 'incorrect chunk number in compute_coord_main_mesh' end select -! write(IMAIN,*) x,' ',y,' ',z + else ! uncomment the following lines to have more regular surface mesh (better aspect ratio for each element) ! note that the ratio bigger_edge_size/smaller_edge_size for the surface mesh is a bit higher (1.43 vs 1.41) -! ratio_xi = ((iproc_xi + offset_x(ignod)/dble(NEX_PER_PROC_XI))/dble(NPROC_XI))*tan(ANGULAR_WIDTH_XI_RAD/2.d0) -! x_ = 2.d0*ratio_xi-tan(ANGULAR_WIDTH_XI_RAD/2.d0) -! ratio_eta = ((iproc_eta + offset_y(ignod)/dble(NEX_PER_PROC_ETA))/dble(NPROC_ETA))*tan(ANGULAR_WIDTH_ETA_RAD/2.d0) -! y_ = 2.d0*ratio_eta-tan(ANGULAR_WIDTH_ETA_RAD/2.d0) +! ratio_xi = ((iproc_xi + offset_x(ignod)/dble(NEX_PER_PROC_XI))/dble(NPROC_XI))*tan(ANGULAR_WIDTH_XI_RAD/2.d0) +! x_ = 2.d0*ratio_xi-tan(ANGULAR_WIDTH_XI_RAD/2.d0) +! ratio_eta = ((iproc_eta + offset_y(ignod)/dble(NEX_PER_PROC_ETA))/dble(NPROC_ETA))*tan(ANGULAR_WIDTH_ETA_RAD/2.d0) +! y_ = 2.d0*ratio_eta-tan(ANGULAR_WIDTH_ETA_RAD/2.d0) ratio_xi = ((iproc_xi + offset_x(ignod)/dble(NEX_PER_PROC_XI))/dble(NPROC_XI)) x = 2.d0*ratio_xi-1 @@ -155,8 +154,8 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, & ! uncomment the following lines to have more regular surface mesh (better aspect ratio for each element) ! note that the ratio bigger_edge_size/smaller_edge_size for the surface mesh is a bit higher (1.43 vs 1.41) -! x= (3.d0*x_+4.d0*x)/7.d0 -! y= (3.d0*y_+4.d0*y)/7.d0 +! x= (3.d0*x_+4.d0*x)/7.d0 +! y= (3.d0*y_+4.d0*y)/7.d0 gamma = ONE / sqrt(ONE + x*x + y*y) @@ -272,12 +271,12 @@ subroutine compute_coord_main_mesh(offset_x,offset_y,offset_z,xelm,yelm,zelm, & endif enddo -! if(ilayer == NUMBER_OF_MESH_LAYERS .and. INCLUDE_CENTRAL_CUBE) write(IMAIN,*) + end subroutine compute_coord_main_mesh !--------------------------------------------------------------------------- -!! DK DK create value of arrays xgrid ygrid and zgrid in the central cube without storing them +! create value of arrays xgrid ygrid and zgrid in the central cube without storing them subroutine compute_coord_central_cube(ix,iy,iz, & xgrid_central_cube,ygrid_central_cube,zgrid_central_cube, & @@ -321,3 +320,139 @@ subroutine compute_coord_central_cube(ix,iy,iz, & end subroutine compute_coord_central_cube + +!--------------------------------------------------------------------------- + +! create observation grid (surface) on which we compute the Roland_Sylvain integrals. +! each processor has a copy of the whole grid and sums its contribution there, +! and at the end of the process an MPI reduction is performed to compute the global sum + + subroutine compute_observation_surface() + + use constants + + use meshfem3D_par, only: myrank,x_observation,y_observation,z_observation,ELLIPTICITY,TOPOGRAPHY,OUTPUT_FILES + + use meshfem3D_models_par, only: nspl,rspl,espl,espl2,ibathy_topo + + implicit none + +! local variables + integer :: ix,iy,ichunk + + double precision :: gamma,x,y,rgt + double precision :: x_top,y_top,z_top + double precision :: ratio_xi, ratio_eta + + double precision :: r,lat,lon,elevation + +! the observation surface is always above the whole Earth, thus each mesh chunk has a size of PI / 2 in each direction + double precision, parameter :: ANGULAR_WIDTH_XI_RAD = PI_OVER_TWO, ANGULAR_WIDTH_ETA_RAD = PI_OVER_TWO + + ! for future GMT display + if(myrank == 0) open(unit=IOUT,file=trim(OUTPUT_FILES)//'/observation_grid_long_lat_topo_for_GMT.txt', & + status='unknown',action='write') + +! loop on all the chunks and then on all the observation nodes in each chunk + do ichunk = 1,NCHUNKS_MAX + do iy = 1,NY_OBSERVATION + do ix = 1,NX_OBSERVATION + + ratio_xi = dble(ix - 1) / dble(NX_OBSERVATION - 1) + x = 2.d0*ratio_xi-1 + + ratio_eta = dble(iy - 1) / dble(NY_OBSERVATION - 1) + y = 2.d0*ratio_eta-1 + + x = tan((ANGULAR_WIDTH_XI_RAD/2.d0) * x) + y = tan((ANGULAR_WIDTH_ETA_RAD/2.d0) * y) + + gamma = ONE / sqrt(ONE + x*x + y*y) + +! first compute the position of the points exactly at the free surface of the Earth (without the oceans) +! keeping in mind that the code non-dimensionalizes the radius of the spherical Earth to one + rgt = R_UNIT_SPHERE*gamma + + ! define the mesh points on the top and the bottom in the six regions of the cubed shpere + select case (ichunk) + + case(CHUNK_AB) + + x_top = -y*rgt + y_top = x*rgt + z_top = rgt + + case(CHUNK_AB_ANTIPODE) + + x_top = -y*rgt + y_top = -x*rgt + z_top = -rgt + + case(CHUNK_AC) + + x_top = -y*rgt + y_top = -rgt + z_top = x*rgt + + case(CHUNK_AC_ANTIPODE) + + x_top = -y*rgt + y_top = rgt + z_top = -x*rgt + + case(CHUNK_BC) + + x_top = -rgt + y_top = y*rgt + z_top = x*rgt + + case(CHUNK_BC_ANTIPODE) + + x_top = rgt + y_top = -y*rgt + z_top = x*rgt + + case default + stop 'incorrect chunk number in compute_coord_main_mesh' + + end select + + ! add ellipticity + if(ELLIPTICITY) call get_ellipticity_single_point(x_top,y_top,z_top,nspl,rspl,espl,espl2) + + ! add topography + if(TOPOGRAPHY) then + ! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees) + call xyz_2_rlatlon_dble(x_top,y_top,z_top,r,lat,lon) + + ! compute elevation at current point + call get_topo_bathy(lat,lon,elevation,ibathy_topo) + + ! store the values obtained for future display with GMT + if(myrank == 0) then + if( lon > 180.0d0 ) lon = lon - 360.0d0 + write(IOUT,*) lon,lat,elevation + endif + + ! non-dimensionalize the elevation, which is in meters + elevation = elevation / R_EARTH + + x_top = x_top*(ONE + elevation) + y_top = y_top*(ONE + elevation) + z_top = z_top*(ONE + elevation) + endif + + ! compute the position of the point + x_observation(ix,iy,ichunk) = x_top * observation_elevation_ratio + y_observation(ix,iy,ichunk) = y_top * observation_elevation_ratio + z_observation(ix,iy,ichunk) = z_top * observation_elevation_ratio + + enddo + enddo + enddo + + ! for future GMT display + if(myrank == 0) close(unit=IOUT) + + end subroutine compute_observation_surface + diff --git a/src/meshfem3D/compute_volumes_and_areas.f90 b/src/meshfem3D/compute_volumes_and_areas.F90 similarity index 68% rename from src/meshfem3D/compute_volumes_and_areas.f90 rename to src/meshfem3D/compute_volumes_and_areas.F90 index ef05c9400..5f5344a4e 100644 --- a/src/meshfem3D/compute_volumes_and_areas.f90 +++ b/src/meshfem3D/compute_volumes_and_areas.F90 @@ -259,18 +259,17 @@ end subroutine compute_Earth_mass ! compute Roland_Sylvain integrals of that part of the slice, and then total integrals for the whole Earth - subroutine compute_Roland_Sylvain_integr(myrank,Roland_Sylvain_integr_total, & - nspec,wxgll,wygll,wzgll,xstore,ystore,zstore,xixstore,xiystore,xizstore, & - etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling) + subroutine compute_Roland_Sylvain_integr(myrank,iregion_code,nspec,wxgll,wygll,wzgll,xstore,ystore,zstore, & + xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling) use constants - implicit none + use meshfem3D_par,only: x_observation,y_observation,z_observation,g_x,g_y,g_z,G_xx,G_yy,G_zz,G_xy,G_xz,G_yz, & + x_observation1D,y_observation1D,z_observation1D,g_x1D,g_y1D,g_z1D,G_xx1D,G_yy1D,G_zz1D,G_xy1D,G_xz1D,G_yz1D,OUTPUT_FILES - double precision, dimension(9) :: Roland_Sylvain_integr_total + implicit none - integer :: myrank - integer :: nspec + integer :: myrank,iregion_code,nspec double precision :: wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ) integer,dimension(nspec) :: idoubling @@ -289,34 +288,40 @@ subroutine compute_Roland_Sylvain_integr(myrank,Roland_Sylvain_integr_total, & integer :: i,j,k,ispec double precision :: xval,yval,zval double precision :: xval_squared,yval_squared,zval_squared - double precision :: x_meshpoint,y_meshpoint,z_meshpoint,rho_meshpoint - double precision :: distance,distance_squared,distance_cubed,distance_fifth_power, & + double precision :: x_meshpoint,y_meshpoint,z_meshpoint + double precision :: distance_squared,distance_cubed, & three_over_distance_squared,one_over_distance_cubed,three_over_distance_fifth_power - double precision :: common_multiplying_factor + double precision :: common_multiplying_factor,common_mult_times_one_over,common_mult_times_three_over - double precision, dimension(9) :: Roland_Sylvain_int_local,Roland_Sylvain_int_total_region - double precision :: elemental_contribution_1,elemental_contribution_2,elemental_contribution_3, & - elemental_contribution_4,elemental_contribution_5,elemental_contribution_6, & - elemental_contribution_7,elemental_contribution_8,elemental_contribution_9 + ! name of the timestamp files + character(len=150) :: outputname - ! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized - ! for the gravity vector force, a distance is involved in the dimensions - double precision, parameter :: nondimensionalizing_factor_gi = RHOAV * R_EARTH - ! for the second-order gravity tensor, no distance is involved in the dimensions - double precision, parameter :: nondimensionalizing_factor_Gij = RHOAV - - double precision, parameter :: scaling_factor_gi = GRAV * nondimensionalizing_factor_gi - double precision, parameter :: scaling_factor_Gij = GRAV * nondimensionalizing_factor_Gij - - ! initializes - Roland_Sylvain_int_local(:) = ZERO +#ifdef FORCE_VECTORIZATION + integer :: ix_iy_ichunk +#else + integer :: ix,iy,ichunk +#endif ! calculates volume of all elements in mesh do ispec = 1,nspec + ! print information about number of elements done so far + if(myrank == 0 .and. mod(ispec,NSPEC_DISPLAY_INTERVAL) == 0) then + write(IMAIN,*) 'for Roland_Sylvain integrals ',ispec,' elements computed out of ',nspec + ! write time stamp file to give information about progression of simulation + write(outputname,"('/timestamp_reg',i1.1,'_ispec',i7.7,'_out_of_',i7.7)") iregion_code,ispec,nspec + ! timestamp file output + open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',action='write') + write(IOUT,*) ispec,' elements done out of ',nspec,' in region ',iregion_code + close(unit=IOUT) + endif + ! suppress fictitious elements in central cube if(idoubling(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle + ! see if we compute the contribution of the crust only + if(COMPUTE_CRUST_CONTRIB_ONLY .and. idoubling(ispec) /= IFLAG_CRUST) cycle + do k = 1,NGLLZ do j = 1,NGLLY do i = 1,NGLLX @@ -334,76 +339,101 @@ subroutine compute_Roland_Sylvain_integr(myrank,Roland_Sylvain_integr_total, & gammayl = gammaystore(i,j,k,ispec) gammazl = gammazstore(i,j,k,ispec) -!! DK DK do this in double precision for accuracy + ! do this in double precision for accuracy jacobianl = 1.d0 / dble(xixl*(etayl*gammazl-etazl*gammayl) & - xiyl*(etaxl*gammazl-etazl*gammaxl) & + xizl*(etaxl*gammayl-etayl*gammaxl)) + if(CHECK_FOR_NEGATIVE_JACOBIANS .and. jacobianl <= ZERO) stop 'error: negative Jacobian found in integral calculation' + x_meshpoint = xstore(i,j,k,ispec) y_meshpoint = ystore(i,j,k,ispec) z_meshpoint = zstore(i,j,k,ispec) - rho_meshpoint = rhostore(i,j,k,ispec) + common_multiplying_factor = jacobianl * weight * rhostore(i,j,k,ispec) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!! beginning of loop on all the data to create !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - xval = x_meshpoint - x_observation - yval = y_meshpoint - y_observation - zval = z_meshpoint - z_observation +! loop on all the chunks and then on all the observation nodes in each chunk + +#ifdef FORCE_VECTORIZATION +! this works only if the arrays are contiguous in memory (which is always the case for static arrays, as used in the code) +! but the code produced is extremely fast because we get a single and fully-vectorized loop on the whole array + do ix_iy_ichunk = 1,NTOTAL_OBSERVATION + + xval = x_meshpoint - x_observation1D(ix_iy_ichunk) + yval = y_meshpoint - y_observation1D(ix_iy_ichunk) + zval = z_meshpoint - z_observation1D(ix_iy_ichunk) xval_squared = xval**2 yval_squared = yval**2 zval_squared = zval**2 distance_squared = xval_squared + yval_squared + zval_squared - distance = sqrt(distance_squared) - distance_cubed = distance_squared*distance - distance_fifth_power = distance_squared*distance_cubed + distance_cubed = distance_squared * sqrt(distance_squared) three_over_distance_squared = 3.d0 / distance_squared one_over_distance_cubed = 1.d0 / distance_cubed three_over_distance_fifth_power = three_over_distance_squared * one_over_distance_cubed -! g_x - elemental_contribution_1 = xval * one_over_distance_cubed + common_mult_times_one_over = common_multiplying_factor * one_over_distance_cubed + common_mult_times_three_over = common_multiplying_factor * three_over_distance_fifth_power -! g_y - elemental_contribution_2 = yval * one_over_distance_cubed + g_x1D(ix_iy_ichunk) = g_x1D(ix_iy_ichunk) + common_mult_times_one_over * xval + g_y1D(ix_iy_ichunk) = g_y1D(ix_iy_ichunk) + common_mult_times_one_over * yval + g_z1D(ix_iy_ichunk) = g_z1D(ix_iy_ichunk) + common_mult_times_one_over * zval -! g_z - elemental_contribution_3 = zval * one_over_distance_cubed + G_xx1D(ix_iy_ichunk) = G_xx1D(ix_iy_ichunk) + common_mult_times_one_over * (xval_squared * three_over_distance_squared - 1.d0) + G_yy1D(ix_iy_ichunk) = G_yy1D(ix_iy_ichunk) + common_mult_times_one_over * (yval_squared * three_over_distance_squared - 1.d0) + G_zz1D(ix_iy_ichunk) = G_zz1D(ix_iy_ichunk) + common_mult_times_one_over * (zval_squared * three_over_distance_squared - 1.d0) -! G_xx - elemental_contribution_4 = (xval_squared * three_over_distance_squared - 1.d0) * one_over_distance_cubed + G_xy1D(ix_iy_ichunk) = G_xy1D(ix_iy_ichunk) + common_mult_times_three_over * xval*yval + G_xz1D(ix_iy_ichunk) = G_xz1D(ix_iy_ichunk) + common_mult_times_three_over * xval*zval + G_yz1D(ix_iy_ichunk) = G_yz1D(ix_iy_ichunk) + common_mult_times_three_over * yval*zval -! G_yy - elemental_contribution_5 = (yval_squared * three_over_distance_squared - 1.d0) * one_over_distance_cubed + enddo -! G_zz - elemental_contribution_6 = (zval_squared * three_over_distance_squared - 1.d0) * one_over_distance_cubed +#else + do ichunk = 1,NCHUNKS_MAX + do iy = 1,NY_OBSERVATION + do ix = 1,NX_OBSERVATION -! G_xy - elemental_contribution_7 = xval*yval * three_over_distance_fifth_power + xval = x_meshpoint - x_observation(ix,iy,ichunk) + yval = y_meshpoint - y_observation(ix,iy,ichunk) + zval = z_meshpoint - z_observation(ix,iy,ichunk) -! G_xz - elemental_contribution_8 = xval*zval * three_over_distance_fifth_power + xval_squared = xval**2 + yval_squared = yval**2 + zval_squared = zval**2 + + distance_squared = xval_squared + yval_squared + zval_squared + distance_cubed = distance_squared * sqrt(distance_squared) + + three_over_distance_squared = 3.d0 / distance_squared + one_over_distance_cubed = 1.d0 / distance_cubed + three_over_distance_fifth_power = three_over_distance_squared * one_over_distance_cubed -! G_yz - elemental_contribution_9 = yval*zval * three_over_distance_fifth_power + common_mult_times_one_over = common_multiplying_factor * one_over_distance_cubed + common_mult_times_three_over = common_multiplying_factor * three_over_distance_fifth_power - common_multiplying_factor = jacobianl * weight * rho_meshpoint + g_x(ix,iy,ichunk) = g_x(ix,iy,ichunk) + common_mult_times_one_over * xval + g_y(ix,iy,ichunk) = g_y(ix,iy,ichunk) + common_mult_times_one_over * yval + g_z(ix,iy,ichunk) = g_z(ix,iy,ichunk) + common_mult_times_one_over * zval - Roland_Sylvain_int_local(1) = Roland_Sylvain_int_local(1) + common_multiplying_factor*elemental_contribution_1 - Roland_Sylvain_int_local(2) = Roland_Sylvain_int_local(2) + common_multiplying_factor*elemental_contribution_2 - Roland_Sylvain_int_local(3) = Roland_Sylvain_int_local(3) + common_multiplying_factor*elemental_contribution_3 - Roland_Sylvain_int_local(4) = Roland_Sylvain_int_local(4) + common_multiplying_factor*elemental_contribution_4 - Roland_Sylvain_int_local(5) = Roland_Sylvain_int_local(5) + common_multiplying_factor*elemental_contribution_5 - Roland_Sylvain_int_local(6) = Roland_Sylvain_int_local(6) + common_multiplying_factor*elemental_contribution_6 - Roland_Sylvain_int_local(7) = Roland_Sylvain_int_local(7) + common_multiplying_factor*elemental_contribution_7 - Roland_Sylvain_int_local(8) = Roland_Sylvain_int_local(8) + common_multiplying_factor*elemental_contribution_8 - Roland_Sylvain_int_local(9) = Roland_Sylvain_int_local(9) + common_multiplying_factor*elemental_contribution_9 + G_xx(ix,iy,ichunk) = G_xx(ix,iy,ichunk) + common_mult_times_one_over * (xval_squared * three_over_distance_squared - 1.d0) + G_yy(ix,iy,ichunk) = G_yy(ix,iy,ichunk) + common_mult_times_one_over * (yval_squared * three_over_distance_squared - 1.d0) + G_zz(ix,iy,ichunk) = G_zz(ix,iy,ichunk) + common_mult_times_one_over * (zval_squared * three_over_distance_squared - 1.d0) + + G_xy(ix,iy,ichunk) = G_xy(ix,iy,ichunk) + common_mult_times_three_over * xval*yval + G_xz(ix,iy,ichunk) = G_xz(ix,iy,ichunk) + common_mult_times_three_over * xval*zval + G_yz(ix,iy,ichunk) = G_yz(ix,iy,ichunk) + common_mult_times_three_over * yval*zval + + enddo + enddo + enddo +#endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!! end of loop on all the data to create @@ -414,26 +444,5 @@ subroutine compute_Roland_Sylvain_integr(myrank,Roland_Sylvain_integr_total, & enddo enddo - ! multiply by the gravitational constant in S.I. units i.e. in m3 kg-1 s-2 - ! and also take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized - Roland_Sylvain_int_local(1:3) = Roland_Sylvain_int_local(1:3) * scaling_factor_gi - Roland_Sylvain_int_local(4:9) = Roland_Sylvain_int_local(4:9) * scaling_factor_Gij - - ! use an MPI reduction to compute the total value of the integral - Roland_Sylvain_int_total_region(:) = ZERO -!! DK DK could use a single MPI call for the nine values - call sum_all_dp(Roland_Sylvain_int_local(1),Roland_Sylvain_int_total_region(1)) - call sum_all_dp(Roland_Sylvain_int_local(2),Roland_Sylvain_int_total_region(2)) - call sum_all_dp(Roland_Sylvain_int_local(3),Roland_Sylvain_int_total_region(3)) - call sum_all_dp(Roland_Sylvain_int_local(4),Roland_Sylvain_int_total_region(4)) - call sum_all_dp(Roland_Sylvain_int_local(5),Roland_Sylvain_int_total_region(5)) - call sum_all_dp(Roland_Sylvain_int_local(6),Roland_Sylvain_int_total_region(6)) - call sum_all_dp(Roland_Sylvain_int_local(7),Roland_Sylvain_int_total_region(7)) - call sum_all_dp(Roland_Sylvain_int_local(8),Roland_Sylvain_int_total_region(8)) - call sum_all_dp(Roland_Sylvain_int_local(9),Roland_Sylvain_int_total_region(9)) - - ! sum volume over all the regions - if(myrank == 0) Roland_Sylvain_integr_total(:) = Roland_Sylvain_integr_total(:) + Roland_Sylvain_int_total_region(:) - end subroutine compute_Roland_Sylvain_integr diff --git a/src/meshfem3D/create_meshes.f90 b/src/meshfem3D/create_meshes.f90 index c159c96e6..9c4e4c391 100644 --- a/src/meshfem3D/create_meshes.f90 +++ b/src/meshfem3D/create_meshes.f90 @@ -43,7 +43,17 @@ subroutine create_meshes() ! and Roland_Sylvain integrals volume_total = ZERO Earth_mass_total = ZERO - Roland_Sylvain_integr_total(:) = ZERO + + g_x(:,:,:) = ZERO + g_y(:,:,:) = ZERO + g_z(:,:,:) = ZERO + + G_xx(:,:,:) = ZERO + G_yy(:,:,:) = ZERO + G_zz(:,:,:) = ZERO + G_xy(:,:,:) = ZERO + G_xz(:,:,:) = ZERO + G_yz(:,:,:) = ZERO ! make sure everybody is synchronized call synchronize_all() diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90 index 9d3a61399..080a04236 100644 --- a/src/meshfem3D/create_regions_mesh.F90 +++ b/src/meshfem3D/create_regions_mesh.F90 @@ -47,7 +47,7 @@ subroutine create_regions_mesh(iregion_code, & use meshfem3D_par,only: & ibool,idoubling,xstore,ystore,zstore, & - IMAIN,volume_total,Earth_mass_total,Roland_Sylvain_integr_total,myrank,LOCAL_PATH, & + IMAIN,volume_total,Earth_mass_total,myrank,LOCAL_PATH, & IREGION_CRUST_MANTLE,IREGION_OUTER_CORE,IREGION_INNER_CORE, & IFLAG_IN_FICTITIOUS_CUBE, & NCHUNKS,SAVE_MESH_FILES,ABSORBING_CONDITIONS, & @@ -152,9 +152,14 @@ subroutine create_regions_mesh(iregion_code, & offset_proc_xi,offset_proc_eta) - ! only create global addressing and the MPI buffers in the first pass select case(ipass) - case( 1 ) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + case( 1 ) !!!!!!!!!!! first pass of the mesher +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! only create global addressing and the MPI buffers in the first pass + ! creates ibool index array for projection from local to global points call synchronize_all() if( myrank == 0) then @@ -176,14 +181,16 @@ subroutine create_regions_mesh(iregion_code, & ! sets up Stacey absorbing boundary indices (nimin,nimax,..) - if(NCHUNKS /= 6) then - call get_absorb(myrank,prname,iregion_code, iboun,nspec,nimin,nimax,& - njmin,njmax, nkmin_xi,nkmin_eta, NSPEC2DMAX_XMIN_XMAX, & - NSPEC2DMAX_YMIN_YMAX, NSPEC2D_BOTTOM) - endif + if(NCHUNKS /= 6) call get_absorb(myrank,prname,iregion_code, iboun,nspec,nimin,nimax,& + njmin,njmax, nkmin_xi,nkmin_eta, NSPEC2DMAX_XMIN_XMAX, & + NSPEC2DMAX_YMIN_YMAX, NSPEC2D_BOTTOM) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + case( 2 ) !!!!!!!!!!! second pass of the mesher +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! only create mass matrix and save all the final arrays in the second pass - case( 2 ) + ! precomputes Jacobian for 2D absorbing boundary surfaces call synchronize_all() if( myrank == 0) then @@ -205,6 +212,10 @@ subroutine create_regions_mesh(iregion_code, & NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,& xigll,yigll,zigll) +!! DK DK for Roland_Sylvain +! creation of the top observation surface if region is the crust_mantle + if(ROLAND_SYLVAIN .and. iregion_code == IREGION_CRUST_MANTLE) call compute_observation_surface() + ! create chunk buffers if more than one chunk call synchronize_all() if( myrank == 0) then @@ -345,7 +356,7 @@ subroutine create_regions_mesh(iregion_code, & if(ier /= 0) stop 'error in allocate 22' ! creating mass matrices in this slice (will be fully assembled in the solver) - ! note: for Stacey boundaries, needs indexing nimin,.. filled in in first pass + ! note: for Stacey boundaries, needs indexing nimin,.. filled in the first pass call create_mass_matrices(nspec,nglob,idoubling,ibool, & iregion_code,xstore,ystore,zstore, & NSPEC2D_TOP,NSPEC2D_BOTTOM) @@ -428,9 +439,8 @@ subroutine create_regions_mesh(iregion_code, & endif ! compute Roland_Sylvain integrals of that part of the slice, and then total integrals for the whole Earth - call compute_Roland_Sylvain_integr(myrank,Roland_Sylvain_integr_total, & - nspec,wxgll,wygll,wzgll,xstore,ystore,zstore,xixstore,xiystore,xizstore, & - etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling) + call compute_Roland_Sylvain_integr(myrank,iregion_code,nspec,wxgll,wygll,wzgll,xstore,ystore,zstore, & + xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,rhostore,idoubling) endif diff --git a/src/meshfem3D/finalize_mesher.f90 b/src/meshfem3D/finalize_mesher.f90 index 7fc4f47e8..e04623b25 100644 --- a/src/meshfem3D/finalize_mesher.f90 +++ b/src/meshfem3D/finalize_mesher.f90 @@ -38,7 +38,66 @@ subroutine finalize_mesher() double precision, external :: wtime ! for Roland_Sylvain integrals - double precision :: g_x, g_y, g_z, G_xx, G_yy, G_zz, G_xy, G_xz, G_yz, real_altitude_of_observ_point + ! take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized + ! for the gravity vector force, a distance is involved in the dimensions + double precision, parameter :: nondimensionalizing_factor_gi = RHOAV * R_EARTH + ! for the second-order gravity tensor, no distance is involved in the dimensions + double precision, parameter :: nondimensionalizing_factor_Gij = RHOAV + + double precision, parameter :: scaling_factor_gi = GRAV * nondimensionalizing_factor_gi + double precision, parameter :: scaling_factor_Gij_Eotvos = GRAV * nondimensionalizing_factor_Gij * SI_UNITS_TO_EOTVOS + + double precision :: real_altitude_of_observ_point + + integer :: ixval,iyval,ichunkval + + if(ROLAND_SYLVAIN) then + + ! multiply by the gravitational constant in S.I. units i.e. in m3 kg-1 s-2 + ! and also take into account the fact that the density and the radius of the Earth have previously been non-dimensionalized + ! the final result is in m.s-2 i.e. in S.I. units + g_x(:,:,:) = g_x(:,:,:) * scaling_factor_gi + g_y(:,:,:) = g_y(:,:,:) * scaling_factor_gi + g_z(:,:,:) = g_z(:,:,:) * scaling_factor_gi + + ! the final result is in Eotvos = 1.e+9 s-2 + G_xx(:,:,:) = G_xx(:,:,:) * scaling_factor_Gij_Eotvos + G_yy(:,:,:) = G_yy(:,:,:) * scaling_factor_Gij_Eotvos + G_zz(:,:,:) = G_zz(:,:,:) * scaling_factor_Gij_Eotvos + G_xy(:,:,:) = G_xy(:,:,:) * scaling_factor_Gij_Eotvos + G_xz(:,:,:) = G_xz(:,:,:) * scaling_factor_Gij_Eotvos + G_yz(:,:,:) = G_yz(:,:,:) * scaling_factor_Gij_Eotvos + + ! use an MPI reduction to compute the total value of the integral into a temporary array + ! and then copy it back into the original array + call sum_all_3Darray_dp(g_x,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) g_x(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(g_y,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) g_y(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(g_z,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) g_z(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(G_xx,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) G_xx(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(G_yy,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) G_yy(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(G_zz,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) G_zz(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(G_xy,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) G_xy(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(G_xz,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) G_xz(:,:,:) = temporary_array_for_sum(:,:,:) + + call sum_all_3Darray_dp(G_yz,temporary_array_for_sum,NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) + if(myrank == 0) G_yz(:,:,:) = temporary_array_for_sum(:,:,:) + + endif !--- print number of points and elements in the mesh for each region if(myrank == 0) then @@ -61,47 +120,94 @@ subroutine finalize_mesher() endif !! DK DK for Roland_Sylvain - ! Roland_Sylvain integrals if(ROLAND_SYLVAIN) then -! in m.s-2 - g_x = Roland_Sylvain_integr_total(1) - g_y = Roland_Sylvain_integr_total(2) - g_z = Roland_Sylvain_integr_total(3) + temporary_array_for_sum(:,:,:) = sqrt(g_x(:,:,:)**2 + g_y(:,:,:)**2 + g_z(:,:,:)**2) + write(IMAIN,*) + write(IMAIN,*) 'minval of norm of g vector on whole observation surface = ',minval(temporary_array_for_sum),' m.s-2' + write(IMAIN,*) 'maxval of norm of g vector on whole observation surface = ',maxval(temporary_array_for_sum),' m.s-2' -! in Eotvos = 1.e+9 s-2 - G_xx = Roland_Sylvain_integr_total(4) * SI_UNITS_TO_EOTVOS - G_yy = Roland_Sylvain_integr_total(5) * SI_UNITS_TO_EOTVOS - G_zz = Roland_Sylvain_integr_total(6) * SI_UNITS_TO_EOTVOS - G_xy = Roland_Sylvain_integr_total(7) * SI_UNITS_TO_EOTVOS - G_xz = Roland_Sylvain_integr_total(8) * SI_UNITS_TO_EOTVOS - G_yz = Roland_Sylvain_integr_total(9) * SI_UNITS_TO_EOTVOS + write(IMAIN,*) + write(IMAIN,*) 'minval of abs(G_xx) on whole observation surface = ',minval(abs(G_xx)),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_xx) on whole observation surface = ',maxval(abs(G_xx)),' Eotvos' write(IMAIN,*) - write(IMAIN,*) 'computed total Roland_Sylvain integral g_x = ',g_x,' m.s-2' - write(IMAIN,*) 'computed total Roland_Sylvain integral g_y = ',g_y,' m.s-2' - write(IMAIN,*) 'computed total Roland_Sylvain integral g_z = ',g_z,' m.s-2' + write(IMAIN,*) 'minval of abs(G_yy) on whole observation surface = ',minval(abs(G_yy)),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_yy) on whole observation surface = ',maxval(abs(G_yy)),' Eotvos' + + write(IMAIN,*) + write(IMAIN,*) 'minval of abs(G_zz) on whole observation surface = ',minval(abs(G_zz)),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_zz) on whole observation surface = ',maxval(abs(G_zz)),' Eotvos' + + write(IMAIN,*) + write(IMAIN,*) 'minval of abs(G_xy) on whole observation surface = ',minval(abs(G_xy)),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_xy) on whole observation surface = ',maxval(abs(G_xy)),' Eotvos' + + write(IMAIN,*) + write(IMAIN,*) 'minval of abs(G_xz) on whole observation surface = ',minval(abs(G_xz)),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_xz) on whole observation surface = ',maxval(abs(G_xz)),' Eotvos' + + write(IMAIN,*) + write(IMAIN,*) 'minval of abs(G_yz) on whole observation surface = ',minval(abs(G_yz)),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_yz) on whole observation surface = ',maxval(abs(G_yz)),' Eotvos' + + write(IMAIN,*) + write(IMAIN,*) 'Minval and maxval of trace of G, which in principle should be zero:' + write(IMAIN,*) + temporary_array_for_sum(:,:,:) = abs(G_xx(:,:,:) + G_yy(:,:,:) + G_zz(:,:,:)) + write(IMAIN,*) 'minval of abs(G_xx + G_yy + G_zz) on whole observation surface = ',minval(temporary_array_for_sum),' Eotvos' + write(IMAIN,*) 'maxval of abs(G_xx + G_yy + G_zz) on whole observation surface = ',maxval(temporary_array_for_sum),' Eotvos' + + write(IMAIN,*) + write(IMAIN,*) '-----------------------------' + write(IMAIN,*) + write(IMAIN,*) 'displaying the fields computed at:' + write(IMAIN,*) ' ix_observation = ',ixr,' out of ',NX_OBSERVATION + write(IMAIN,*) ' iy_observation = ',iyr,' out of ',NY_OBSERVATION + write(IMAIN,*) ' of mesh chunk ',ichunkr + write(IMAIN,*) + write(IMAIN,*) 'computed g_x = ',g_x(ixr,iyr,ichunkr),' m.s-2' + write(IMAIN,*) 'computed g_y = ',g_y(ixr,iyr,ichunkr),' m.s-2' + write(IMAIN,*) 'computed g_z = ',g_z(ixr,iyr,ichunkr),' m.s-2' write(IMAIN,*) - write(IMAIN,*) 'computed norm of g vector = ',sqrt(g_x**2 + g_y**2 + g_z**2),' m.s-2' + write(IMAIN,*) 'computed norm of g vector = ',sqrt(g_x(ixr,iyr,ichunkr)**2 + g_y(ixr,iyr,ichunkr)**2 + & + g_z(ixr,iyr,ichunkr)**2),' m.s-2' - real_altitude_of_observ_point = sqrt(x_observation**2 + y_observation**2 + z_observation**2) + real_altitude_of_observ_point = sqrt(x_observation(ixr,iyr,ichunkr)**2 + y_observation(ixr,iyr,ichunkr)**2 + & + z_observation(ixr,iyr,ichunkr)**2) ! gravity force vector norm decays approximately as (r / r_prime)^2 above the surface of the Earth write(IMAIN,*) ' (should be not too far from ', & sngl(STANDARD_GRAVITY_EARTH * (R_UNIT_SPHERE / real_altitude_of_observ_point)**2),' m.s-2)' write(IMAIN,*) - write(IMAIN,*) 'computed total Roland_Sylvain integral G_xx = ',G_xx,' Eotvos' - write(IMAIN,*) 'computed total Roland_Sylvain integral G_yy = ',G_yy,' Eotvos' - write(IMAIN,*) 'computed total Roland_Sylvain integral G_zz = ',G_zz,' Eotvos' + write(IMAIN,*) 'computed G_xx = ',G_xx(ixr,iyr,ichunkr),' Eotvos' + write(IMAIN,*) 'computed G_yy = ',G_yy(ixr,iyr,ichunkr),' Eotvos' + write(IMAIN,*) 'computed G_zz = ',G_zz(ixr,iyr,ichunkr),' Eotvos' write(IMAIN,*) write(IMAIN,*) 'G tensor should be traceless, G_xx + G_yy + G_zz = 0.' - write(IMAIN,*) 'Actual sum obtained = ',G_xx + G_yy + G_zz - if(max(abs(G_xx),abs(G_yy),abs(G_zz)) > TINYVAL) write(IMAIN,*) ' i.e., ', & - sngl(100.d0*abs(G_xx + G_yy + G_zz) / max(abs(G_xx),abs(G_yy),abs(G_zz))),'% of max(abs(Gxx),abs(Gyy),abs(Gzz))' + write(IMAIN,*) 'Actual sum obtained = ',G_xx(ixr,iyr,ichunkr) + G_yy(ixr,iyr,ichunkr) + G_zz(ixr,iyr,ichunkr) + if(max(abs(G_xx(ixr,iyr,ichunkr)),abs(G_yy(ixr,iyr,ichunkr)),abs(G_zz(ixr,iyr,ichunkr))) > TINYVAL) & + write(IMAIN,*) ' i.e., ',sngl(100.d0*abs(G_xx(ixr,iyr,ichunkr) + G_yy(ixr,iyr,ichunkr) + G_zz(ixr,iyr,ichunkr)) / & + max(abs(G_xx(ixr,iyr,ichunkr)),abs(G_yy(ixr,iyr,ichunkr)),abs(G_zz(ixr,iyr,ichunkr)))), & + '% of max(abs(G_xx),abs(G_yy),abs(G_zz))' write(IMAIN,*) - write(IMAIN,*) 'computed total Roland_Sylvain integral G_xy = ',G_xy,' Eotvos' - write(IMAIN,*) 'computed total Roland_Sylvain integral G_xz = ',G_xz,' Eotvos' - write(IMAIN,*) 'computed total Roland_Sylvain integral G_yz = ',G_yz,' Eotvos' + write(IMAIN,*) 'computed G_xy = ',G_xy(ixr,iyr,ichunkr),' Eotvos' + write(IMAIN,*) 'computed G_xz = ',G_xz(ixr,iyr,ichunkr),' Eotvos' + write(IMAIN,*) 'computed G_yz = ',G_yz(ixr,iyr,ichunkr),' Eotvos' + + ! for future GMT display + open(unit=IOUT,file=trim(OUTPUT_FILES)//'/results_g_and_G_for_GMT.txt',status='unknown',action='write') + ! loop on all the chunks and then on all the observation nodes in each chunk + do ichunkval = 1,NCHUNKS_MAX + do iyval = 1,NY_OBSERVATION + do ixval = 1,NX_OBSERVATION + write(IOUT,*) g_x(ixval,iyval,ichunkval),g_y(ixval,iyval,ichunkval),g_z(ixval,iyval,ichunkval), & + G_xx(ixval,iyval,ichunkval),G_yy(ixval,iyval,ichunkval),G_zz(ixval,iyval,ichunkval), & + G_xy(ixval,iyval,ichunkval),G_xz(ixval,iyval,ichunkval),G_yz(ixval,iyval,ichunkval) + enddo + enddo + enddo + close(unit=IOUT) endif diff --git a/src/meshfem3D/get_ellipticity.f90 b/src/meshfem3D/get_ellipticity.f90 index cab10d5ba..8d51937ce 100644 --- a/src/meshfem3D/get_ellipticity.f90 +++ b/src/meshfem3D/get_ellipticity.f90 @@ -123,3 +123,39 @@ subroutine get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl,rspl,espl,e end subroutine get_ellipticity_gll + +! +!------------------------------------------------------------------------------------------------- +! + + subroutine get_ellipticity_single_point(x,y,z,nspl,rspl,espl,espl2) + + use constants + + implicit none + + integer :: nspl + double precision :: x,y,z + double precision :: rspl(NR),espl(NR),espl2(NR) + + ! local parameters + double precision :: ell + double precision :: r,theta,phi,factor + double precision :: cost,p20 + + call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi) + + cost = dcos(theta) + p20 = 0.5d0*(3.0d0*cost*cost-1.0d0) + + ! get ellipticity using spline evaluation + call spline_evaluation(rspl,espl,espl2,nspl,r,ell) + + factor = ONE-(TWO/3.0d0)*ell*p20 + + x = x*factor + y = y*factor + z = z*factor + + end subroutine get_ellipticity_single_point + diff --git a/src/meshfem3D/meshfem3D_models.f90 b/src/meshfem3D/meshfem3D_models.f90 index 250b0edaf..087e89515 100644 --- a/src/meshfem3D/meshfem3D_models.f90 +++ b/src/meshfem3D/meshfem3D_models.f90 @@ -62,12 +62,10 @@ subroutine meshfem3D_models_broadcast(myrank,NSPEC, & endif ! sets up spline coefficients for ellipticity - if(ELLIPTICITY) & - call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST) + if(ELLIPTICITY) call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST) ! read topography and bathymetry file - if(TOPOGRAPHY) & - call model_topo_bathy_broadcast(myrank,ibathy_topo,LOCAL_PATH) + if(TOPOGRAPHY) call model_topo_bathy_broadcast(myrank,ibathy_topo,LOCAL_PATH) ! reads 1D reference models ! re-defines/initializes models 1066a and ak135 and ref diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90 index 7bf097c06..dcdc49377 100644 --- a/src/meshfem3D/meshfem3D_par.f90 +++ b/src/meshfem3D/meshfem3D_par.f90 @@ -159,8 +159,28 @@ module meshfem3D_par ! check Earth mass computed in the final mesh double precision :: Earth_mass_total - ! compute Roland_Sylvain integrals in the final mesh - double precision, dimension(9) :: Roland_Sylvain_integr_total + ! arrays containing the positions of the observation points in non-dimensionalized value for Roland_Sylvain integrals + ! the 1D equivalenced versions are for the FORCE_VECTORIZATION version of the loops + double precision, dimension(NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) :: x_observation,y_observation,z_observation + double precision, dimension(NTOTAL_OBSERVATION) :: x_observation1D,y_observation1D,z_observation1D + equivalence(x_observation,x_observation1D) + equivalence(y_observation,y_observation1D) + equivalence(z_observation,z_observation1D) + + ! arrays containing the computed fields for Roland_Sylvain integrals at the observation points + ! the 1D equivalenced versions are for the FORCE_VECTORIZATION version of the loops + double precision, dimension(NX_OBSERVATION,NY_OBSERVATION,NCHUNKS_MAX) :: g_x,g_y,g_z,G_xx,G_yy,G_zz,G_xy,G_xz,G_yz, & + temporary_array_for_sum + double precision, dimension(NTOTAL_OBSERVATION) :: g_x1D,g_y1D,g_z1D,G_xx1D,G_yy1D,G_zz1D,G_xy1D,G_xz1D,G_yz1D + equivalence(g_x,g_x1D) + equivalence(g_y,g_y1D) + equivalence(g_z,g_z1D) + equivalence(G_xx,G_xx1D) + equivalence(G_yy,G_yy1D) + equivalence(G_zz,G_zz1D) + equivalence(G_xy,G_xy1D) + equivalence(G_xz,G_xz1D) + equivalence(G_yz,G_yz1D) ! for loop on all the slices integer :: iregion_code @@ -185,17 +205,17 @@ module meshfem3D_par double precision :: static_memory_size integer :: NSPECMAX_ANISO_IC,NSPECMAX_ISO_MANTLE,NSPECMAX_TISO_MANTLE, & - NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUATION, & - NSPEC_INNER_CORE_ATTENUATION, & - NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, & - NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, & - NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, & - NSPEC_CRUST_MANTLE_ADJOINT, & - NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, & - NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, & - NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, & - NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, & - NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION + NSPECMAX_ANISO_MANTLE,NSPEC_CRUST_MANTLE_ATTENUATION, & + NSPEC_INNER_CORE_ATTENUATION, & + NSPEC_CRUST_MANTLE_STR_OR_ATT,NSPEC_INNER_CORE_STR_OR_ATT, & + NSPEC_CRUST_MANTLE_STR_AND_ATT,NSPEC_INNER_CORE_STR_AND_ATT, & + NSPEC_CRUST_MANTLE_STRAIN_ONLY,NSPEC_INNER_CORE_STRAIN_ONLY, & + NSPEC_CRUST_MANTLE_ADJOINT, & + NSPEC_OUTER_CORE_ADJOINT,NSPEC_INNER_CORE_ADJOINT, & + NGLOB_CRUST_MANTLE_ADJOINT,NGLOB_OUTER_CORE_ADJOINT, & + NGLOB_INNER_CORE_ADJOINT,NSPEC_OUTER_CORE_ROT_ADJOINT, & + NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, & + NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION ! this for the different corners of the slice (which are different if the superbrick is cut) ! 1 : xi_min, eta_min diff --git a/src/shared/model_topo_bathy.f90 b/src/shared/model_topo_bathy.f90 index db473a13a..d13ec541c 100644 --- a/src/shared/model_topo_bathy.f90 +++ b/src/shared/model_topo_bathy.f90 @@ -56,7 +56,7 @@ subroutine model_topo_bathy_broadcast(myrank,ibathy_topo,LOCAL_PATH) ! read/save topo file on master call read_topo_bathy_file(ibathy_topo) - call save_topo_bathy_database(ibathy_topo,LOCAL_PATH) + if(.not. ROLAND_SYLVAIN) call save_topo_bathy_database(ibathy_topo,LOCAL_PATH) endif ! broadcast the information read on the master to the nodes @@ -176,6 +176,7 @@ subroutine save_topo_bathy_database(ibathy_topo,LOCAL_PATH) ! saves topography and bathymetry file for solver open(unit=27,file=prname(1:len_trim(prname))//'topo.bin', & status='unknown',form='unformatted',action='write',iostat=ier) + if( ier /= 0 ) then ! inform about missing database topo file print*,'TOPOGRAPHY problem:' @@ -214,6 +215,7 @@ subroutine read_topo_bathy_database(ibathy_topo,LOCAL_PATH) ! reads topography and bathymetry file from saved database file open(unit=27,file=prname(1:len_trim(prname))//'topo.bin', & status='unknown',form='unformatted',action='read',iostat=ier) + if( ier /= 0 ) then ! inform user print*,'TOPOGRAPHY problem:' @@ -226,7 +228,8 @@ subroutine read_topo_bathy_database(ibathy_topo,LOCAL_PATH) call read_topo_bathy_file(ibathy_topo) ! saves database topo file for next time - call save_topo_bathy_database(ibathy_topo,LOCAL_PATH) + if(.not. ROLAND_SYLVAIN) call save_topo_bathy_database(ibathy_topo,LOCAL_PATH) + else ! database topo file exists read(27) ibathy_topo @@ -235,6 +238,7 @@ subroutine read_topo_bathy_database(ibathy_topo,LOCAL_PATH) ! user output write(IMAIN,*) " topography/bathymetry: min/max = ",minval(ibathy_topo),maxval(ibathy_topo) call flush_IMAIN() + endif end subroutine read_topo_bathy_database @@ -310,12 +314,14 @@ subroutine get_topo_bathy(xlat,xlon,value,ibathy_topo) + dble(ibathy_topo(iel1+1,iadd1))*ratio_lon*(1.-ratio_lat) & + dble(ibathy_topo(iel1+1,iadd1+1))*ratio_lon*ratio_lat & + dble(ibathy_topo(iel1,iadd1+1))*(1.-ratio_lon)*ratio_lat + else if( iadd1 <= NY_BATHY-1 .and. iel1 == NX_BATHY ) then ! interpolates for points on longitude border value = dble(ibathy_topo(iel1,iadd1))*(1-ratio_lon)*(1.-ratio_lat) & + dble(ibathy_topo(1,iadd1))*ratio_lon*(1.-ratio_lat) & + dble(ibathy_topo(1,iadd1+1))*ratio_lon*ratio_lat & + dble(ibathy_topo(iel1,iadd1+1))*(1.-ratio_lon)*ratio_lat + else ! for points on latitude boundaries value = dble(ibathy_topo(iel1,iadd1)) @@ -323,6 +329,3 @@ subroutine get_topo_bathy(xlat,xlon,value,ibathy_topo) end subroutine get_topo_bathy -! ------------------------------------------- - - diff --git a/src/shared/parallel.f90 b/src/shared/parallel.f90 index 85754c9d4..4aff699ce 100644 --- a/src/shared/parallel.f90 +++ b/src/shared/parallel.f90 @@ -25,7 +25,6 @@ ! !===================================================================== - !------------------------------------------------------------------------------------------------- ! ! MPI wrapper functions @@ -62,7 +61,6 @@ subroutine finalize_mpi() end subroutine finalize_mpi - ! !------------------------------------------------------------------------------------------------- ! @@ -81,8 +79,6 @@ subroutine abort_mpi() end subroutine abort_mpi - - ! !------------------------------------------------------------------------------------------------- ! @@ -122,7 +118,6 @@ subroutine synchronize_all_comm(comm) end subroutine synchronize_all_comm - ! !------------------------------------------------------------------------------------------------- ! @@ -158,6 +153,7 @@ subroutine test_request(request,flag_result_test) call MPI_TEST(request,flag_result_test,msg_status,ier) end subroutine test_request + ! !------------------------------------------------------------------------------------------------- ! @@ -350,6 +346,9 @@ subroutine max_allreduce_i(buffer,countval) integer,dimension(countval) :: send ! seems not to be supported on all kind of MPI implementations... + !! DK DK: yes, I confirm, using MPI_IN_PLACE is tricky + !! DK DK (see the answer at http://stackoverflow.com/questions/17741574/in-place-mpi-reduce-crashes-with-openmpi + !! DK DK for how to use it right) !call MPI_ALLREDUCE(MPI_IN_PLACE, buffer, countval, MPI_INTEGER, MPI_MAX, MPI_COMM_WORLD, ier) send(:) = buffer(:) @@ -359,7 +358,6 @@ subroutine max_allreduce_i(buffer,countval) end subroutine max_allreduce_i - ! !------------------------------------------------------------------------------------------------- ! @@ -414,6 +412,24 @@ subroutine sum_all_dp(sendbuf, recvbuf) end subroutine sum_all_dp +! +!------------------------------------------------------------------------------------------------- +! + + subroutine sum_all_3Darray_dp(sendbuf, recvbuf, nx,ny,nz) + + use mpi + + implicit none + + integer :: nx,ny,nz + double precision, dimension(nx,ny,nz) :: sendbuf, recvbuf + integer :: ier + + ! this works only if the arrays are contiguous in memory (which is always the case for static arrays, as used in the code) + call MPI_REDUCE(sendbuf,recvbuf,nx*ny*nz,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ier) + + end subroutine sum_all_3Darray_dp ! !------------------------------------------------------------------------------------------------- @@ -530,7 +546,6 @@ subroutine bcast_all_singler(buffer) end subroutine bcast_all_singler - ! !------------------------------------------------------------------------------------------------- ! @@ -568,7 +583,6 @@ subroutine bcast_all_singledp(buffer) end subroutine bcast_all_singledp - ! !------------------------------------------------------------------------------------------------- ! @@ -645,12 +659,10 @@ subroutine bcast_all_l(buffer, countval) end subroutine bcast_all_l - ! !------------------------------------------------------------------------------------------------- ! - subroutine recv_singlei(recvbuf, dest, recvtag) use mpi @@ -672,7 +684,6 @@ end subroutine recv_singlei !------------------------------------------------------------------------------------------------- ! - subroutine recv_singlel(recvbuf, dest, recvtag) use mpi @@ -694,7 +705,6 @@ end subroutine recv_singlel !------------------------------------------------------------------------------------------------- ! - subroutine recv_i(recvbuf, recvcount, dest, recvtag) use mpi @@ -717,7 +727,6 @@ end subroutine recv_i !------------------------------------------------------------------------------------------------- ! - subroutine recv_cr(recvbuf, recvcount, dest, recvtag) use constants @@ -743,7 +752,6 @@ end subroutine recv_cr !------------------------------------------------------------------------------------------------- ! - subroutine recv_dp(recvbuf, recvcount, dest, recvtag) use mpi @@ -762,7 +770,6 @@ subroutine recv_dp(recvbuf, recvcount, dest, recvtag) end subroutine recv_dp - ! !------------------------------------------------------------------------------------------------- ! @@ -862,7 +869,6 @@ subroutine send_dp(sendbuf, sendcount, dest, sendtag) end subroutine send_dp - ! !------------------------------------------------------------------------------------------------- ! @@ -1009,7 +1015,6 @@ subroutine gather_all_dp(sendbuf, sendcnt, recvbuf, recvcount, NPROC) end subroutine gather_all_dp - ! !------------------------------------------------------------------------------------------------- ! @@ -1090,8 +1095,6 @@ end subroutine gatherv_all_r !------------------------------------------------------------------------------------------------- ! - - subroutine world_size(sizeval) use mpi diff --git a/src/shared/rthetaphi_xyz.f90 b/src/shared/rthetaphi_xyz.f90 index e41927b82..1e2b946bf 100644 --- a/src/shared/rthetaphi_xyz.f90 +++ b/src/shared/rthetaphi_xyz.f90 @@ -33,7 +33,9 @@ subroutine xyz_2_rthetaphi(x,y,z,r,theta,phi) implicit none - real(kind=CUSTOM_REAL) x,y,z,r,theta,phi + real(kind=CUSTOM_REAL), intent(in) :: x,y,z + real(kind=CUSTOM_REAL), intent(out) :: r,theta,phi + double precision xmesh,ymesh,zmesh ! distinguish between single and double precision for reals @@ -81,7 +83,9 @@ subroutine xyz_2_rthetaphi_dble(x,y,z,r,theta,phi) implicit none - double precision x,y,z,r,theta,phi + double precision, intent(in) :: x,y,z + double precision, intent(out) :: r,theta,phi + double precision xmesh,ymesh,zmesh xmesh = x diff --git a/src/specfem3D/check_stability.f90 b/src/specfem3D/check_stability.f90 index e770fec0f..df024b11e 100644 --- a/src/specfem3D/check_stability.f90 +++ b/src/specfem3D/check_stability.f90 @@ -25,7 +25,6 @@ ! !===================================================================== - subroutine check_stability() ! computes the maximum of the norm of the displacement @@ -553,7 +552,7 @@ subroutine write_timestamp_file(Usolidnorm_all,Ufluidnorm_all,b_Usolidnorm_all,b write(outputname,"('/timestamp_backward_and_adjoint',i6.6)") it endif - ! file output + ! timestamp file output open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',action='write') write(IOUT,*) 'Time step # ',it