Merge branch 'master' of https://github.com/sstaehler/kerner

inparam_basic

@@ -29,6 +29,10 @@ MESH_FILE_ABAQUS     'Meshes/flat_triangles.inp'
 # or for absolute values of the model parameters like Vp (false)
 KERNEL_FOR_RELATIVE_PERTURBATIONS  true
 
+# For plotting reasons one may wish to skip the integration over cell-volume.
+# Resulting kernels bear the unit [s/m^3]
+INTEGRATE_OVER_VOLUME false
+
 # Prefix of output file names.
 # Kernel files are called $OUTPUT_FILE_kernel.xdmf
 # Wavefield movies are called $OUTPUT_FILE_wavefield.xdmf

mc_integration.f90

@@ -103,12 +103,13 @@ end subroutine check_montecarlo_integral
 !-------------------------------------------------------------------------------
 
 !-------------------------------------------------------------------------------
-subroutine initialize_montecarlo(this, nfuncs, volume, allowed_error, allowed_relerror)
+subroutine initialize_montecarlo(this, nfuncs, volume, allowed_error, allowed_relerror, int_over_volume)
     class(integrated_type), intent(inout) :: this 
     integer, intent(in)                   :: nfuncs
     real(kind=dp), intent(in)             :: volume
     real(kind=dp), intent(in)             :: allowed_error
     real(kind=dp), intent(in), optional   :: allowed_relerror
+    logical, intent(in), optional         :: int_over_volume
 
     if(allocated(this%fsum)) deallocate(this%fsum)
     if(allocated(this%f2sum)) deallocate(this%f2sum)
@@ -133,12 +134,21 @@ subroutine initialize_montecarlo(this, nfuncs, volume, allowed_error, allowed_re
     this%allowed_variance = allowed_error ** 2
     this%converged        = .false.
     this%nmodels          = 0
+
     if(present(allowed_relerror)) then
        this%allowed_relerror = allowed_relerror
     else
        this%allowed_relerror = 1d-100
     end if
 
+    ! for plotting one may wish to skip volume integration
+    if(present(int_over_volume)) then
+       if (.not.int_over_volume) then
+          this%volume = 1.d0
+       end if
+    end if
+
+
     this%isinitialized    = .true.
 
 end subroutine initialize_montecarlo

slave.f90

@@ -364,10 +364,11 @@ function slave_work(parameters, sem_data, inv_mesh, fft_data) result(slave_resul
         ! Initialize basis kernel Monte Carlo integrals for current element
         iclockold = tick()
         do ibasisfunc = 1, nbasisfuncs_per_elem
-           call int_kernel(ibasisfunc)%initialize_montecarlo(parameters%nkernel,             &
-                                                             volume,                         &
-                                                             parameters%allowed_error,       &
-                                                             parameters%allowed_relative_error) 
+           call int_kernel(ibasisfunc)%initialize_montecarlo(parameters%nkernel,               &
+                                                             volume,                           &
+                                                             parameters%allowed_error,         &
+                                                             parameters%allowed_relative_error,&
+                                                             parameters%int_over_volume) 
         end do
 
         do while (.not.allallconverged(int_kernel)) ! Beginning of Monte Carlo loop

type_parameter.f90

@@ -56,6 +56,7 @@ module type_parameter
         logical                              :: write_smgr           = .true.
         logical                              :: quasirandom          = .true.
         logical                              :: relative_kernel      = .true.
+        logical                              :: int_over_volume      = .true.
         contains
            procedure, pass                   :: read_parameters
            procedure, pass                   :: read_receiver
@@ -171,6 +172,9 @@ subroutine read_parameters(this, input_file_in)
         case('KERNEL_FOR_RELATIVE_PERTURBATIONS')
            read(keyvalue, *) this%relative_kernel
 
+        case('INTEGRATE_OVER_VOLUME')
+           read(keyvalue, *) this%int_over_volume
+
         case('DECONVOLVE_STF')
            read(keyvalue, *) this%deconv_stf
 
@@ -221,6 +225,7 @@ subroutine read_parameters(this, input_file_in)
   call pbroadcast_log(this%write_smgr, 0)
   call pbroadcast_log(this%quasirandom, 0)
   call pbroadcast_log(this%relative_kernel, 0)
+  call pbroadcast_log(this%int_over_volume, 0)
   call pbroadcast_char(this%fftw_plan, 0)
   call pbroadcast_char(this%whattodo, 0)
   call pbroadcast_char(this%int_type, 0)