writing out forces instead of acceleration

tested in the toy mesh for Los Alamos
elifo · Sep 19, 2023 · 0b02c89 · 0b02c89
1 parent 8f7ceae
commit 0b02c89
Show file tree

Hide file tree

Showing 5 changed files with 21 additions and 4 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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