save forces instead of displacement

I do this to write out acceleration and forces for Los Alamos group simulations -- quick and dirty style.
elifo · Dec 27, 2023 · d032c9a · d032c9a
1 parent 0479958
commit d032c9a
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 10 deletions.
diff --git a/.gitignore b/.gitignore
@@ -42,3 +42,4 @@ doc/USER_MANUAL/Schedule
 .*swn
 *~
 .DS_Store
+*.pyc
diff --git a/src/gpu/kernels/compute_elastic_seismogram_kernel.cu b/src/gpu/kernels/compute_elastic_seismogram_kernel.cu
@@ -51,6 +51,7 @@ __global__ void compute_elastic_seismogram_kernel(int nrec_local,
                                                   realw_const_p rmassx,
                                                   realw_const_p rmassy,
                                                   realw_const_p rmassz,
+                                                  int seismotype,
                                                   int it){
 
   int irec_local = blockIdx.x + blockIdx.y*gridDim.x;
@@ -86,15 +87,22 @@ __global__ void compute_elastic_seismogram_kernel(int nrec_local,
         realw hlagrange = hxir * hetar * hgammar;
         int iglob = d_ibool[INDEX4_PADDED(NGLLX,NGLLX,NGLLX,I,J,K,ispec)] - 1;
 
-        //sh_dxd[tx] = hlagrange * displ[NDIM*iglob];
-        // sh_dyd[tx] = hlagrange * displ[NDIM*iglob+1];
-        // sh_dzd[tx] = hlagrange * displ[NDIM*iglob+2];
-        // Elif: write forces for Los Alamos, use only save_accel!!!
-        sh_dxd[tx] = hlagrange * displ[NDIM*iglob]/ rmassx[iglob];
-        sh_dyd[tx] = hlagrange * displ[NDIM*iglob+1]/rmassy[iglob];
-        sh_dzd[tx] = hlagrange * displ[NDIM*iglob+2]/ rmassz[iglob];
-
 
+        // default way of writing out displ/veloc/accel
+        if (seismotype != 1) {
+        sh_dxd[tx] = hlagrange * displ[NDIM*iglob];
+        sh_dyd[tx] = hlagrange * displ[NDIM*iglob+1];
+        sh_dzd[tx] = hlagrange * displ[NDIM*iglob+2];
+        }
+        // Elif: write forces for Los Alamos, use save_displacement=True!!!
+        if ( seismotype == 1) {
+           sh_dxd[tx] = hlagrange * displ[NDIM*iglob]/ rmassx[iglob];
+           sh_dyd[tx] = hlagrange * displ[NDIM*iglob+1]/rmassy[iglob];
+           sh_dzd[tx] = hlagrange * displ[NDIM*iglob+2]/ rmassz[iglob];
+           // Elif: test I/O
+           if (tx == 0) printf("ELIF: writing out forces instead of displacement!");
+           if (tx == 0) printf("%d \n", seismotype);
+        }
         //debug
         //if (tx == 0) printf("thread %d %d %d - %f %f %f\n",ispec,iglob,irec_local,hlagrange,displ[0 + 2*iglob],displ[1 + 2*iglob]);
       }

diff --git a/src/gpu/kernels/kernel_proto.cu.h b/src/gpu/kernels/kernel_proto.cu.h
@@ -610,6 +610,7 @@ __global__ void compute_elastic_seismogram_kernel(int nrec_local,
                                                   realw_const_p rmassx,
                                                   realw_const_p rmassy,
                                                   realw_const_p rmassz,
+                                                  int seismotype,
                                                   int it);
 
 

diff --git a/src/gpu/write_seismograms_cuda.cu b/src/gpu/write_seismograms_cuda.cu
@@ -93,12 +93,15 @@ void FC_FUNC_(compute_seismograms_cuda,
     if (seismotype == 1){
       // deplacement
       if (mp->simulation_type == 1 || mp->simulation_type == 2){
-        displ = mp->d_displ;
+        //displ = mp->d_displ;
+        // ELIF: write out forces instead of displacement (quick & dirty :S)
+        displ = mp->d_accel;
         potential = mp->d_potential_acoustic;
       }else{
         // kernel simulations
         // reconstructed forward wavefield stored in b_displ, b_veloc, b_accel
-        displ = mp->d_b_displ;
+        //displ = mp->d_b_displ;
+        displ = mp->d_b_accel;
         potential = mp->d_b_potential_acoustic;
       }
       d_seismo = mp->d_seismograms_d;
@@ -176,6 +179,7 @@ void FC_FUNC_(compute_seismograms_cuda,
                                                                                  mp->d_rmassx,
                                                                                  mp->d_rmassy,
                                                                                  mp->d_rmassz,
+                                                                                 seismotype,
                                                                                  seismo_current);
       }
 #endif