Merge branch 'better-Gauss' into 'development'

improved sampling from normal distribution

See merge request damask/DAMASK!483

PRIVATE

@@ -1 +1 @@
-Subproject commit 2ad27552c43316735b6ef425737fe3c8a5231598
+Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493

src/math.f90

@@ -895,7 +895,7 @@ pure function math_33toVoigt6_stress(sigma) result(sigma_tilde)
 
 
   sigma_tilde = [sigma(1,1), sigma(2,2), sigma(3,3), &
-                 sigma(3,2), sigma(3,1), sigma(1,2)] 
+                 sigma(3,2), sigma(3,1), sigma(1,2)]
 
 end function math_33toVoigt6_stress
 
@@ -910,7 +910,7 @@ pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde)
 
 
   epsilon_tilde = [          epsilon(1,1),           epsilon(2,2),           epsilon(3,3), &
-                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)] 
+                   2.0_pReal*epsilon(3,2), 2.0_pReal*epsilon(3,1), 2.0_pReal*epsilon(1,2)]
 
 end function math_33toVoigt6_strain
 
@@ -961,39 +961,36 @@ pure function math_3333toVoigt66(m3333)
 end function math_3333toVoigt66
 
 
-
 !--------------------------------------------------------------------------------------------------
-!> @brief draw a random sample from Gauss variable
+!> @brief Draw a sample from a normal distribution.
+!> @details https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
+!>          https://masuday.github.io/fortran_tutorial/random.html
 !--------------------------------------------------------------------------------------------------
-real(pReal) function math_sampleGaussVar(mu, sigma, width)
+impure elemental subroutine math_normal(x,mu,sigma)
 
-  real(pReal), intent(in) ::            mu, &                                                       !< mean
+  real(pReal), intent(out) :: x
-                                        sigma                                                       !< standard deviation
+  real(pReal), intent(in), optional :: mu, sigma
-  real(pReal), intent(in), optional ::  width                                                       !< cut off as multiples of standard deviation
 
-  real(pReal), dimension(2) ::          rnd                                                         ! random numbers
+  real(pReal) :: sigma_, mu_
-  real(pReal) ::                        scatter, &                                                  ! normalized scatter around mean
+  real(pReal), dimension(2) :: rnd
-                                        width_
 
-  if (abs(sigma) < tol_math_check) then
+
-    math_sampleGaussVar = mu
+  if (present(mu)) then
+    mu_ = mu
   else
-    if (present(width)) then
+    mu_ = 0.0_pReal
-      width_ = width
+  end if
-    else
-      width_ = 3.0_pReal                                                                            ! use +-3*sigma as default scatter
-    endif
 
-    do
+  if (present(sigma)) then
-      call random_number(rnd)
+    sigma_ = sigma
-      scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal)
+  else
-      if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit                                              ! test if scattered value is drawn
+    sigma_ = 1.0_pReal
-    enddo
+  end if
 
-    math_sampleGaussVar = scatter * sigma
+  call random_number(rnd)
-  endif
+  x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(2.0_pReal*PI*(1.0_pReal - rnd(2)))
 
-end function math_sampleGaussVar
+end subroutine math_normal
 
 
 !--------------------------------------------------------------------------------------------------
@@ -1434,6 +1431,26 @@ subroutine selfTest
   if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) &
     error stop 'math_LeviCivita'
 
+  normal_distribution: block
+    real(pReal), dimension(500000) :: r
+    real(pReal) :: mu, sigma
+
+    call random_number(mu)
+    call random_number(sigma)
+
+    sigma = 1.0_pReal + sigma*5.0_pReal
+    mu = (mu-0.5_pReal)*10_pReal
+
+    call math_normal(r,mu,sigma)
+
+    if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) &
+      error stop 'math_normal(mu)'
+
+    mu = sum(r)/real(size(r),pReal)
+    if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) &
+      error stop 'math_normal(sigma)'
+  end block normal_distribution
+
 end subroutine selfTest
 
 end module math

src/phase_mechanical_plastic_nonlocal.f90

@@ -1592,21 +1592,15 @@ subroutine stateInit(ini,phase,Nentries)
         stt%rhoSglMobile(s,e) = densityBinning
       end do
     else                                ! homogeneous distribution with noise
-      do e = 1, Nentries
+      do f = 1,size(ini%N_sl,1)
-        do f = 1,size(ini%N_sl,1)
+        from = 1 + sum(ini%N_sl(1:f-1))
-          from = 1 + sum(ini%N_sl(1:f-1))
+        upto = sum(ini%N_sl(1:f))
-          upto = sum(ini%N_sl(1:f))
+        call math_normal(stt%rho_sgl_mob_edg_pos(from:upto,:),ini%rho_u_ed_pos_0(f),ini%sigma_rho_u)
-          do s = from,upto
+        call math_normal(stt%rho_sgl_mob_edg_neg(from:upto,:),ini%rho_u_ed_neg_0(f),ini%sigma_rho_u)
-            noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), &
+        call math_normal(stt%rho_sgl_mob_scr_pos(from:upto,:),ini%rho_u_sc_pos_0(f),ini%sigma_rho_u)
-                     math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)]
+        call math_normal(stt%rho_sgl_mob_scr_neg(from:upto,:),ini%rho_u_sc_neg_0(f),ini%sigma_rho_u)
-            stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1)
+        stt%rho_dip_edg(from:upto,:) = ini%rho_d_ed_0(f)
-            stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1)
+        stt%rho_dip_scr(from:upto,:) = ini%rho_d_sc_0(f)
-            stt%rho_sgl_mob_scr_pos(s,e) = ini%rho_u_sc_pos_0(f) + noise(2)
-            stt%rho_sgl_mob_scr_neg(s,e) = ini%rho_u_sc_neg_0(f) + noise(2)
-          end do
-          stt%rho_dip_edg(from:upto,e)   = ini%rho_d_ed_0(f)
-          stt%rho_dip_scr(from:upto,e)   = ini%rho_d_sc_0(f)
-        end do
       end do
     end if