From 0b02c89ccba9fe5db990eb2212fe8c5a79967d45 Mon Sep 17 00:00:00 2001 From: Elif Oral Date: Tue, 19 Sep 2023 13:44:42 -0700 Subject: [PATCH] writing out forces instead of acceleration tested in the toy mesh for Los Alamos --- .../compute_forces_viscoelastic_calling_routine.F90 | 3 +++ src/specfem3D/compute_seismograms.f90 | 10 ++++++++-- src/specfem3D/read_mesh_databases.F90 | 4 ++++ src/specfem3D/setup_sources_receivers.f90 | 6 +++++- src/specfem3D/specfem3D_par.F90 | 2 +- 5 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 index 9cd1e8a79..bb7b03cb8 100644 --- a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 +++ b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 @@ -360,6 +360,9 @@ subroutine compute_forces_viscoelastic_calling() if (.not. GPU_MODE) then ! on CPU do iglob = 1,NGLOB_AB + ! Elif: assign before overwriting to write out for Los Alamos + force_losalamos(:,iglob) = accel(:,iglob) + accel(1,iglob) = accel(1,iglob)*rmassx(iglob) accel(2,iglob) = accel(2,iglob)*rmassy(iglob) accel(3,iglob) = accel(3,iglob)*rmassz(iglob) diff --git a/src/specfem3D/compute_seismograms.f90 b/src/specfem3D/compute_seismograms.f90 index 279187a6c..62abbf399 100644 --- a/src/specfem3D/compute_seismograms.f90 +++ b/src/specfem3D/compute_seismograms.f90 @@ -45,7 +45,7 @@ subroutine compute_seismograms() ! wavefields use specfem_par_acoustic, only: ispec_is_acoustic,potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, & b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic - use specfem_par_elastic, only: ispec_is_elastic,displ,veloc,accel, & + use specfem_par_elastic, only: ispec_is_elastic,displ,veloc,accel,force_losalamos, & b_displ,b_veloc,b_accel use specfem_par_poroelastic, only: ispec_is_poroelastic,displs_poroelastic,velocs_poroelastic,accels_poroelastic, & b_displs_poroelastic,b_velocs_poroelastic,b_accels_poroelastic @@ -118,10 +118,16 @@ subroutine compute_seismograms() ! elastic wave field if (ispec_is_elastic(ispec)) then ! interpolates displ/veloc/accel at receiver locations - call compute_interpolated_dva_viscoelast(displ,veloc,accel,NGLOB_AB, & + !call compute_interpolated_dva_viscoelast(displ,veloc,accel,NGLOB_AB, & + ! ispec,NSPEC_AB,ibool, & + ! hxir,hetar,hgammar, & + ! dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd,pd) + ! Elif: write out forces instead of accel + call compute_interpolated_dva_viscoelast(displ,veloc,force_losalamos,NGLOB_AB, & ispec,NSPEC_AB,ibool, & hxir,hetar,hgammar, & dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd,pd) + endif ! elastic ! acoustic wave field diff --git a/src/specfem3D/read_mesh_databases.F90 b/src/specfem3D/read_mesh_databases.F90 index 301885aef..f0deb63a3 100644 --- a/src/specfem3D/read_mesh_databases.F90 +++ b/src/specfem3D/read_mesh_databases.F90 @@ -268,7 +268,11 @@ subroutine read_mesh_databases() allocate(accel(NDIM,NGLOB_AB),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 1429') if (ier /= 0) stop 'Error allocating array accel' + allocate(force_losalamos(NDIM,NGLOB_AB),stat=ier) + if (ier /= 0) call exit_MPI_without_rank('error allocating array 1430') + if (ier /= 0) stop 'Error allocating array force_losalamos' displ(:,:) = 0.0_CUSTOM_REAL; veloc(:,:) = 0.0_CUSTOM_REAL; accel(:,:) = 0.0_CUSTOM_REAL + force_losalamos(:,:) = 0.0_CUSTOM_REAL if (SIMULATION_TYPE /= 1) then allocate(accel_adj_coupling(NDIM,NGLOB_AB),stat=ier) diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90 index ce11d8112..70c0474c5 100644 --- a/src/specfem3D/setup_sources_receivers.f90 +++ b/src/specfem3D/setup_sources_receivers.f90 @@ -43,15 +43,19 @@ subroutine setup_sources_receivers() ! reads in stations file and locates receivers call setup_receivers() + write(*,*) 'ELIF:: DEBUG:: done: setup_receivers' ! pre-compute source arrays call setup_sources_precompute_arrays() + write(*,*) 'ELIF:: DEBUG:: done: setup_sources_precompute_arrays' ! pre-compute receiver interpolation factors call setup_receivers_precompute_intp() + write(*,*) 'ELIF:: DEBUG:: done: setup_receivers_precompute_intp' ! write source and receiver VTK files for Paraview - call setup_sources_receivers_VTKfile() + write(*,*) 'ELIF:: DEBUG:: WRITE_SAVE_SOURCES_RECEIVERS_VTK', WRITE_SAVE_SOURCES_RECEIVERS_VTK + !if (WRITE_SAVE_SOURCES_RECEIVERS_VTK) & call setup_sources_receivers_VTKfile() ! user output if (myrank == 0) then diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90 index f8b65264b..d7d9ef3a1 100644 --- a/src/specfem3D/specfem3D_par.F90 +++ b/src/specfem3D/specfem3D_par.F90 @@ -396,7 +396,7 @@ module specfem_par_elastic real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: epsilon_trace_over_3 ! displacement, velocity, acceleration - real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: displ,veloc,accel + real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: displ,veloc,accel,force_losalamos real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_adj_coupling ! mass matrix