diff --git a/PRIVATE b/PRIVATE index 2ad27552c..96c32ba42 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2ad27552c43316735b6ef425737fe3c8a5231598 +Subproject commit 96c32ba4237a51eaad92cd139e1a716ee5b32493 diff --git a/python/damask/VERSION b/python/damask/VERSION index 7ba8f4262..42ae7a5db 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-283-gdacd08f39 +v3.0.0-alpha5-297-g5ecfba1e5 diff --git a/python/damask/_colormap.py b/python/damask/_colormap.py index 5a13f440d..ba7617fbb 100644 --- a/python/damask/_colormap.py +++ b/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], - high: Sequence[float], + def from_range(low: ArrayLike, + 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 \ - (min(bounds[:2]),max(bounds[:2])) + l,r = (field[mask].min(),field[mask].max()) if bounds is None else \ + 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 - hi,lo = hi+0.5*avg,lo-0.5*avg # extend range to have actual data centered within + if abs(delta) * 1e8 <= avg: # delta is similar to numerical noise + 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 diff --git a/src/math.f90 b/src/math.f90 index 207a27e78..2d8f8fb3b 100644 --- a/src/math.f90 +++ b/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 - sigma !< standard deviation - real(pReal), intent(in), optional :: width !< cut off as multiples of standard deviation + real(pReal), intent(out) :: x + real(pReal), intent(in), optional :: mu, sigma - real(pReal), dimension(2) :: rnd ! random numbers - real(pReal) :: scatter, & ! normalized scatter around mean - width_ + real(pReal) :: sigma_, mu_ + real(pReal), dimension(2) :: rnd - if (abs(sigma) < tol_math_check) then - math_sampleGaussVar = mu + + if (present(mu)) then + mu_ = mu else - if (present(width)) then - width_ = width - else - width_ = 3.0_pReal ! use +-3*sigma as default scatter - endif + mu_ = 0.0_pReal + end if - do - call random_number(rnd) - scatter = width_ * (2.0_pReal * rnd(1) - 1.0_pReal) - if (rnd(2) <= exp(-0.5_pReal * scatter**2)) exit ! test if scattered value is drawn - enddo + if (present(sigma)) then + sigma_ = sigma + else + sigma_ = 1.0_pReal + end if - math_sampleGaussVar = scatter * sigma - endif + call random_number(rnd) + 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 diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index c762039d4..bf0dc3b19 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/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) - from = 1 + sum(ini%N_sl(1:f-1)) - upto = sum(ini%N_sl(1:f)) - do s = from,upto - noise = [math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u), & - math_sampleGaussVar(0.0_pReal, ini%sigma_rho_u)] - stt%rho_sgl_mob_edg_pos(s,e) = ini%rho_u_ed_pos_0(f) + noise(1) - stt%rho_sgl_mob_edg_neg(s,e) = ini%rho_u_ed_neg_0(f) + noise(1) - 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 + do f = 1,size(ini%N_sl,1) + from = 1 + sum(ini%N_sl(1:f-1)) + 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) + call math_normal(stt%rho_sgl_mob_edg_neg(from:upto,:),ini%rho_u_ed_neg_0(f),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) + call math_normal(stt%rho_sgl_mob_scr_neg(from:upto,:),ini%rho_u_sc_neg_0(f),ini%sigma_rho_u) + stt%rho_dip_edg(from:upto,:) = ini%rho_d_ed_0(f) + stt%rho_dip_scr(from:upto,:) = ini%rho_d_sc_0(f) end do end if