Merge remote-tracking branch 'origin/development' into keyword-view

PRIVATE

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

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha5-283-gdacd08f39
+v3.0.0-alpha5-297-g5ecfba1e5

python/damask/_colormap.py

@@ -6,6 +6,10 @@ from pathlib import Path
 from typing import Sequence, Union, TextIO
 
 import numpy as np
+try:
+    from numpy.typing import ArrayLike
+except ImportError:
+    ArrayLike = Union[np.ndarray,Sequence[float]] # type: ignore
 import scipy.interpolate as interp
 import matplotlib as mpl
 if os.name == 'posix' and 'DISPLAY' not in os.environ:
@@ -78,8 +82,8 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def from_range(low: Sequence[float],
+    def from_range(low: ArrayLike,
-                   high: Sequence[float],
+                   high: ArrayLike,
                    name: str = 'DAMASK colormap',
                    N: int = 256,
                    model: str = 'rgb') -> 'Colormap':
@@ -129,7 +133,7 @@ class Colormap(mpl.colors.ListedColormap):
         if model.lower() not in toMsh:
             raise ValueError(f'Invalid color model: {model}.')
 
-        low_high = np.vstack((low,high))
+        low_high = np.vstack((low,high)).astype(float)
         out_of_bounds = np.bool_(False)
 
         if   model.lower() == 'rgb':
@@ -142,7 +146,7 @@ class Colormap(mpl.colors.ListedColormap):
             out_of_bounds = np.any(low_high[:,0]<0)
 
         if out_of_bounds:
-            raise ValueError(f'{model.upper()} colors {low} | {high} are out of bounds.')
+            raise ValueError(f'{model.upper()} colors {low_high[0]} | {low_high[1]} are out of bounds.')
 
         low_,high_ = map(toMsh[model.lower()],low_high)
         msh = map(functools.partial(Colormap._interpolate_msh,low=low_,high=high_),np.linspace(0,1,N))
@@ -225,7 +229,7 @@ class Colormap(mpl.colors.ListedColormap):
 
     def shade(self,
               field: np.ndarray,
-              bounds: Sequence[float] = None,
+              bounds: ArrayLike = None,
               gap: float = None) -> Image:
         """
         Generate PIL image of 2D field using colormap.
@@ -235,7 +239,7 @@ class Colormap(mpl.colors.ListedColormap):
         field : numpy.array, shape (:,:)
             Data to be shaded.
         bounds : sequence of float, len (2), optional
-            Value range (low,high) spanned by colormap.
+            Value range (left,right) spanned by colormap.
         gap : field.dtype, optional
             Transparent value. NaN will always be rendered transparent.
 
@@ -248,17 +252,17 @@ class Colormap(mpl.colors.ListedColormap):
         mask = np.logical_not(np.isnan(field) if gap is None else \
                np.logical_or (np.isnan(field), field == gap))                                       # mask NaN (and gap if present)
 
-        lo,hi = (field[mask].min(),field[mask].max()) if bounds is None else \
+        l,r = (field[mask].min(),field[mask].max()) if bounds is None else \
-                (min(bounds[:2]),max(bounds[:2]))
+              np.array(bounds,float)[:2]
 
-        delta,avg = hi-lo,0.5*(hi+lo)
+        delta,avg = r-l,0.5*abs(r+l)
 
-        if delta * 1e8 <= avg:                                                                      # delta is similar to numerical noise
+        if abs(delta) * 1e8 <= avg:                                                                 # delta is similar to numerical noise
-            hi,lo = hi+0.5*avg,lo-0.5*avg                                                           # extend range to have actual data centered within
+            l,r = l-0.5*avg*np.sign(delta),r+0.5*avg*np.sign(delta),                                # extend range to have actual data centered within
 
         return Image.fromarray(
             (np.dstack((
-                        self.colors[(np.round(np.clip((field-lo)/(hi-lo),0.0,1.0)*(self.N-1))).astype(np.uint16),:3],
+                        self.colors[(np.round(np.clip((field-l)/delta,0.0,1.0)*(self.N-1))).astype(np.uint16),:3],
                         mask.astype(float)
                        )
                       )*255

src/math.f90

@@ -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