From 4a84c42112dab832fedd0f354747f7f7aaa5d8dc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 13:49:41 +0100 Subject: [PATCH 01/26] improving constants IUPAC does not specify an official ambient temperature, T=273.15K from standard temperature and pressure (STP) is not particularly suited for our purposes --- src/constants.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 3b99763de..43098188b 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -9,8 +9,8 @@ module constants public real(pReal), parameter :: & - T_ROOM = 300.0_pReal, & !< Room temperature in K. ToDo: IUPAC: 298.15 - K_B = 1.38e-23_pReal, & !< Boltzmann constant in J/Kelvin - N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol + T_ROOM = 298.15_pReal, & !< Room temperature in K (25°C)/Standard Ambient Temperaure and Pressure (SATP) + K_B = 1.380649e-23_pReal, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook) + N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook) end module constants From dd0f2cfa3c96bd4d8b94930b4940af6dd6093a2d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 14:39:32 +0100 Subject: [PATCH 02/26] tau is often easier to use than pi https://tauday.com/tau-manifesto --- src/grid/spectral_utilities.f90 | 40 +++++++++---------- src/math.f90 | 18 +++++---- ...phase_mechanical_plastic_dislotungsten.f90 | 2 +- src/rotations.f90 | 16 ++++---- 4 files changed, 39 insertions(+), 37 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index ee5fd4c82..25de88e3f 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -936,36 +936,36 @@ pure function utilities_getFreqDerivative(k_s) select case (spectral_derivative_ID) case (DERIVATIVE_CONTINUOUS_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, 2.0_pReal*PI*real(k_s,pReal)/geomSize,pReal) + utilities_getFreqDerivative = cmplx(0.0_pReal, TAU*real(k_s,pReal)/geomSize,pReal) case (DERIVATIVE_CENTRAL_DIFF_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, sin(2.0_pReal*PI*real(k_s,pReal)/real(grid,pReal)), pReal)/ & + utilities_getFreqDerivative = cmplx(0.0_pReal, sin(TAU*real(k_s,pReal)/real(grid,pReal)), pReal)/ & cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) case (DERIVATIVE_FWBW_DIFF_ID) utilities_getFreqDerivative(1) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & + cmplx(cos(TAU*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(2) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & + cmplx(cos(TAU*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(3) = & - cmplx(cos(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & - sin(2.0_pReal*PI*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & + cmplx(cos(TAU*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) end select diff --git a/src/math.f90 b/src/math.f90 index db0666ab0..1eecd7497 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -21,10 +21,12 @@ module math config #endif - real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree - real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian - complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) + real(pReal), parameter :: & + PI = acos(-1.0_pReal), & !< ratio of a circle's circumference to its diameter + TAU = 2.0_pReal*PI, & !< ratio of a circle's circumference to its radius + INDEG = 360.0_pReal/TAU, & !< conversion from radian to degree + INRAD = TAU/360.0_pReal !< conversion from degree to radian + complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,TAU) !< Re(0.0), Im(Tau) real(pReal), dimension(3,3), parameter :: & math_I3 = reshape([& @@ -984,7 +986,7 @@ impure elemental subroutine math_normal(x,mu,sigma) end if 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))) + x = mu_ + sigma_ * sqrt(-2.0_pReal*log(1.0_pReal-rnd(1)))*cos(TAU*(1.0_pReal - rnd(2))) end subroutine math_normal @@ -1088,7 +1090,7 @@ pure function math_rotationalPart(F) result(R) if (dNeq0(x)) then Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal)) lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) & - *cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) + *cos((Phi-TAU*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal) else lambda = sqrt(I_C(1)/3.0_pReal) @@ -1154,8 +1156,8 @@ pure function math_eigvalsh33(m) phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) math_eigvalsh33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & [cos( phi /3.0_pReal), & - cos((phi+2.0_pReal*PI)/3.0_pReal), & - cos((phi+4.0_pReal*PI)/3.0_pReal) & + cos((phi+TAU)/3.0_pReal), & + cos((phi+2.0_pReal*TAU)/3.0_pReal) & ] & + I(1)/3.0_pReal endif diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index cd71a7fd7..9f0465158 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -341,7 +341,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & 0.0_pReal, & prm%dipoleformation) - v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) & + v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(TAU*K_B*T)) & * (1.0_pReal/(d_hat+prm%d_caron)) dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? end where diff --git a/src/rotations.f90 b/src/rotations.f90 index b73fbe8da..7c96d87e1 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -198,7 +198,7 @@ subroutine fromEulers(self,eu,degrees) Eulers = merge(eu*INRAD,eu,degrees) endif - if (any(Eulers<0.0_pReal) .or. any(Eulers>2.0_pReal*PI) .or. Eulers(2) > PI) & + if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) & call IO_error(402,ext_msg='fromEulers') self%q = eu2qu(Eulers) @@ -480,7 +480,7 @@ pure function qu2eu(qu) result(eu) atan2( 2.0_pReal*chi, q03-q12 ), & atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] endif degenerated - where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU]) end function qu2eu @@ -628,7 +628,7 @@ pure function om2eu(om) result(eu) eu = [atan2(om(1,2),om(1,1)), 0.5_pReal*PI*(1.0_pReal-om(3,3)),0.0_pReal ] end if where(abs(eu) < 1.e-8_pReal) eu = 0.0_pReal - where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU]) end function om2eu @@ -1437,7 +1437,7 @@ subroutine selfTest() elseif(i==2) then qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal]) elseif(i==3) then - qu = eu2qu([2.0_pReal*PI,PI,2.0_pReal*PI]) + qu = eu2qu([TAU,PI,TAU]) elseif(i==4) then qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal] elseif(i==5) then @@ -1448,10 +1448,10 @@ subroutine selfTest() call random_number(x) A = sqrt(x(3)) B = sqrt(1-0_pReal -x(3)) - qu = [cos(2.0_pReal*PI*x(1))*A,& - sin(2.0_pReal*PI*x(2))*B,& - cos(2.0_pReal*PI*x(2))*B,& - sin(2.0_pReal*PI*x(1))*A] + qu = [cos(TAU*x(1))*A,& + sin(TAU*x(2))*B,& + cos(TAU*x(2))*B,& + sin(TAU*x(1))*A] if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal) endif From 774ec4e8c4e37124c1e7941c1e8eb6638f06b1e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 14:47:14 +0100 Subject: [PATCH 03/26] not used --- src/math.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/math.f90 b/src/math.f90 index 1eecd7497..0343d190a 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -26,7 +26,6 @@ module math TAU = 2.0_pReal*PI, & !< ratio of a circle's circumference to its radius INDEG = 360.0_pReal/TAU, & !< conversion from radian to degree INRAD = TAU/360.0_pReal !< conversion from degree to radian - complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,TAU) !< Re(0.0), Im(Tau) real(pReal), dimension(3,3), parameter :: & math_I3 = reshape([& From 477f2d81678d416d1025814156bb880c5da6e699 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 14:49:45 +0100 Subject: [PATCH 04/26] capitalize constants --- src/rotations.f90 | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 7c96d87e1..ec000af60 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -81,7 +81,6 @@ module rotations end type rotation real(pReal), parameter :: & - SPI = sqrt(PI), & PREF = sqrt(6.0_pReal/PI), & A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & AP = PI**(2.0_pReal/3.0_pReal), & @@ -1209,7 +1208,7 @@ pure function ho2cu(ho) result(cu) else special q2 = qxy + maxval(abs(xyz2))**2 sq2 = sqrt(q2) - q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) + q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & @@ -1217,7 +1216,7 @@ pure function ho2cu(ho) result(cu) endif special ! inverse M_1 - xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc + xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC ! reverse the coordinates back to order according to the original pyramid number cu = xyz1(p(:,2)) @@ -1323,26 +1322,26 @@ pure function cu2ho(cu) result(ho) else center ! get pyramide and scale by grid parameter ratio p = GetPyramidOrder(cu) - XYZ = cu(p(:,1)) * sc + XYZ = cu(p(:,1)) * SC ! intercept all the points along the z-axis special: if (all(dEq0(XYZ(1:2)))) then - LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ] + LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ] else special order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger c = cos(q) s = sin(q) - q = prek * XYZ(order(2))/ sqrt(R2-c) + q = PREK * XYZ(order(2))/ sqrt(R2-c) T = [ (R2*c - 1.0), R2 * s] * q ! transform to sphere grid (inverse Lambert) ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] c = sum(T**2) s = Pi * c/(24.0*XYZ(3)**2) - c = sPi * c / sqrt(24.0_pReal) / XYZ(3) + c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3) q = sqrt( 1.0 - s ) - LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ] + LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ] endif special ! reverse the coordinates back to order according to the original pyramid number From 487912cfb0be72b18a4797fc81af750960c33eb1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 15:05:07 +0100 Subject: [PATCH 05/26] following Python notation --- src/grid/discretization_grid.f90 | 140 ++++++------- src/grid/grid_damage_spectral.f90 | 36 ++-- src/grid/grid_mech_FEM.f90 | 54 ++--- src/grid/grid_mech_spectral_basic.f90 | 32 +-- src/grid/grid_mech_spectral_polarisation.f90 | 50 ++--- src/grid/grid_thermal_spectral.f90 | 30 +-- src/grid/spectral_utilities.f90 | 200 +++++++++---------- 7 files changed, 271 insertions(+), 271 deletions(-) diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 8de2de6ef..67f735a0d 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -27,10 +27,10 @@ module discretization_grid private integer, dimension(3), public, protected :: & - grid !< (global) grid + cells !< (global) cells integer, public, protected :: & - grid3, & !< (local) grid in 3rd direction - grid3Offset !< (local) grid offset in 3rd direction + cells3, & !< (local) cells in 3rd direction + cells3Offset !< (local) cells offset in 3rd direction real(pReal), dimension(3), public, protected :: & geomSize !< (global) physical size real(pReal), public, protected :: & @@ -55,7 +55,7 @@ subroutine discretization_grid_init(restart) mySize, & !< domain size of this process origin !< (global) distance to origin integer, dimension(3) :: & - myGrid !< domain grid of this process + myGrid !< domain cells of this process integer, dimension(:), allocatable :: & materialAt, materialAt_global @@ -77,7 +77,7 @@ subroutine discretization_grid_init(restart) if (worldrank == 0) then fileContent = IO_read(interface_geomFile) - call readVTI(grid,geomSize,origin,materialAt_global,fileContent) + call readVTI(cells,geomSize,origin,materialAt_global,fileContent) fname = interface_geomFile if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.) @@ -88,37 +88,37 @@ subroutine discretization_grid_init(restart) end if - call MPI_Bcast(grid,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) + call MPI_Bcast(cells,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1') + if (cells(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1') call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - print'(/,1x,a,3(i12,1x))', 'cells a b c: ', grid + print'(/,1x,a,3(i12,1x))', 'cells a b c: ', cells print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin - if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)') + if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)') call fftw_mpi_init - devNull = fftw_mpi_local_size_3d(int(grid(3),C_INTPTR_T), & - int(grid(2),C_INTPTR_T), & - int(grid(1),C_INTPTR_T)/2+1, & + devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T), & + int(cells(2),C_INTPTR_T), & + int(cells(1),C_INTPTR_T)/2+1, & PETSC_COMM_WORLD, & - z, & ! domain grid size along z - z_offset) ! domain grid offset along z + z, & ! domain cells size along z + z_offset) ! domain cells offset along z if (z==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes') - grid3 = int(z) - grid3Offset = int(z_offset) - size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal) - size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal) - myGrid = [grid(1:2),grid3] + cells3 = int(z) + cells3Offset = int(z_offset) + size3 = geomSize(3)*real(cells3,pReal) /real(cells(3),pReal) + size3Offset = geomSize(3)*real(cells3Offset,pReal)/real(cells(3),pReal) + myGrid = [cells(1:2),cells3] mySize = [geomSize(1:2),size3] - call MPI_Gather(product(grid(1:2))*grid3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,& + call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,& 1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_Gather(product(myGrid), 1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,& @@ -131,10 +131,10 @@ subroutine discretization_grid_init(restart) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call discretization_init(materialAt, & - IPcoordinates0(myGrid,mySize,grid3Offset), & - Nodes0(myGrid,mySize,grid3Offset),& - merge((grid(1)+1) * (grid(2)+1) * (grid3+1),& ! write top layer... - (grid(1)+1) * (grid(2)+1) * grid3,& ! ...unless not last process + IPcoordinates0(myGrid,mySize,cells3Offset), & + Nodes0(myGrid,mySize,cells3Offset),& + merge((cells(1)+1) * (cells(2)+1) * (cells3+1),& ! write top layer... + (cells(1)+1) * (cells(2)+1) * cells3,& ! ...unless not last process worldrank+1==worldsize)) !-------------------------------------------------------------------------------------------------- @@ -142,7 +142,7 @@ subroutine discretization_grid_init(restart) if (.not. restart) then call results_openJobFile call results_closeGroup(results_addGroup('geometry')) - call results_addAttribute('cells', grid, '/geometry') + call results_addAttribute('cells', cells, '/geometry') call results_addAttribute('size', geomSize,'/geometry') call results_addAttribute('origin',origin, '/geometry') call results_closeJobFile @@ -170,13 +170,13 @@ end subroutine discretization_grid_init !> @brief Parse vtk image data (.vti) !> @details https://vtk.org/Wiki/VTK_XML_Formats !-------------------------------------------------------------------------------------------------- -subroutine readVTI(grid,geomSize,origin,material, & +subroutine readVTI(cells,geomSize,origin,material, & fileContent) integer, dimension(3), intent(out) :: & - grid ! grid (across all processes!) + cells ! cells (across all processes!) real(pReal), dimension(3), intent(out) :: & - geomSize, & ! size (across all processes!) + geomSize, & ! size (across all processes!) origin ! origin (across all processes!) integer, dimension(:), intent(out), allocatable :: & material @@ -190,7 +190,7 @@ subroutine readVTI(grid,geomSize,origin,material, & s - grid = -1 + cells = -1 geomSize = -1.0_pReal inFile = .false. @@ -215,7 +215,7 @@ subroutine readVTI(grid,geomSize,origin,material, & if (.not. inImage) then if (index(fileContent(startPos:endPos),' @brief Calculate undeformed position of IPs/cell centers (pretend to be an element) !--------------------------------------------------------------------------------------------------- -function IPcoordinates0(grid,geomSize,grid3Offset) +function IPcoordinates0(cells,geomSize,cells3Offset) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, intent(in) :: grid3Offset ! grid(3) offset + integer, intent(in) :: cells3Offset ! cells(3) offset - real(pReal), dimension(3,product(grid)) :: ipCoordinates0 + real(pReal), dimension(3,product(cells)) :: ipCoordinates0 integer :: & a,b,c, & @@ -516,9 +516,9 @@ function IPcoordinates0(grid,geomSize,grid3Offset) i = 0 - do c = 1, grid(3); do b = 1, grid(2); do a = 1, grid(1) + do c = 1, cells(3); do b = 1, cells(2); do a = 1, cells(1) i = i + 1 - IPcoordinates0(1:3,i) = geomSize/real(grid,pReal) * (real([a,b,grid3Offset+c],pReal) -0.5_pReal) + IPcoordinates0(1:3,i) = geomSize/real(cells,pReal) * (real([a,b,cells3Offset+c],pReal) -0.5_pReal) end do; end do; end do end function IPcoordinates0 @@ -527,22 +527,22 @@ end function IPcoordinates0 !--------------------------------------------------------------------------------------------------- !> @brief Calculate position of undeformed nodes (pretend to be an element) !--------------------------------------------------------------------------------------------------- -pure function nodes0(grid,geomSize,grid3Offset) +pure function nodes0(cells,geomSize,cells3Offset) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, intent(in) :: grid3Offset ! grid(3) offset + integer, intent(in) :: cells3Offset ! cells(3) offset - real(pReal), dimension(3,product(grid+1)) :: nodes0 + real(pReal), dimension(3,product(cells+1)) :: nodes0 integer :: & a,b,c, & n n = 0 - do c = 0, grid3; do b = 0, grid(2); do a = 0, grid(1) + do c = 0, cells3; do b = 0, cells(2); do a = 0, cells(1) n = n + 1 - nodes0(1:3,n) = geomSize/real(grid,pReal) * real([a,b,grid3Offset+c],pReal) + nodes0(1:3,n) = geomSize/real(cells,pReal) * real([a,b,cells3Offset+c],pReal) end do; end do; end do end function nodes0 @@ -551,17 +551,17 @@ end function nodes0 !-------------------------------------------------------------------------------------------------- !> @brief Calculate IP interface areas !-------------------------------------------------------------------------------------------------- -pure function cellSurfaceArea(geomSize,grid) +pure function cellSurfaceArea(geomSize,cells) real(pReal), dimension(3), intent(in) :: geomSize ! size (for this process!) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) - real(pReal), dimension(6,1,product(grid)) :: cellSurfaceArea + real(pReal), dimension(6,1,product(cells)) :: cellSurfaceArea - cellSurfaceArea(1:2,1,:) = geomSize(2)/real(grid(2)) * geomSize(3)/real(grid(3)) - cellSurfaceArea(3:4,1,:) = geomSize(3)/real(grid(3)) * geomSize(1)/real(grid(1)) - cellSurfaceArea(5:6,1,:) = geomSize(1)/real(grid(1)) * geomSize(2)/real(grid(2)) + cellSurfaceArea(1:2,1,:) = geomSize(2)/real(cells(2)) * geomSize(3)/real(cells(3)) + cellSurfaceArea(3:4,1,:) = geomSize(3)/real(cells(3)) * geomSize(1)/real(cells(1)) + cellSurfaceArea(5:6,1,:) = geomSize(1)/real(cells(1)) * geomSize(2)/real(cells(2)) end function cellSurfaceArea @@ -588,42 +588,42 @@ end function cellSurfaceNormal !-------------------------------------------------------------------------------------------------- !> @brief Build IP neighborhood relations !-------------------------------------------------------------------------------------------------- -pure function IPneighborhood(grid) +pure function IPneighborhood(cells) - integer, dimension(3), intent(in) :: grid ! grid (for this process!) + integer, dimension(3), intent(in) :: cells ! cells (for this process!) - integer, dimension(3,6,1,product(grid)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID] + integer, dimension(3,6,1,product(cells)) :: IPneighborhood !< 6 neighboring IPs as [element ID, IP ID, face ID] integer :: & x,y,z, & e e = 0 - do z = 0,grid(3)-1; do y = 0,grid(2)-1; do x = 0,grid(1)-1 + do z = 0,cells(3)-1; do y = 0,cells(2)-1; do x = 0,cells(1)-1 e = e + 1 ! element ID - IPneighborhood(1,1,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x+1,grid(1)) & + IPneighborhood(1,1,1,e) = z * cells(1) * cells(2) & + + y * cells(1) & + + modulo(x+1,cells(1)) & + 1 - IPneighborhood(1,2,1,e) = z * grid(1) * grid(2) & - + y * grid(1) & - + modulo(x-1,grid(1)) & + IPneighborhood(1,2,1,e) = z * cells(1) * cells(2) & + + y * cells(1) & + + modulo(x-1,cells(1)) & + 1 - IPneighborhood(1,3,1,e) = z * grid(1) * grid(2) & - + modulo(y+1,grid(2)) * grid(1) & + IPneighborhood(1,3,1,e) = z * cells(1) * cells(2) & + + modulo(y+1,cells(2)) * cells(1) & + x & + 1 - IPneighborhood(1,4,1,e) = z * grid(1) * grid(2) & - + modulo(y-1,grid(2)) * grid(1) & + IPneighborhood(1,4,1,e) = z * cells(1) * cells(2) & + + modulo(y-1,cells(2)) * cells(1) & + x & + 1 - IPneighborhood(1,5,1,e) = modulo(z+1,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & + IPneighborhood(1,5,1,e) = modulo(z+1,cells(3)) * cells(1) * cells(2) & + + y * cells(1) & + x & + 1 - IPneighborhood(1,6,1,e) = modulo(z-1,grid(3)) * grid(1) * grid(2) & - + y * grid(1) & + IPneighborhood(1,6,1,e) = modulo(z-1,cells(3)) * cells(1) * cells(2) & + + y * cells(1) & + x & + 1 ! IP ID diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 96b72dbcc..dc21a1a2c 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -106,9 +106,9 @@ subroutine grid_damage_spectral_init() !-------------------------------------------------------------------------------------------------- ! init fields - allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal) - allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) - allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) + allocate(phi_current(cells(1),cells(2),cells3), source=1.0_pReal) + allocate(phi_lastInc(cells(1),cells(2),cells3), source=1.0_pReal) + allocate(phi_stagInc(cells(1),cells(2),cells3), source=1.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -117,23 +117,23 @@ subroutine grid_damage_spectral_init() call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 1_pPetscInt, 0_pPetscInt, & ! #dof (phi, scalar), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells damage_grid,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(damage_grid,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(damage_grid,err_PETSc) CHKERRQ(err_PETSc) - call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -213,7 +213,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution) !-------------------------------------------------------------------------------------------------- ! updating damage state ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_set_phi(phi_current(i,j,k),ce) end do; end do; end do @@ -255,7 +255,7 @@ subroutine grid_damage_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_set_phi(phi_current(i,j,k),ce) end do; end do; end do @@ -289,12 +289,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current + scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_FFTvectorBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) end do; end do; end do @@ -302,7 +302,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field call utilities_FFTscalarBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & @@ -315,14 +315,14 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) call utilities_FFTscalarBackward - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness + where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) & + scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc + where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) & + scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness !-------------------------------------------------------------------------------------------------- ! constructing residual - r = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current + r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current err_PETSc = 0 end subroutine formResidual @@ -339,7 +339,7 @@ subroutine updateReference() K_ref = 0.0_pReal mu_ref = 0.0_pReal - do ce = 1, product(grid(1:2))*grid3 + do ce = 1, product(cells(1:2))*cells3 K_ref = K_ref + homogenization_K_phi(ce) mu_ref = mu_ref + homogenization_mu_phi(ce) end do diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 6a3f1323d..cc082d2f5 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -153,9 +153,9 @@ subroutine grid_mechanical_FEM_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(P_current (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -164,16 +164,16 @@ subroutine grid_mechanical_FEM_init call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if(err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, & DMDA_STENCIL_BOX, & - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 3_pPetscInt, 1_pPetscInt, & ! #dof (u, vector), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells mechanical_grid,err_PETSc) CHKERRQ(err_PETSc) call DMsetFromOptions(mechanical_grid,err_PETSc) @@ -214,7 +214,7 @@ subroutine grid_mechanical_FEM_init call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc) CHKERRQ(err_PETSc) - delta = geomSize/real(grid,pReal) ! grid spacing + delta = geomSize/real(cells,pReal) ! grid spacing detJ = product(delta) ! cell volume BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), & @@ -255,11 +255,11 @@ subroutine grid_mechanical_FEM_init call HDF5_read(u_lastInc,groupHandle,'u_lastInc') elseif (interface_restartInc == 0) then restartRead - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) + F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity + F = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response + homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 F, & ! target F @@ -386,7 +386,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai F_lastInc = F - homogenization_F0 = reshape(F, [3,3,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3]) endif !-------------------------------------------------------------------------------------------------- @@ -556,13 +556,13 @@ subroutine formResidual(da_local,x_local, & ! get deformation gradient call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo - F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) + F(1:3,1:3,i,j,k-cells3Offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) enddo; enddo; enddo call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) @@ -589,14 +589,14 @@ subroutine formResidual(da_local,x_local, & call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc) CHKERRQ(err_PETSc) ele = 0 - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk) enddo; enddo; enddo ele = ele + 1 - f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-grid3offset)))*detJ + & + f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-cells3Offset)))*detJ + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + & homogenization_dPdF(2,2,2,2,ele) + & homogenization_dPdF(3,3,3,3,ele))/3.0_pReal @@ -615,17 +615,17 @@ subroutine formResidual(da_local,x_local, & ! applying boundary conditions call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc) CHKERRQ(err_PETSc) - if (grid3offset == 0) then - r(0:2,0, 0, 0) = 0.0_pReal - r(0:2,grid(1),0, 0) = 0.0_pReal - r(0:2,0, grid(2),0) = 0.0_pReal - r(0:2,grid(1),grid(2),0) = 0.0_pReal + if (cells3Offset == 0) then + r(0:2,0, 0, 0) = 0.0_pReal + r(0:2,cells(1),0, 0) = 0.0_pReal + r(0:2,0, cells(2),0) = 0.0_pReal + r(0:2,cells(1),cells(2),0) = 0.0_pReal end if - if (grid3+grid3offset == grid(3)) then - r(0:2,0, 0, grid(3)) = 0.0_pReal - r(0:2,grid(1),0, grid(3)) = 0.0_pReal - r(0:2,0, grid(2),grid(3)) = 0.0_pReal - r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal + if (cells3+cells3Offset == cells(3)) then + r(0:2,0, 0, cells(3)) = 0.0_pReal + r(0:2,cells(1),0, cells(3)) = 0.0_pReal + r(0:2,0, cells(2),cells(3)) = 0.0_pReal + r(0:2,cells(1),cells(2),cells(3)) = 0.0_pReal end if call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc) CHKERRQ(err_PETSc) @@ -663,7 +663,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc) call MatZeroEntries(Jac,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ctr = 0 do kk = -1, 0; do jj = -1, 0; do ii = -1, 0 ctr = ctr + 1 @@ -719,7 +719,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc) call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, cells(1) ce = ce + 1 x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce) enddo; enddo; enddo diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index c60c542d0..2afe9edc7 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -96,7 +96,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_mechanical_spectral_basic_init - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -153,8 +153,8 @@ subroutine grid_mechanical_spectral_basic_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate(F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -163,23 +163,23 @@ subroutine grid_mechanical_spectral_basic_init call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 9_pPetscInt, 0_pPetscInt, & ! #dof (F, tensor), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells da,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(da,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(da,err_PETSc) CHKERRQ(err_PETSc) - call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 9, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -217,11 +217,11 @@ subroutine grid_mechanical_spectral_basic_init call HDF5_read(F_lastInc,groupHandle,'F_lastInc') elseif (interface_restartInc == 0) then restartRead - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity + F = reshape(F_lastInc,[9,cells(1),cells(2),cells3]) end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response + homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -343,11 +343,11 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ end if Fdot = utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & + F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) - F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) + F_lastInc = reshape(F,[3,3,cells(1),cells(2),cells3]) - homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- @@ -359,7 +359,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average - rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) + rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc) CHKERRQ(err_PETSc) @@ -530,7 +530,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r ! store fPK field for subsequent FFT forward transform + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier @@ -538,7 +538,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! constructing residual - r = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 810047403..998b94e8a 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -107,7 +107,7 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_mechanical_spectral_polarisation_init - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: P PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI PetscScalar, pointer, dimension(:,:,:,:) :: & @@ -171,10 +171,10 @@ subroutine grid_mechanical_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - allocate(F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate(F_lastInc (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(Fdot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(F_tau_lastInc(3,3,cells(1),cells(2),cells3),source = 0.0_pReal) + allocate(F_tauDot (3,3,cells(1),cells(2),cells3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -183,23 +183,23 @@ subroutine grid_mechanical_spectral_polarisation_init call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 18_pPetscInt, 0_pPetscInt, & ! #dof (2xtensor), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells da,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(da,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(da,err_PETSc) CHKERRQ(err_PETSc) - call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor) + call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (cells x 18, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -241,13 +241,13 @@ subroutine grid_mechanical_spectral_polarisation_init call HDF5_read(F_tau_lastInc,groupHandle,'F_tau_lastInc') elseif (interface_restartInc == 0) then restartRead - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + F_lastInc = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3) ! initialize to identity + F = reshape(F_lastInc,[9,cells(1),cells(2),cells3]) F_tau = 2.0_pReal*F F_tau_lastInc = 2.0_pReal*F_lastInc end if restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3]) ! set starting condition for homogenization_mechanical_response + homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3]) ! set starting condition for homogenization_mechanical_response call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F @@ -379,15 +379,15 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D end if Fdot = utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, & + F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) F_tauDot = utilities_calculateRate(guess, & - F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), Delta_t_old, & + F_tau_lastInc,reshape(F_tau,[3,3,cells(1),cells(2),cells3]), Delta_t_old, & rotation_BC%rotate(F_aimDot,active=.true.)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) - F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) + F_lastInc = reshape(F, [3,3,cells(1),cells(2),cells3]) + F_tau_lastInc = reshape(F_tau,[3,3,cells(1),cells(2),cells3]) - homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3]) + homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3]) end if !-------------------------------------------------------------------------------------------------- @@ -400,12 +400,12 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),& - [9,grid(1),grid(2),grid3]) + [9,cells(1),cells(2),cells3]) if (guess) then F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), & - [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition + [9,cells(1),cells(2),cells3]) ! does not have any average value as boundary condition else - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = math_I3 & + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, & @@ -597,7 +597,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! tensorField_real = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) tensorField_real(1:3,1:3,i,j,k) = & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%alpha*matmul(F(1:3,1:3,i,j,k), & @@ -612,7 +612,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! constructing residual - r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) + r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -629,14 +629,14 @@ subroutine formResidual(in, FandF_tau, & params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r_F !< stress field in disguise + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise call utilities_FFTtensorForward err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual e = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) e = e + 1 r_F(1:3,1:3,i,j,k) = & math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & @@ -648,7 +648,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! calculating curl tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F call utilities_FFTtensorForward err_curl = utilities_curlRMS() diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index f102cf2c1..193d0ecd6 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -101,12 +101,12 @@ subroutine grid_thermal_spectral_init(T_0) !-------------------------------------------------------------------------------------------------- ! init fields - allocate(T_current(grid(1),grid(2),grid3), source=T_0) - allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0) - allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0) + allocate(T_current(cells(1),cells(2),cells3), source=T_0) + allocate(T_lastInc(cells(1),cells(2),cells3), source=T_0) + allocate(T_stagInc(cells(1),cells(2),cells3), source=T_0) ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_thermal_setField(T_0,0.0_pReal,ce) end do; end do; end do @@ -118,23 +118,23 @@ subroutine grid_thermal_spectral_init(T_0) call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc) CHKERRQ(err_PETSc) localK = 0_pPetscInt - localK(worldrank) = int(grid3,pPetscInt) + localK(worldrank) = int(cells3,pPetscInt) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - int(grid(1),pPetscInt),int(grid(2),pPetscInt),int(grid(3),pPetscInt), & ! global grid + int(cells(1),pPetscInt),int(cells(2),pPetscInt),int(cells(3),pPetscInt), & ! global cells 1_pPetscInt, 1_pPetscInt, int(worldsize,pPetscInt), & 1_pPetscInt, 0_pPetscInt, & ! #dof (T, scalar), ghost boundary width (domain overlap) - [int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid + [int(cells(1),pPetscInt)],[int(cells(2),pPetscInt)],localK, & ! local cells thermal_grid,err_PETSc) ! handle, error CHKERRQ(err_PETSc) call DMsetFromOptions(thermal_grid,err_PETSc) CHKERRQ(err_PETSc) call DMsetUp(thermal_grid,err_PETSc) CHKERRQ(err_PETSc) - call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (cells x 1, i.e. every def grad tensor) CHKERRQ(err_PETSc) call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) @@ -198,7 +198,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution) !-------------------------------------------------------------------------------------------------- ! updating thermal state ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) end do; end do; end do @@ -241,7 +241,7 @@ subroutine grid_thermal_spectral_forward(cutBack) call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc) CHKERRQ(err_PETSc) ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce) end do; end do; end do @@ -274,12 +274,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = T_current + scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of temperature field call utilities_FFTvectorBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) end do; end do; end do @@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field call utilities_FFTscalarBackward ce = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & @@ -302,7 +302,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! constructing residual - r = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) + r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3) err_PETSc = 0 end subroutine formResidual @@ -319,7 +319,7 @@ subroutine updateReference() K_ref = 0.0_pReal mu_ref = 0.0_pReal - do ce = 1, product(grid(1:2))*grid3 + do ce = 1, product(cells(1:2))*cells3 K_ref = K_ref + homogenization_K_T(ce) mu_ref = mu_ref + homogenization_mu_T(ce) end do diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 25de88e3f..9e3da5e8d 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -29,9 +29,9 @@ module spectral_utilities include 'fftw3-mpi.f03' !-------------------------------------------------------------------------------------------------- -! grid related information information +! grid related information real(pReal), protected, public :: wgt !< weighting factor 1/Nelems - integer, protected, public :: grid1Red !< grid(1)/2 + integer, protected, public :: grid1Red !< cells(1)/2 real(pReal), protected, public, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence !-------------------------------------------------------------------------------------------------- @@ -201,8 +201,8 @@ subroutine spectral_utilities_init num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc) CHKERRQ(err_PETSc) - grid1Red = grid(1)/2 + 1 - wgt = 1.0/real(product(grid),pReal) + grid1Red = cells(1)/2 + 1 + wgt = 1.0/real(product(cells),pReal) num%memory_efficient = num_grid%get_asInt('memory_efficient', defaultVal=1) > 0 ! ToDo: should be logical in YAML file num%divergence_correction = num_grid%get_asInt('divergence_correction', defaultVal=2) @@ -231,9 +231,9 @@ subroutine spectral_utilities_init enddo elseif (num%divergence_correction == 2) then do j = 1, 3 - if ( j /= int(minloc(geomSize/real(grid,pReal),1)) & - .and. j /= int(maxloc(geomSize/real(grid,pReal),1))) & - scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal) + if ( j /= int(minloc(geomSize/real(cells,pReal),1)) & + .and. j /= int(maxloc(geomSize/real(cells,pReal),1))) & + scaledGeomSize = geomSize/geomSize(j)*real(cells(j),pReal) enddo else scaledGeomSize = geomSize @@ -262,11 +262,11 @@ subroutine spectral_utilities_init !-------------------------------------------------------------------------------------------------- ! MPI allocation - gridFFTW = int(grid,C_INTPTR_T) + gridFFTW = int(cells,C_INTPTR_T) alloc_local = fftw_mpi_local_size_3d(gridFFTW(3), gridFFTW(2), gridFFTW(1)/2 +1, & PETSC_COMM_WORLD, local_K, local_K_offset) - allocate (xi1st (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension - allocate (xi2nd (3,grid1Red,grid(2),grid3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension + allocate (xi1st (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension + allocate (xi2nd (3,grid1Red,cells(2),cells3),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension tensorField = fftw_alloc_complex(tensorSize*alloc_local) call c_f_pointer(tensorField, tensorField_real, [3_C_INTPTR_T,3_C_INTPTR_T, & @@ -327,27 +327,27 @@ subroutine spectral_utilities_init !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - do k = grid3Offset+1, grid3Offset+grid3 + do k = cells3Offset+1, cells3Offset+cells3 k_s(3) = k - 1 - if (k > grid(3)/2 + 1) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 - do j = 1, grid(2) + if (k > cells(3)/2 + 1) k_s(3) = k_s(3) - cells(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + do j = 1, cells(2) k_s(2) = j - 1 - if (j > grid(2)/2 + 1) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 + if (j > cells(2)/2 + 1) k_s(2) = k_s(2) - cells(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do i = 1, grid1Red k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1 - xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s) - where(mod(grid,2)==0 .and. [i,j,k] == grid/2+1 .and. & + xi2nd(1:3,i,j,k-cells3Offset) = utilities_getFreqDerivative(k_s) + where(mod(cells,2)==0 .and. [i,j,k] == cells/2+1 .and. & spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID) ! for even grids, set the Nyquist Freq component to 0.0 - xi1st(1:3,i,j,k-grid3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal) + xi1st(1:3,i,j,k-cells3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal) elsewhere - xi1st(1:3,i,j,k-grid3Offset) = xi2nd(1:3,i,j,k-grid3Offset) + xi1st(1:3,i,j,k-cells3Offset) = xi2nd(1:3,i,j,k-cells3Offset) endwhere enddo; enddo; enddo - if (num%memory_efficient) then ! allocate just single fourth order tensor + if (num%memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal)) else ! precalculation of gamma_hat field - allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) + allocate (gamma_hat(3,3,3,3,grid1Red,cells(2),cells3), source = cmplx(0.0_pReal,0.0_pReal,pReal)) endif end subroutine spectral_utilities_init @@ -373,10 +373,10 @@ subroutine utilities_updateGamma(C) if (.not. num%memory_efficient) then gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A - do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red + do k = cells3Offset+1, cells3Offset+cells3; do j = 1, cells(2); do i = 1, grid1Red if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 do concurrent (l = 1:3, m = 1:3) - xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) + xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) end do do concurrent(l = 1:3, m = 1:3) temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) @@ -387,8 +387,8 @@ subroutine utilities_updateGamma(C) call math_invert(A_inv, err, A) temp33_complex = cmplx(A_inv(1:3,1:3),A_inv(1:3,4:6),pReal) do concurrent(l=1:3, m=1:3, n=1:3, o=1:3) - gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* & - conjg(-xi1st(o,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) + gamma_hat(l,m,n,o,i,j,k-cells3Offset) = temp33_complex(l,n)* & + conjg(-xi1st(o,i,j,k-cells3Offset))*xi1st(m,i,j,k-cells3Offset) end do end if end if @@ -405,7 +405,7 @@ end subroutine utilities_updateGamma !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTtensorForward - tensorField_real(1:3,1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal + tensorField_real(1:3,1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) end subroutine utilities_FFTtensorForward @@ -429,7 +429,7 @@ end subroutine utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTscalarForward - scalarField_real(grid(1)+1:grid1Red*2,:,:) = 0.0_pReal + scalarField_real(cells(1)+1:grid1Red*2,:,:) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) end subroutine utilities_FFTscalarForward @@ -454,7 +454,7 @@ end subroutine utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- subroutine utilities_FFTvectorForward - vectorField_real(1:3,grid(1)+1:grid1Red*2,:,:) = 0.0_pReal + vectorField_real(1:3,cells(1)+1:grid1Red*2,:,:) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) end subroutine utilities_FFTvectorForward @@ -493,8 +493,8 @@ subroutine utilities_fourierGammaConvolution(fieldAim) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) memoryEfficient: if (num%memory_efficient) then - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red - if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red + if (any([i,j,k+cells3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 do concurrent(l = 1:3, m = 1:3) xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) end do @@ -519,7 +519,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) end if end do; end do; end do else memoryEfficient - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red do concurrent(l = 1:3, m = 1:3) temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3,i,j,k) * tensorField_fourier(1:3,1:3,i,j,k)) end do @@ -527,7 +527,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) end do; end do; end do end if memoryEfficient - if (grid3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) + if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) end subroutine utilities_fourierGammaConvolution @@ -544,7 +544,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation - do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1Red + do k = 1, cells3; do j = 1, cells(2) ;do i = 1, grid1Red GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) & / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal) & * sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) @@ -571,7 +571,7 @@ real(pReal) function utilities_divergenceRMS() !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space utilities_divergenceRMS = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2) + do k = 1, cells3; do j = 1, cells(2) do i = 2, grid1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice. utilities_divergenceRMS = utilities_divergenceRMS & + 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k), & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2, i.e. do not take square root and square again @@ -579,7 +579,7 @@ real(pReal) function utilities_divergenceRMS() +sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),& conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2)) enddo - utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1) + utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if cells(1) /= 1) + sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2) & + sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), & @@ -589,7 +589,7 @@ real(pReal) function utilities_divergenceRMS() + sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), & conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2) enddo; enddo - if (grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if (cells(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of cells(1) == 1 call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space @@ -616,7 +616,7 @@ real(pReal) function utilities_curlRMS() ! calculating max curl criterion in Fourier space utilities_curlRMS = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); + do k = 1, cells3; do j = 1, cells(2); do i = 2, grid1Red - 1 do l = 1, 3 curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) & @@ -638,7 +638,7 @@ real(pReal) function utilities_curlRMS() -tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2)) enddo utilities_curlRMS = utilities_curlRMS & - + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1) + + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (DC) does not have a conjugate complex counterpart (if cells(1) /= 1) do l = 1, 3 curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) & -tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3)) @@ -648,13 +648,13 @@ real(pReal) function utilities_curlRMS() -tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2)) enddo utilities_curlRMS = utilities_curlRMS & - + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1) + + sum(curl_fourier%re**2 + curl_fourier%im**2) ! this layer (Nyquist) does not have a conjugate complex counterpart (if cells(1) /= 1) enddo; enddo call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' utilities_curlRMS = sqrt(utilities_curlRMS) * wgt - if (grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1 + if (cells(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of cells(1) == 1 end function utilities_curlRMS @@ -736,7 +736,7 @@ subroutine utilities_fourierScalarGradient() integer :: i, j, k - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k) ! ToDo: no -conjg? enddo; enddo; enddo @@ -750,7 +750,7 @@ subroutine utilities_fourierVectorDivergence() integer :: i, j, k - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red scalarField_fourier(i,j,k) = sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k))) enddo; enddo; enddo @@ -764,7 +764,7 @@ subroutine utilities_fourierVectorGradient() integer :: i, j, k, m, n - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red do m = 1, 3; do n = 1, 3 tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k) enddo; enddo @@ -780,7 +780,7 @@ subroutine utilities_fourierTensorDivergence() integer :: i, j, k - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red + do k = 1, cells3; do j = 1, cells(2); do i = 1,grid1Red vectorField_fourier(:,i,j,k) = matmul(tensorField_fourier(:,:,i,j,k),conjg(-xi1st(:,i,j,k))) enddo; enddo; enddo @@ -795,8 +795,8 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress - real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target + real(pReal), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target real(pReal), intent(in) :: Delta_t !< loading time type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame @@ -810,15 +810,15 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& print'(/,1x,a)', '... evaluating constitutive response ......................................' flush(IO_STDOUT) - homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field + homogenization_F = reshape(F,[3,3,product(cells(1:2))*cells3]) ! set materialpoint target F to estimated field - call homogenization_mechanical_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field + call homogenization_mechanical_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) ! calculate P field if (.not. terminallyIll) & - call homogenization_thermal_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_thermal_response(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) if (.not. terminallyIll) & - call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) + call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3]) - P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) + P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3]) P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' @@ -833,7 +833,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& dPdF_norm_max = 0.0_pReal dPdF_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal) - do i = 1, product(grid(1:2))*grid3 + do i = 1, product(cells(1:2))*cells3 if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i) dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2) @@ -878,16 +878,16 @@ pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) dt !< Delta_t between field0 and field logical, intent(in) :: & heterogeneous !< calculate field of rates - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & field0, & !< data of previous step field !< data of current step - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: & utilities_calculateRate if (heterogeneous) then utilities_calculateRate = (field-field0) / dt else - utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid3) + utilities_calculateRate = spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3) endif end function utilities_calculateRate @@ -901,12 +901,12 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim) real(pReal), intent(in) :: & Delta_t !< Delta_t of current step - real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: & field_lastInc, & !< initial field rate !< rate by which to forward real(pReal), intent(in), optional, dimension(3,3) :: & aim !< average field value aim - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: & + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: & utilities_forwardField real(pReal), dimension(3,3) :: fieldDiff !< - aim integer(MPI_INTEGER_KIND) :: err_MPI @@ -918,7 +918,7 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' fieldDiff = fieldDiff - aim utilities_forwardField = utilities_forwardField - & - spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid3) + spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3) endif end function utilities_forwardField @@ -939,34 +939,34 @@ pure function utilities_getFreqDerivative(k_s) utilities_getFreqDerivative = cmplx(0.0_pReal, TAU*real(k_s,pReal)/geomSize,pReal) case (DERIVATIVE_CENTRAL_DIFF_ID) - utilities_getFreqDerivative = cmplx(0.0_pReal, sin(TAU*real(k_s,pReal)/real(grid,pReal)), pReal)/ & - cmplx(2.0_pReal*geomSize/real(grid,pReal), 0.0_pReal, pReal) + utilities_getFreqDerivative = cmplx(0.0_pReal, sin(TAU*real(k_s,pReal)/real(cells,pReal)), pReal)/ & + cmplx(2.0_pReal*geomSize/real(cells,pReal), 0.0_pReal, pReal) case (DERIVATIVE_FWBW_DIFF_ID) utilities_getFreqDerivative(1) = & - cmplx(cos(TAU*real(k_s(1),pReal)/real(grid(1),pReal)) - 1.0_pReal, & - sin(TAU*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(TAU*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(TAU*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(TAU*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(TAU*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(1)/real(grid(1),pReal), 0.0_pReal, pReal) + cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(1)/real(cells(1),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(2) = & - cmplx(cos(TAU*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(TAU*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(TAU*real(k_s(2),pReal)/real(grid(2),pReal)) - 1.0_pReal, & - sin(TAU*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(TAU*real(k_s(3),pReal)/real(grid(3),pReal)) + 1.0_pReal, & - sin(TAU*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(2)/real(grid(2),pReal), 0.0_pReal, pReal) + cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(2)/real(cells(2),pReal), 0.0_pReal, pReal) utilities_getFreqDerivative(3) = & - cmplx(cos(TAU*real(k_s(1),pReal)/real(grid(1),pReal)) + 1.0_pReal, & - sin(TAU*real(k_s(1),pReal)/real(grid(1),pReal)), pReal)* & - cmplx(cos(TAU*real(k_s(2),pReal)/real(grid(2),pReal)) + 1.0_pReal, & - sin(TAU*real(k_s(2),pReal)/real(grid(2),pReal)), pReal)* & - cmplx(cos(TAU*real(k_s(3),pReal)/real(grid(3),pReal)) - 1.0_pReal, & - sin(TAU*real(k_s(3),pReal)/real(grid(3),pReal)), pReal)/ & - cmplx(4.0_pReal*geomSize(3)/real(grid(3),pReal), 0.0_pReal, pReal) + cmplx(cos(TAU*real(k_s(1),pReal)/real(cells(1),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(1),pReal)/real(cells(1),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(2),pReal)/real(cells(2),pReal)) + 1.0_pReal, & + sin(TAU*real(k_s(2),pReal)/real(cells(2),pReal)), pReal)* & + cmplx(cos(TAU*real(k_s(3),pReal)/real(cells(3),pReal)) - 1.0_pReal, & + sin(TAU*real(k_s(3),pReal)/real(cells(3),pReal)), pReal)/ & + cmplx(4.0_pReal*geomSize(3)/real(cells(3),pReal), 0.0_pReal, pReal) end select end function utilities_getFreqDerivative @@ -979,10 +979,10 @@ end function utilities_getFreqDerivative !-------------------------------------------------------------------------------------------------- subroutine utilities_updateCoords(F) - real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F - real(pReal), dimension(3, grid(1),grid(2),grid3) :: IPcoords - real(pReal), dimension(3, grid(1),grid(2),grid3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI) - real(pReal), dimension(3, grid(1)+1,grid(2)+1,grid3+1) :: nodeCoords + real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: F + real(pReal), dimension(3, cells(1),cells(2),cells3) :: IPcoords + real(pReal), dimension(3, cells(1),cells(2),cells3+2) :: IPfluct_padded ! Fluctuations of cell center displacement (padded along z for MPI) + real(pReal), dimension(3, cells(1)+1,cells(2)+1,cells3+1) :: nodeCoords integer :: & i,j,k,n, & c @@ -1010,14 +1010,14 @@ subroutine utilities_updateCoords(F) 1, 1, 1, & 0, 1, 1 ], [3,8]) - step = geomSize/real(grid, pReal) + step = geomSize/real(cells, pReal) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space to get fluctuations of cell center discplacements - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F call utilities_FFTtensorForward() - do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red - if (any([i,j,k+grid3Offset] /= 1)) then + do k = 1, cells3; do j = 1, cells(2); do i = 1, grid1Red + if (any([i,j,k+cells3Offset] /= 1)) then vectorField_fourier(1:3,i,j,k) = matmul(tensorField_fourier(1:3,1:3,i,j,k),xi2nd(1:3,i,j,k)) & / sum(conjg(-xi2nd(1:3,i,j,k))*xi2nd(1:3,i,j,k)) * cmplx(wgt,0.0,pReal) else @@ -1029,25 +1029,25 @@ subroutine utilities_updateCoords(F) !-------------------------------------------------------------------------------------------------- ! average F - if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt + if (cells3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) - IPfluct_padded(1:3,1:grid(1),1:grid(2),2:grid3+1) = vectorField_real(1:3,1:grid(1),1:grid(2),1:grid3) - c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer + IPfluct_padded(1:3,1:cells(1),1:cells(2),2:cells3+1) = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3) + c = product(shape(IPfluct_padded(:,:,:,1))) !< amount of data to transfer rank_t = modulo(worldrank+1_MPI_INTEGER_KIND,worldsize) rank_b = modulo(worldrank-1_MPI_INTEGER_KIND,worldsize) ! send bottom layer to process below call MPI_Isend(IPfluct_padded(:,:,:,2), c,MPI_DOUBLE,rank_b,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(1),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_Irecv(IPfluct_padded(:,:,:,grid3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI) + call MPI_Irecv(IPfluct_padded(:,:,:,cells3+2),c,MPI_DOUBLE,rank_t,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(2),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' ! send top layer to process above - call MPI_Isend(IPfluct_padded(:,:,:,grid3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI) + call MPI_Isend(IPfluct_padded(:,:,:,cells3+1),c,MPI_DOUBLE,rank_t,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(3),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_Irecv(IPfluct_padded(:,:,:,1), c,MPI_DOUBLE,rank_b,1_MPI_INTEGER_KIND,MPI_COMM_WORLD,request(4),err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' @@ -1063,24 +1063,24 @@ subroutine utilities_updateCoords(F) !-------------------------------------------------------------------------------------------------- ! calculate nodal displacements nodeCoords = 0.0_pReal - do k = 0,grid3; do j = 0,grid(2); do i = 0,grid(1) - nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+grid3Offset],pReal))) + do k = 0,cells3; do j = 0,cells(2); do i = 0,cells(1) + nodeCoords(1:3,i+1,j+1,k+1) = matmul(Favg,step*(real([i,j,k+cells3Offset],pReal))) averageFluct: do n = 1,8 me = [i+neighbor(1,n),j+neighbor(2,n),k+neighbor(3,n)] nodeCoords(1:3,i+1,j+1,k+1) = nodeCoords(1:3,i+1,j+1,k+1) & - + IPfluct_padded(1:3,modulo(me(1)-1,grid(1))+1,modulo(me(2)-1,grid(2))+1,me(3)+1)*0.125_pReal + + IPfluct_padded(1:3,modulo(me(1)-1,cells(1))+1,modulo(me(2)-1,cells(2))+1,me(3)+1)*0.125_pReal enddo averageFluct enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! calculate cell center displacements - do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1) + do k = 1,cells3; do j = 1,cells(2); do i = 1,cells(1) IPcoords(1:3,i,j,k) = vectorField_real(1:3,i,j,k) & - + matmul(Favg,step*(real([i,j,k+grid3Offset],pReal)-0.5_pReal)) + + matmul(Favg,step*(real([i,j,k+cells3Offset],pReal)-0.5_pReal)) enddo; enddo; enddo - call discretization_setNodeCoords(reshape(NodeCoords,[3,(grid(1)+1)*(grid(2)+1)*(grid3+1)])) - call discretization_setIPcoords (reshape(IPcoords, [3,grid(1)*grid(2)*grid3])) + call discretization_setNodeCoords(reshape(NodeCoords,[3,(cells(1)+1)*(cells(2)+1)*(cells3+1)])) + call discretization_setIPcoords (reshape(IPcoords, [3,cells(1)*cells(2)*cells3])) end subroutine utilities_updateCoords From 762f93d724cdebc3a7b1e667d3de83b8c149a163 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 15:30:59 +0100 Subject: [PATCH 06/26] following naming convention --- src/grid/DAMASK_grid.f90 | 2 +- src/grid/grid_mech_FEM.f90 | 2 +- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarisation.f90 | 2 +- src/grid/spectral_utilities.f90 | 6 +- src/lattice.f90 | 8 +-- src/material.f90 | 2 +- src/phase_mechanical.f90 | 2 +- src/phase_mechanical_plastic_nonlocal.f90 | 2 +- src/rotations.f90 | 58 ++++++++++---------- 10 files changed, 43 insertions(+), 43 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 6377748be..901132d99 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -32,7 +32,7 @@ program DAMASK_grid implicit none type :: tLoadCase - type(rotation) :: rot !< rotation of BC + type(tRotation) :: rot !< rotation of BC type(tBoundaryCondition) :: stress, & !< stress BC deformation !< deformation BC (dot_F, F, or L) real(pReal) :: t, & !< length of increment diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index cc082d2f5..6ad035fc2 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -339,7 +339,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: & diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 2afe9edc7..fa2e17bd9 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -310,7 +310,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_ type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: F diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 998b94e8a..ee87c77a6 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -342,7 +342,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 9e3da5e8d..d95580145 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -86,7 +86,7 @@ module spectral_utilities type, public :: tSolutionParams real(pReal), dimension(3,3) :: stress_BC logical, dimension(3,3) :: stress_mask - type(rotation) :: rotation_BC + type(tRotation) :: rotation_BC real(pReal) :: Delta_t end type tSolutionParams @@ -666,7 +666,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C) real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance real(pReal), intent(in), dimension(3,3,3,3) :: C !< current average stiffness - type(rotation), intent(in) :: rot_BC !< rotation of load frame + type(tRotation), intent(in) :: rot_BC !< rotation of load frame logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC integer :: i, j @@ -798,7 +798,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& real(pReal), intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P !< PK stress real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: F !< deformation gradient target real(pReal), intent(in) :: Delta_t !< loading time - type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame + type(tRotation), intent(in), optional :: rotation_BC !< rotation of load frame integer :: i diff --git a/src/lattice.f90 b/src/lattice.f90 index e18c71edf..2c2840bca 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -498,7 +498,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA) real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(3,3,sum(Ntwin)):: coordinateSystem - type(rotation) :: R + type(tRotation) :: R integer :: i @@ -537,7 +537,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S - type(rotation) :: R + type(tRotation) :: R real(pReal) :: a_bcc, a_fcc, cOverA_trans integer :: i @@ -599,7 +599,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system real(pReal), dimension(3) :: direction, normal, np - type(rotation) :: R + type(tRotation) :: R integer :: i @@ -1976,7 +1976,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) a_bcc, & !< lattice parameter a for bcc target lattice a_fcc !< lattice parameter a for fcc parent lattice - type(rotation) :: & + type(tRotation) :: & R, & !< Pitsch rotation B !< Rotation of fcc to Bain coordinate system real(pReal), dimension(3,3) :: & diff --git a/src/material.f90 b/src/material.f90 index 1e3a4b4ec..b2beb0635 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -18,7 +18,7 @@ module material private type :: tRotationContainer - type(Rotation), dimension(:), allocatable :: data + type(tRotation), dimension(:), allocatable :: data end type type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 1917d81c9..bef8dbb18 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -958,7 +958,7 @@ subroutine crystallite_results(group,ph) !-------------------------------------------------------------------------------------------------- function to_quaternion(dataset) - type(rotation), dimension(:), intent(in) :: dataset + type(tRotation), dimension(:), intent(in) :: dataset real(pReal), dimension(4,size(dataset,1)) :: to_quaternion integer :: i diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 1715af2d9..0dc3dc4bd 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -1395,7 +1395,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e) nThresholdValues logical, dimension(param(ph)%sum_N_sl) :: & belowThreshold - type(rotation) :: mis + type(tRotation) :: mis associate(prm => param(ph)) diff --git a/src/rotations.f90 b/src/rotations.f90 index ec000af60..fcc7da0d2 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -55,7 +55,7 @@ module rotations real(pReal), parameter :: P = -1.0_pReal !< parameter for orientation conversion. - type, public :: rotation + type, public :: tRotation real(pReal), dimension(4) :: q contains procedure, public :: asQuaternion @@ -78,7 +78,7 @@ module rotations procedure, public :: rotStiffness procedure, public :: misorientation procedure, public :: standardize - end type rotation + end type tRotation real(pReal), parameter :: & PREF = sqrt(6.0_pReal/PI), & @@ -117,8 +117,8 @@ end subroutine rotations_init !--------------------------------------------------------------------------------------------------- pure function asQuaternion(self) - class(rotation), intent(in) :: self - real(pReal), dimension(4) :: asQuaternion + class(tRotation), intent(in) :: self + real(pReal), dimension(4) :: asQuaternion asQuaternion = self%q @@ -126,8 +126,8 @@ end function asQuaternion !--------------------------------------------------------------------------------------------------- pure function asEulers(self) - class(rotation), intent(in) :: self - real(pReal), dimension(3) :: asEulers + class(tRotation), intent(in) :: self + real(pReal), dimension(3) :: asEulers asEulers = qu2eu(self%q) @@ -135,8 +135,8 @@ end function asEulers !--------------------------------------------------------------------------------------------------- pure function asAxisAngle(self) - class(rotation), intent(in) :: self - real(pReal), dimension(4) :: asAxisAngle + class(tRotation), intent(in) :: self + real(pReal), dimension(4) :: asAxisAngle asAxisAngle = qu2ax(self%q) @@ -144,8 +144,8 @@ end function asAxisAngle !--------------------------------------------------------------------------------------------------- pure function asMatrix(self) - class(rotation), intent(in) :: self - real(pReal), dimension(3,3) :: asMatrix + class(tRotation), intent(in) :: self + real(pReal), dimension(3,3) :: asMatrix asMatrix = qu2om(self%q) @@ -153,8 +153,8 @@ end function asMatrix !--------------------------------------------------------------------------------------------------- pure function asRodrigues(self) - class(rotation), intent(in) :: self - real(pReal), dimension(4) :: asRodrigues + class(tRotation), intent(in) :: self + real(pReal), dimension(4) :: asRodrigues asRodrigues = qu2ro(self%q) @@ -162,8 +162,8 @@ end function asRodrigues !--------------------------------------------------------------------------------------------------- pure function asHomochoric(self) - class(rotation), intent(in) :: self - real(pReal), dimension(3) :: asHomochoric + class(tRotation), intent(in) :: self + real(pReal), dimension(3) :: asHomochoric asHomochoric = qu2ho(self%q) @@ -174,7 +174,7 @@ end function asHomochoric !--------------------------------------------------------------------------------------------------- subroutine fromQuaternion(self,qu) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(4), intent(in) :: qu if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion') @@ -185,7 +185,7 @@ end subroutine fromQuaternion !--------------------------------------------------------------------------------------------------- subroutine fromEulers(self,eu,degrees) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(3), intent(in) :: eu logical, intent(in), optional :: degrees @@ -206,7 +206,7 @@ end subroutine fromEulers !--------------------------------------------------------------------------------------------------- subroutine fromAxisAngle(self,ax,degrees,P) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(4), intent(in) :: ax logical, intent(in), optional :: degrees integer, intent(in), optional :: P @@ -236,7 +236,7 @@ end subroutine fromAxisAngle !--------------------------------------------------------------------------------------------------- subroutine fromMatrix(self,om) - class(rotation), intent(out) :: self + class(tRotation), intent(out) :: self real(pReal), dimension(3,3), intent(in) :: om if (dNeq(math_det33(om),1.0_pReal,tol=1.0e-5_pReal)) & @@ -253,10 +253,10 @@ end subroutine fromMatrix !--------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) - type(rotation) :: rRot - class(rotation), intent(in) :: self,R + type(tRotation) :: rRot + class(tRotation), intent(in) :: self,R - rRot = rotation(multiply_quaternion(self%q,R%q)) + rRot = tRotation(multiply_quaternion(self%q,R%q)) call rRot%standardize() end function rotRot__ @@ -267,7 +267,7 @@ end function rotRot__ !--------------------------------------------------------------------------------------------------- pure elemental subroutine standardize(self) - class(rotation), intent(inout) :: self + class(tRotation), intent(inout) :: self if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q @@ -281,7 +281,7 @@ end subroutine standardize pure function rotVector(self,v,active) result(vRot) real(pReal), dimension(3) :: vRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(3) :: v logical, intent(in), optional :: active @@ -317,7 +317,7 @@ end function rotVector pure function rotTensor2(self,T,active) result(tRot) real(pReal), dimension(3,3) :: tRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(3,3) :: T logical, intent(in), optional :: active @@ -346,7 +346,7 @@ end function rotTensor2 pure function rotTensor4(self,T,active) result(tRot) real(pReal), dimension(3,3,3,3) :: tRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(3,3,3,3) :: T logical, intent(in), optional :: active @@ -378,7 +378,7 @@ end function rotTensor4 pure function rotStiffness(self,C,active) result(cRot) real(pReal), dimension(6,6) :: cRot - class(rotation), intent(in) :: self + class(tRotation), intent(in) :: self real(pReal), intent(in), dimension(6,6) :: C logical, intent(in), optional :: active @@ -415,8 +415,8 @@ end function rotStiffness !--------------------------------------------------------------------------------------------------- pure elemental function misorientation(self,other) - type(rotation) :: misorientation - class(rotation), intent(in) :: self, other + type(tRotation) :: misorientation + class(tRotation), intent(in) :: self, other misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) @@ -1415,7 +1415,7 @@ end function conjugate_quaternion !-------------------------------------------------------------------------------------------------- subroutine selfTest() - type(rotation) :: R + type(tRotation) :: R real(pReal), dimension(4) :: qu, ax, ro real(pReal), dimension(3) :: x, eu, ho, v3 real(pReal), dimension(3,3) :: om, t33 From f70df11b673f2f1b400ccb855f8010389852bda1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 29 Jan 2022 15:59:14 +0100 Subject: [PATCH 07/26] consistent reporting --- python/damask/_grid.py | 4 ++-- src/grid/discretization_grid.f90 | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 94f7eefd0..dadfcced5 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -62,8 +62,8 @@ class Grid: mat_max = np.nanmax(self.material) mat_N = self.N_materials return util.srepr([ - f'cells : {util.srepr(self.cells, " x ")}', - f'size : {util.srepr(self.size, " x ")} / m³', + f'cells: {util.srepr(self.cells, " × ")}', + f'size: {util.srepr(self.size, " × ")} / m³', f'origin: {util.srepr(self.origin," ")} / m', f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else f' (min: {mat_min}, max: {mat_max})') diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 67f735a0d..6e8c154ef 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -96,9 +96,9 @@ subroutine discretization_grid_init(restart) call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - print'(/,1x,a,3(i12,1x))', 'cells a b c: ', cells - print '(1x,a,3(es12.5,1x))', 'size x y z: ', geomSize - print '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin + print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3) + print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' / m³' + print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' / m' if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)') From a673abb413579ecd7f54f2fb966511c58e56a290 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 30 Jan 2022 06:24:50 +0100 Subject: [PATCH 08/26] indicate that this mapping should be used only for C --- src/math.f90 | 38 ++++++++++++++++---------------- src/phase_mechanical_elastic.f90 | 10 ++++----- src/rotations.f90 | 3 ++- 3 files changed, 26 insertions(+), 25 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 0343d190a..5c6b10fdd 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -883,7 +883,7 @@ end function math_Voigt6to33_strain !-------------------------------------------------------------------------------------------------- -!> @brief Convert 3x3 tensor into 6 Voigt stress vector. +!> @brief Convert 3x3 stress tensor into 6 Voigt vector. !-------------------------------------------------------------------------------------------------- pure function math_33toVoigt6_stress(sigma) result(sigma_tilde) @@ -898,7 +898,7 @@ end function math_33toVoigt6_stress !-------------------------------------------------------------------------------------------------- -!> @brief Convert 3x3 tensor into 6 Voigt strain vector. +!> @brief Convert 3x3 stress tensor into 6 Voigt vector. !-------------------------------------------------------------------------------------------------- pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde) @@ -914,48 +914,48 @@ end function math_33toVoigt6_strain !-------------------------------------------------------------------------------------------------- -!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix. +!> @brief Convert 6x6 Voigt stiffness matrix into symmetric 3x3x3x3 tensor. !-------------------------------------------------------------------------------------------------- -pure function math_Voigt66to3333(m66) +pure function math_Voigt66to3333_stiffness(C_tilde) result(C) - real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 - real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix + real(pReal), dimension(3,3,3,3) :: C + real(pReal), dimension(6,6), intent(in) :: C_tilde integer :: i,j do i=1,6; do j=1,6 - math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) - math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) - math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) - math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = m66(i,j) + C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j) + C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = C_tilde(i,j) + C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j) + C(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(2,j),MAPVOIGT(1,j)) = C_tilde(i,j) end do; end do -end function math_Voigt66to3333 +end function math_Voigt66to3333_stiffness !-------------------------------------------------------------------------------------------------- -!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix. +!> @brief Convert 3x3x3x3 stiffness tensor into 6x6 Voigt matrix. !-------------------------------------------------------------------------------------------------- -pure function math_3333toVoigt66(m3333) +pure function math_3333toVoigt66_stiffness(C) result(C_tilde) - real(pReal), dimension(6,6) :: math_3333toVoigt66 - real(pReal), dimension(3,3,3,3), intent(in) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check) + real(pReal), dimension(6,6) :: C_tilde + real(pReal), dimension(3,3,3,3), intent(in) :: C integer :: i,j #ifndef __INTEL_COMPILER do concurrent(i=1:6, j=1:6) - math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) end do #else do i=1,6; do j=1,6 - math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + C_tilde(i,j) = C(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) end do; end do #endif -end function math_3333toVoigt66 +end function math_3333toVoigt66_stiffness !-------------------------------------------------------------------------------------------------- @@ -1344,7 +1344,7 @@ subroutine selfTest if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' - if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_3333toVoigt66_stiffness(math_Voigt66to3333_stiffness(t66)),t66,1.0e-15_pReal))) & error stop 'math_3333toVoigt66/math_Voigt66to3333' call random_number(v6) diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 8972367d5..0922cf41d 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -52,20 +52,20 @@ module subroutine elastic_init(phases) prm%C_11(1) = elastic%get_asFloat('C_11') prm%C_11(2) = elastic%get_asFloat('C_11,T', defaultVal=0.0_pReal) prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal) - + prm%C_12(1) = elastic%get_asFloat('C_12') prm%C_12(2) = elastic%get_asFloat('C_12,T', defaultVal=0.0_pReal) prm%C_12(3) = elastic%get_asFloat('C_12,T^2',defaultVal=0.0_pReal) - + prm%C_44(1) = elastic%get_asFloat('C_44') prm%C_44(2) = elastic%get_asFloat('C_44,T', defaultVal=0.0_pReal) prm%C_44(3) = elastic%get_asFloat('C_44,T^2',defaultVal=0.0_pReal) - + if (any(phase_lattice(ph) == ['hP','tI'])) then prm%C_13(1) = elastic%get_asFloat('C_13') prm%C_13(2) = elastic%get_asFloat('C_13,T', defaultVal=0.0_pReal) prm%C_13(3) = elastic%get_asFloat('C_13,T^2',defaultVal=0.0_pReal) - + prm%C_33(1) = elastic%get_asFloat('C_33') prm%C_33(2) = elastic%get_asFloat('C_33,T', defaultVal=0.0_pReal) prm%C_33(3) = elastic%get_asFloat('C_33,T^2',defaultVal=0.0_pReal) @@ -200,7 +200,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & C66 = phase_damage_C66(phase_homogenizedC66(ph,en),ph,en) - C = math_Voigt66to3333(C66) + C = math_Voigt66to3333_stiffness(C66) E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_Voigt6to33_stress(matmul(C66,math_33toVoigt6_strain(matmul(matmul(transpose(Fi),E),Fi))))!< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration diff --git a/src/rotations.f90 b/src/rotations.f90 index fcc7da0d2..27bf8c227 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -1503,7 +1503,8 @@ subroutine selfTest() call random_number(C) C = C+transpose(C) - if (any(dNeq(R%rotStiffness(C),math_3333toVoigt66(R%rotate(math_Voigt66to3333(C))),1.0e-12_pReal))) & + if (any(dNeq(R%rotStiffness(C), & + math_3333toVoigt66_stiffness(R%rotate(math_Voigt66to3333_stiffness(C))),1.0e-12_pReal))) & error stop 'rotStiffness' call R%fromQuaternion(qu * (1.0_pReal + merge(+5.e-9_pReal,-5.e-9_pReal, mod(i,2) == 0))) ! allow reasonable tolerance for ASCII/YAML From 08d3e22bcda2cc2de5dd657ea74a570071f5bb30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Feb 2022 08:35:00 +0100 Subject: [PATCH 09/26] size is known --- src/phase_mechanical.f90 | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index bef8dbb18..fb836e2e6 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -40,8 +40,6 @@ submodule(phase) mechanical integer(kind(PLASTIC_undefined_ID)), dimension(:), allocatable :: & phase_plasticity !< plasticity of each phase - integer :: phase_plasticity_maxSizeDotState - interface module subroutine eigen_init(phases) @@ -296,8 +294,6 @@ module subroutine mechanical_init(phases) do ph = 1,phases%length plasticState(ph)%state0 = plasticState(ph)%state enddo - phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -588,9 +584,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul real(pReal), intent(in),dimension(:) :: subState0 real(pReal), intent(in) :: Delta_t integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop logical :: & broken @@ -601,9 +597,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul sizeDotState real(pReal) :: & zeta - real(pReal), dimension(phase_plasticity_maxSizeDotState) :: & + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & r ! state residuum - real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: & + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: & dotState @@ -615,7 +611,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul sizeDotState = plasticState(ph)%sizeDotState plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t + + plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t iteration: do NiterationState = 1, num%nState @@ -729,7 +725,8 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip ph, & en, & sizeDotState - real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + r ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) @@ -740,7 +737,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip sizeDotState = plasticState(ph)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t + r = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t plasticState(ph)%state(1:sizeDotState,en) = subState0 & + plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t @@ -753,7 +750,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) return - broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, & + broken = .not. converged(r + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState)) @@ -847,7 +844,8 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D ph, & en, & sizeDotState - real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: & + plastic_RKdotState ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) From 1871dcd75ae1fe3f87fb5a3287e6861023f08fe7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Feb 2022 09:03:39 +0100 Subject: [PATCH 10/26] small polishing --- src/phase_mechanical.f90 | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index fb836e2e6..01ca685c2 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -628,12 +628,12 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul dotState(1:sizeDotState,2)) plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta & + dotState(1:sizeDotState,1) * (1.0_pReal - zeta) - r(1:sizeDotState) = plasticState(ph)%state(1:sizeDotState,en) & - - subState0 & - - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t - plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) & - - r(1:sizeDotState) - if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then + r = plasticState(ph)%state(1:sizeDotState,en) & + - subState0 & + - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r + + if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then broken = plastic_deltaState(ph,en) exit iteration endif @@ -648,19 +648,18 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul !-------------------------------------------------------------------------------------------------- real(pReal) pure function damper(omega_0,omega_1,omega_2) - real(pReal), dimension(:), intent(in) :: & - omega_0, omega_1, omega_2 + real(pReal), dimension(:), intent(in) :: & + omega_0, omega_1, omega_2 - real(pReal) :: dot_prod12, dot_prod22 + real(pReal) :: dot_prod12, dot_prod22 - dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2) - dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2) - if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif + dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2) + dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2) + + damper = merge(0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22), & + 1.0_pReal, & + min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) end function damper From 49ca62f3af1cb794652667306c4f62454c5c817b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Feb 2022 11:22:27 +0100 Subject: [PATCH 11/26] only Chuck Norris can divide by zero --- src/phase_mechanical.f90 | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 01ca685c2..6645cd36b 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -657,9 +657,11 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul dot_prod12 = dot_product(omega_0-omega_1, omega_1-omega_2) dot_prod22 = dot_product(omega_1-omega_2, omega_1-omega_2) - damper = merge(0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22), & - 1.0_pReal, & - min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) + if (min(dot_product(omega_0,omega_1),dot_prod12) < 0.0_pReal .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif end function damper From 4eb9c3dea724d2e671da9cf7a8e3378f9938c15d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Feb 2022 09:23:03 +0100 Subject: [PATCH 12/26] correct name --- src/phase_mechanical.f90 | 12 ++++++------ src/phase_mechanical_plastic.f90 | 2 +- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 6645cd36b..68b069aab 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -600,7 +600,7 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & r ! state residuum real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: & - dotState + dotState_last ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) @@ -615,8 +615,8 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul iteration: do NiterationState = 1, num%nState - dotState(1:sizeDotState,2) = merge(dotState(1:sizeDotState,1),0.0, nIterationState > 1) - dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,en) + dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) + dotState_last(1:sizeDotState,1) = plasticState(ph)%dotState(:,en) broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) exit iteration @@ -624,10 +624,10 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul broken = plastic_dotState(Delta_t, co,ip,el,ph,en) if(broken) exit iteration - zeta = damper(plasticState(ph)%dotState(:,en),dotState(1:sizeDotState,1),& - dotState(1:sizeDotState,2)) + zeta = damper(plasticState(ph)%dotState(:,en),dotState_last(1:sizeDotState,1),& + dotState_last(1:sizeDotState,2)) plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta & - + dotState(1:sizeDotState,1) * (1.0_pReal - zeta) + + dotState_last(1:sizeDotState,1) * (1.0_pReal - zeta) r = plasticState(ph)%state(1:sizeDotState,en) & - subState0 & - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 0a24b0cbf..c83db109d 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -285,7 +285,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & do i=1,3; do j=1,3 dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) end do; end do end if From a14735f3db53ba65857d61231eb31631c4ad3873 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Feb 2022 09:39:13 +0100 Subject: [PATCH 13/26] first step towards removal of global dot state less state => better code --- src/phase_mechanical.f90 | 85 ++++++++++++++++++-------------- src/phase_mechanical_plastic.f90 | 7 +-- 2 files changed, 51 insertions(+), 41 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 68b069aab..857251abb 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -79,7 +79,7 @@ submodule(phase) mechanical en end subroutine plastic_isotropic_LiAndItsTangent - module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) + module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point @@ -88,7 +88,8 @@ submodule(phase) mechanical en real(pReal), intent(in) :: & subdt !< timestep - logical :: broken + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + dotState end function plastic_dotState module function plastic_deltaState(ph, en) result(broken) @@ -598,39 +599,39 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul real(pReal) :: & zeta real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & - r ! state residuum + r, & ! state residuum + dotState real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,2) :: & dotState_last ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState - plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t + plasticState(ph)%state(1:sizeDotState,en) = subState0 + dotState * Delta_t iteration: do NiterationState = 1, num%nState dotState_last(1:sizeDotState,2) = merge(dotState_last(1:sizeDotState,1),0.0, nIterationState > 1) - dotState_last(1:sizeDotState,1) = plasticState(ph)%dotState(:,en) + dotState_last(1:sizeDotState,1) = dotState broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) exit iteration - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) exit iteration + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) exit iteration - zeta = damper(plasticState(ph)%dotState(:,en),dotState_last(1:sizeDotState,1),& - dotState_last(1:sizeDotState,2)) - plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) * zeta & - + dotState_last(1:sizeDotState,1) * (1.0_pReal - zeta) + zeta = damper(dotState,dotState_last(1:sizeDotState,1),dotState_last(1:sizeDotState,2)) + dotState = dotState * zeta & + + dotState_last(1:sizeDotState,1) * (1.0_pReal - zeta) r = plasticState(ph)%state(1:sizeDotState,en) & - subState0 & - - plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t + - dotState * Delta_t plasticState(ph)%state(1:sizeDotState,en) = plasticState(ph)%state(1:sizeDotState,en) - r if (converged(r,plasticState(ph)%state(1:sizeDotState,en),plasticState(ph)%atol(1:sizeDotState))) then @@ -683,6 +684,8 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res logical :: & broken + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + dotState integer :: & ph, & en, & @@ -691,13 +694,14 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState(1:sizeDotState,en) * Delta_t + + dotState * Delta_t broken = plastic_deltaState(ph,en) if(broken) return @@ -727,20 +731,22 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip en, & sizeDotState real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & - r + r, & + dotState ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState - r = - plasticState(ph)%dotstate(1:sizeDotState,en) * 0.5_pReal * Delta_t + r = - dotState * 0.5_pReal * Delta_t plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotstate(1:sizeDotState,en) * Delta_t + + dotState * Delta_t broken = plastic_deltaState(ph,en) if(broken) return @@ -748,10 +754,10 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip broken = integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) if(broken) return - broken = plastic_dotState(Delta_t, co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return - broken = .not. converged(r + 0.5_pReal * plasticState(ph)%dotState(:,en) * Delta_t, & + broken = .not. converged(r + 0.5_pReal * dotState * Delta_t, & plasticState(ph)%state(1:sizeDotState,en), & plasticState(ph)%atol(1:sizeDotState)) @@ -845,45 +851,48 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D ph, & en, & sizeDotState + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + dotState real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState,size(B)) :: & plastic_RKdotState ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) + broken = .true. - broken = plastic_dotState(Delta_t,co,ip,el,ph,en) - if(broken) return + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) return sizeDotState = plasticState(ph)%sizeDotState do stage = 1, size(A,1) - plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,en) - plasticState(ph)%dotState(:,en) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + plastic_RKdotState(1:sizeDotState,stage) = dotState + dotState = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) do n = 2, stage - plasticState(ph)%dotState(:,en) = plasticState(ph)%dotState(:,en) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + dotState = dotState & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) enddo plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t + + dotState * Delta_t broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) if(broken) exit - broken = plastic_dotState(Delta_t*C(stage),co,ip,el,ph,en) - if(broken) exit + dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + if (any(IEEE_is_NaN(dotState))) exit enddo if(broken) return - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,en) - plasticState(ph)%dotState(:,en) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plastic_RKdotState(1:sizeDotState,size(B)) = dotState + dotState = matmul(plastic_RKdotState,B) plasticState(ph)%state(1:sizeDotState,en) = subState0 & - + plasticState(ph)%dotState (1:sizeDotState,en) * Delta_t + + dotState * Delta_t if(present(DB)) & broken = .not. converged(matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) * Delta_t, & diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index c83db109d..6d2aac827 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -296,7 +296,7 @@ end subroutine plastic_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) +module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) integer, intent(in) :: & co, & !< component-ID of integration point @@ -308,7 +308,8 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) subdt !< timestep real(pReal), dimension(3,3) :: & Mp - logical :: broken + real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + dotState if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then @@ -337,7 +338,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken) end select plasticType end if - broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en))) + dotState = plasticState(ph)%dotState(:,en) end function plastic_dotState From 1eda18b9e55a33f470dd952c842feeec5814902f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 1 Feb 2022 14:47:18 +0100 Subject: [PATCH 14/26] fine tuning optimization flag using AGRESSIVE with gcc reduced runtime of the small example on my laptop from 2:00 min to 0:44 min --- cmake/Compiler-GNU.cmake | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index cafcfabca..5e2bf2ce6 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -12,9 +12,9 @@ endif () if (OPTIMIZATION STREQUAL "OFF") set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") - set (OPTIMIZATION_FLAGS "-O2 -mtune=generic") + set (OPTIMIZATION_FLAGS "-O2 -mtune=generic -flto") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") - set (OPTIMIZATION_FLAGS "-O3 -march=native -ffast-math -funroll-loops -ftree-vectorize") + set (OPTIMIZATION_FLAGS "-O3 -march=native -mtune=native -ffast-math -funroll-loops -ftree-vectorize -flto") endif () set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) From 0807112a71085a7832f1969225712361361ebd66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 3 Feb 2022 14:37:44 +0100 Subject: [PATCH 15/26] =?UTF-8?q?new=20tests=20with=20T=5Fref=20=3D=20293.?= =?UTF-8?q?15K=20(20=C2=B0C)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit seems to be the most common choice for room temperature --- PRIVATE | 2 +- src/constants.f90 | 2 +- src/rotations.f90 | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/PRIVATE b/PRIVATE index ebb7f0ce7..5774122bf 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit ebb7f0ce78d11275020af0ba60f929f95b446932 +Subproject commit 5774122bf48d637704bb4afb10b87c34a4dbcaba diff --git a/src/constants.f90 b/src/constants.f90 index 43098188b..1effb8d95 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -9,7 +9,7 @@ module constants public real(pReal), parameter :: & - T_ROOM = 298.15_pReal, & !< Room temperature in K (25°C)/Standard Ambient Temperaure and Pressure (SATP) + T_ROOM = 293.15_pReal, & !< Room temperature in K (20°C) K_B = 1.380649e-23_pReal, & !< Boltzmann constant in J/Kelvin (https://doi.org/10.1351/goldbook) N_A = 6.02214076e23_pReal !< Avogadro constant in 1/mol (https://doi.org/10.1351/goldbook) diff --git a/src/rotations.f90 b/src/rotations.f90 index 27bf8c227..53f1fe976 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -1338,16 +1338,16 @@ pure function cu2ho(cu) result(ho) ! transform to sphere grid (inverse Lambert) ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] c = sum(T**2) - s = Pi * c/(24.0*XYZ(3)**2) + s = PI * c/(24.0*XYZ(3)**2) c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3) q = sqrt( 1.0 - s ) LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ] - endif special + end if special ! reverse the coordinates back to order according to the original pyramid number ho = LamXYZ(p(:,2)) - endif center + end if center end function cu2ho From 85ca3a3f24b159993ff3c2e0a0f2a8ab1a7c4ef9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 3 Feb 2022 22:57:32 +0100 Subject: [PATCH 16/26] simplified/standardized --- src/phase.f90 | 12 ++++++++---- src/phase_mechanical.f90 | 2 +- src/phase_mechanical_plastic.f90 | 8 ++++---- src/phase_mechanical_plastic_nonlocal.f90 | 16 ++++++---------- 4 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/phase.f90 b/src/phase.f90 index bea1a5543..b7bace81a 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -408,7 +408,7 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - NEntries,sizeState,sizeDotState,sizeDeltaState) + NEntries,sizeState,sizeDotState,sizeDeltaState,offsetDeltaState) class(tState), intent(inout) :: & state @@ -417,12 +417,17 @@ subroutine phase_allocateState(state, & sizeState, & sizeDotState, & sizeDeltaState - + integer, intent(in), optional :: & + offsetDeltaState state%sizeState = sizeState state%sizeDotState = sizeDotState state%sizeDeltaState = sizeDeltaState - state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + if (present(offsetDeltaState)) then + state%offsetDeltaState = offsetDeltaState ! ToDo: this is a fix for broken nonlocal + else + state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + end if allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%state0 (sizeState,NEntries), source=0.0_pReal) @@ -432,7 +437,6 @@ subroutine phase_allocateState(state, & allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal) - end subroutine phase_allocateState diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 234bee8e9..04cc4b946 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -88,7 +88,7 @@ submodule(phase) mechanical en real(pReal), intent(in) :: & subdt !< timestep - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState end function plastic_dotState diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 6d2aac827..912aadc03 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -154,11 +154,10 @@ submodule(phase:mechanical) plastic en end subroutine dislotungsten_dotState - module subroutine nonlocal_dotState(Mp,Temperature,timestep,ph,en,ip,el) + module subroutine nonlocal_dotState(Mp,timestep,ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & @@ -308,7 +307,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) subdt !< timestep real(pReal), dimension(3,3) :: & Mp - real(pReal), dimension(plasticState(material_phaseID(co,(el-1)*discretization_nIPs+ip))%sizeDotState) :: & + real(pReal), dimension(plasticState(ph)%sizeDotState) :: & dotState @@ -334,7 +333,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(dotState) call dislotungsten_dotState(Mp,thermal_T(ph,en),ph,en) case (PLASTIC_NONLOCAL_ID) plasticType - call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el) + call nonlocal_dotState(Mp,subdt,ph,en,ip,el) end select plasticType end if @@ -394,6 +393,7 @@ module function plastic_deltaState(ph, en) result(broken) myOffset, & mySize + broken = .false. select case (phase_plasticity(ph)) diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 0dc3dc4bd..a8f56b5c5 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -403,7 +403,7 @@ module function plastic_nonlocal_init() result(myPlasticity) 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState) + call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0) ! ToDo: state structure does not follow convention allocate(geom(ph)%V_0(Nmembers)) call storeGeometry(ph) @@ -411,9 +411,6 @@ module function plastic_nonlocal_init() result(myPlasticity) if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) & call IO_error(212,ext_msg='IPneighborhood does not exist') - - plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention - st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) @@ -941,13 +938,12 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -module subroutine nonlocal_dotState(Mp, Temperature,timestep, & +module subroutine nonlocal_dotState(Mp,timestep, & ph,en,ip,el) real(pReal), dimension(3,3), intent(in) :: & Mp !< MandelStress real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & ph, & @@ -984,7 +980,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & real(pReal) :: & D_SD, & mu, & - nu + nu, Temperature if (timestep <= 0.0_pReal) then plasticState(ph)%dotState = 0.0_pReal @@ -995,6 +991,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & mu = elastic_mu(ph,en) nu = elastic_nu(ph,en) + Temperature = thermal_T(ph,en) tau = 0.0_pReal dot_gamma = 0.0_pReal @@ -1195,7 +1192,6 @@ function rhoDotFlux(timestep,ph,en,ip,el) associate(prm => param(ph), & dst => dependentState(ph), & - dot => dotState(ph), & stt => state(ph)) ns = prm%sum_N_sl @@ -1206,7 +1202,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) rho0 = getRho0(ph,en) my_rhoSgl0 = rho0(:,sgl) - forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here + forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en) !ToDo: MD: I think we should use state0 here dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en) @@ -1217,7 +1213,7 @@ function rhoDotFlux(timestep,ph,en,ip,el) if (plasticState(ph)%nonlocal) then !*** check CFL (Courant-Friedrichs-Lewy) condition for flux - if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... + if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... .and. prm%C_CFL * abs(v0) * timestep & > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) #ifdef DEBUG From c25c17d0b6fe2d3eaef9fb6c431bc7049623cb0f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 10:57:18 +0100 Subject: [PATCH 17/26] reading the gcc manual -mtune is implied by -march, -ffast-math is not recommended, and -mtune=native is save --- cmake/Compiler-GNU.cmake | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index 5e2bf2ce6..a06db2c4e 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -12,9 +12,9 @@ endif () if (OPTIMIZATION STREQUAL "OFF") set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") - set (OPTIMIZATION_FLAGS "-O2 -mtune=generic -flto") + set (OPTIMIZATION_FLAGS "-O2 -mtune=native -flto") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") - set (OPTIMIZATION_FLAGS "-O3 -march=native -mtune=native -ffast-math -funroll-loops -ftree-vectorize -flto") + set (OPTIMIZATION_FLAGS "-O3 -march=native -funroll-loops -ftree-vectorize -flto") endif () set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) @@ -123,6 +123,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -ffpe-trap=invalid,zero,overflow") set (DEBUG_FLAGS "${DEBUG_FLAGS} -g") # Generate symbolic debugging information in the object file +set (DEBUG_FLAGS "${DEBUG_FLAGS} -Og") +# Optimize debugging experience + set (DEBUG_FLAGS "${DEBUG_FLAGS} -fbacktrace") set (DEBUG_FLAGS "${DEBUG_FLAGS} -fdump-core") set (DEBUG_FLAGS "${DEBUG_FLAGS} -fcheck=all") From f36db86b3c2e4e2f6737e4083a019e8a536ed2d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 18:24:24 +0100 Subject: [PATCH 18/26] better debug for gfortran, removing options not needed --- CMakeLists.txt | 2 +- cmake/Compiler-GNU.cmake | 12 +++++------- cmake/Compiler-Intel.cmake | 7 +++++-- cmake/Compiler-IntelLLVM.cmake | 5 ++++- 4 files changed, 15 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8c7b129f0..b4c405319 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,7 +42,7 @@ string(TOUPPER "${CMAKE_BUILD_TYPE}" CMAKE_BUILD_TYPE) if(CMAKE_BUILD_TYPE STREQUAL "DEBUG" OR CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") set(DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG") set(PARALLEL "OFF") - set(OPTI "OFF") + set(OPTI "DEBUG") elseif(CMAKE_BUILD_TYPE STREQUAL "RELEASE") set(PARALLEL "ON") set(OPTI "DEFENSIVE") diff --git a/cmake/Compiler-GNU.cmake b/cmake/Compiler-GNU.cmake index a06db2c4e..c204d1af2 100644 --- a/cmake/Compiler-GNU.cmake +++ b/cmake/Compiler-GNU.cmake @@ -9,26 +9,24 @@ if (OPENMP) set (OPENMP_FLAGS "-fopenmp") endif () -if (OPTIMIZATION STREQUAL "OFF") +if (OPTIMIZATION STREQUAL "DEBUG") + set (OPTIMIZATION_FLAGS "-Og") +elseif (OPTIMIZATION STREQUAL "OFF") set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") - set (OPTIMIZATION_FLAGS "-O2 -mtune=native -flto") + set (OPTIMIZATION_FLAGS "-O2 -mtune=native") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") set (OPTIMIZATION_FLAGS "-O3 -march=native -funroll-loops -ftree-vectorize -flto") endif () set (STANDARD_CHECK "-std=f2018 -pedantic-errors" ) -set (LINKER_FLAGS "${LINKER_FLAGS} -Wl") -# options parsed directly to the linker -set (LINKER_FLAGS "${LINKER_FLAGS},-undefined,dynamic_lookup" ) -# ensure to link against dynamic libraries #------------------------------------------------------------------------------------------------ # Fine tuning compilation options set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor -set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIC -fPIE") +set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE") # position independent code set (COMPILE_FLAGS "${COMPILE_FLAGS} -ffree-line-length-132") diff --git a/cmake/Compiler-Intel.cmake b/cmake/Compiler-Intel.cmake index 3afffd2be..7f34e4a13 100644 --- a/cmake/Compiler-Intel.cmake +++ b/cmake/Compiler-Intel.cmake @@ -9,12 +9,12 @@ if (OPENMP) set (OPENMP_FLAGS "-qopenmp -parallel") endif () -if (OPTIMIZATION STREQUAL "OFF") +if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG") set (OPTIMIZATION_FLAGS "-O0 -no-ip") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") set (OPTIMIZATION_FLAGS "-O2") elseif (OPTIMIZATION STREQUAL "AGGRESSIVE") - set (OPTIMIZATION_FLAGS "-ipo -O3 -no-prec-div -fp-model fast=2 -xHost") + set (OPTIMIZATION_FLAGS "-ipo -O3 -fp-model fast=2 -xHost") # -fast = -ipo, -O3, -no-prec-div, -static, -fp-model fast=2, and -xHost" endif () @@ -110,6 +110,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") # generate debug information for parameters # Disabled due to ICE when compiling phase_damage.f90 (not understandable, there is no parameter in there) +set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all") +# generate complete debugging information + # Additional options # -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits # -check: Checks at runtime, where diff --git a/cmake/Compiler-IntelLLVM.cmake b/cmake/Compiler-IntelLLVM.cmake index 326cfe319..883873e1c 100644 --- a/cmake/Compiler-IntelLLVM.cmake +++ b/cmake/Compiler-IntelLLVM.cmake @@ -9,7 +9,7 @@ if (OPENMP) set (OPENMP_FLAGS "-qopenmp") endif () -if (OPTIMIZATION STREQUAL "OFF") +if (OPTIMIZATION STREQUAL "OFF" OR OPTIMIZATION STREQUAL "DEBUG") set (OPTIMIZATION_FLAGS "-O0") elseif (OPTIMIZATION STREQUAL "DEFENSIVE") set (OPTIMIZATION_FLAGS "-O2") @@ -109,6 +109,9 @@ set (DEBUG_FLAGS "${DEBUG_FLAGS} -fpe-all=0") set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug-parameters all") # generate debug information for parameters +set (DEBUG_FLAGS "${DEBUG_FLAGS} -debug all") +# generate complete debugging information + # Additional options # -heap-arrays: Should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits # -check: Checks at runtime, where From 12e7922fafa6a2b818d437456e14fbf591d6d112 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 18:32:17 +0100 Subject: [PATCH 19/26] use modern Fortran interface not possible for HDF5... --- src/HDF5_utilities.f90 | 33 ++++++++++--------- ...hase_mechanical_eigen_thermalexpansion.f90 | 2 +- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index d2076c2cc..9dc53ee17 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -10,8 +10,9 @@ module HDF5_utilities #include use PETScSys #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) - use MPI + use MPI_f08 #endif + use MPI, only: MPI_INFO_NULL_F90 => MPI_INFO_NULL #endif use prec @@ -162,9 +163,9 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel) character, intent(in), optional :: mode logical, intent(in), optional :: parallel - character :: m - integer(HID_T) :: plist_id - integer :: hdferr + character :: m + integer(HID_T) :: plist_id + integer :: hdferr if (present(mode)) then @@ -178,9 +179,9 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel) #ifdef PETSC if (present(parallel)) then - if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr) else - call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr) end if if(hdferr < 0) error stop 'HDF5 error' #endif @@ -1860,7 +1861,7 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_ globalShape !< shape of the dataset (all processes) integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id, aplist_id - integer, dimension(worldsize) :: & + integer(MPI_INTEGER_KIND), dimension(worldsize) :: & readSize !< contribution of all processes integer :: hdferr integer(MPI_INTEGER_KIND) :: err_MPI @@ -1871,13 +1872,13 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_ if(hdferr < 0) error stop 'HDF5 error' !-------------------------------------------------------------------------------------------------- - readSize = 0 - readSize(worldrank+1) = int(localShape(ubound(localShape,1))) + readSize = 0_MPI_INTEGER_KIND + readSize(worldrank+1) = int(localShape(ubound(localShape,1)),MPI_INTEGER_KIND) #ifdef PETSC if (parallel) then call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr) if(hdferr < 0) error stop 'HDF5 error' - call MPI_allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process + call MPI_Allreduce(MPI_IN_PLACE,readSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' end if #endif @@ -1954,8 +1955,8 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & totalShape !< shape of the dataset (all processes) integer(HID_T), intent(out) :: dset_id, filespace_id, memspace_id, plist_id - integer, dimension(worldsize) :: writeSize !< contribution of all processes - integer(HID_T) :: dcpl + integer(MPI_INTEGER_KIND), dimension(worldsize) :: writeSize !< contribution of all processes + integer(HID_T) :: dcpl integer :: hdferr integer(MPI_INTEGER_KIND) :: err_MPI integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T @@ -1974,11 +1975,11 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & !-------------------------------------------------------------------------------------------------- ! determine the global data layout among all processes - writeSize = 0 - writeSize(worldrank+1) = int(myShape(ubound(myShape,1))) + writeSize = 0_MPI_INTEGER_KIND + writeSize(worldrank+1) = int(myShape(ubound(myShape,1)),MPI_INTEGER_KIND) #ifdef PETSC if (parallel) then - call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,err_MPI) ! get total output size over each process + call MPI_Allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI) ! get total output size over each process if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' end if #endif @@ -2009,7 +2010,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, & if (hdferr < 0) error stop 'HDF5 error' end if end if - + !-------------------------------------------------------------------------------------------------- ! create dataspace in memory (local shape) and in file (global shape) call H5Screate_simple_f(size(myShape), myShape, memspace_id, hdferr, myShape) diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index a5d9868a8..3c422616b 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -27,7 +27,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length logical, dimension(:,:), allocatable :: myKinematics - integer :: Ninstances,p,i,k + integer :: Ninstances, p, k class(tNode), pointer :: & phases, & phase, & From 4ca0ea6af238bd1aaef9339d8039fdf06d71328b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 18:32:53 +0100 Subject: [PATCH 20/26] avoid linking issues with gfortran+MPI most likely related to the fact that HDF5 uses the old Fortran inferface, not MPI_f08 as DAMASK --- src/grid/grid_mech_spectral_basic.f90 | 8 +++++++- src/grid/grid_mech_spectral_polarisation.f90 | 8 +++++++- 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index fa2e17bd9..2f2b73f01 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -79,6 +79,12 @@ module grid_mechanical_spectral_basic err_BC, & !< deviation from stress BC err_div !< RMS of div of P +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + type(MPI_Status) :: status +#else + integer, dimension(MPI_STATUS_SIZE) :: status +#endif + integer :: & totalIter = 0 !< total iteration in current increment @@ -244,7 +250,7 @@ subroutine grid_mechanical_spectral_basic_init call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', & MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI) + call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_File_close(fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index ee87c77a6..b72cc4232 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -90,6 +90,12 @@ module grid_mechanical_spectral_polarisation err_curl, & !< RMS of curl of F err_div !< RMS of div of P +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) + type(MPI_Status) :: status +#else + integer, dimension(MPI_STATUS_SIZE) :: status +#endif + integer :: & totalIter = 0 !< total iteration in current increment @@ -270,7 +276,7 @@ subroutine grid_mechanical_spectral_polarisation_init call MPI_File_open(MPI_COMM_WORLD, trim(getSolverJobName())//'.C_ref', & MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_STATUS_IGNORE,err_MPI) + call MPI_File_read(fileUnit,C_minMaxAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,status,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' call MPI_File_close(fileUnit,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' From 519ca6f99044a566a671b2b732644c16ad9852ea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Feb 2022 20:00:55 +0100 Subject: [PATCH 21/26] ifort does not allow rename of already imported symbols --- src/HDF5_utilities.f90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90 index 9dc53ee17..93de4053b 100644 --- a/src/HDF5_utilities.f90 +++ b/src/HDF5_utilities.f90 @@ -11,8 +11,8 @@ module HDF5_utilities use PETScSys #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) use MPI_f08 -#endif use MPI, only: MPI_INFO_NULL_F90 => MPI_INFO_NULL +#endif #endif use prec @@ -179,9 +179,15 @@ integer(HID_T) function HDF5_openFile(fileName,mode,parallel) #ifdef PETSC if (present(parallel)) then +#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr) else call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL_F90, hdferr) +#else + if (parallel) call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) + else + call H5Pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) +#endif end if if(hdferr < 0) error stop 'HDF5 error' #endif From 0cf01120295e3de772196cfd849eb7915d1d1979 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 15:20:26 +0100 Subject: [PATCH 22/26] simplified --- src/system_routines.f90 | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/src/system_routines.f90 b/src/system_routines.f90 index 7a4e41a57..2eb0b7958 100644 --- a/src/system_routines.f90 +++ b/src/system_routines.f90 @@ -24,7 +24,7 @@ module system_routines function setCWD_C(cwd) bind(C) use, intrinsic :: ISO_C_Binding, only: C_INT, C_CHAR - + integer(C_INT) :: setCWD_C character(kind=C_CHAR), dimension(*), intent(in) :: cwd end function setCWD_C @@ -150,14 +150,14 @@ function getUserName() getUserName = c_f_string(getUserName_Cstring) else getUserName = 'n/a (Error!)' - endif + end if end function getUserName !-------------------------------------------------------------------------------------------------- -!> @brief convert C string to Fortran string -!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string +!> @brief Convert C string to Fortran string. +!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string. !-------------------------------------------------------------------------------------------------- pure function c_f_string(c_string) result(f_string) @@ -174,28 +174,23 @@ pure function c_f_string(c_string) result(f_string) else f_string = f_string(:i-1) exit - endif - enddo arrayToString + end if + end do arrayToString end function c_f_string !-------------------------------------------------------------------------------------------------- -!> @brief convert Fortran string to C string -!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string +!> @brief Convert Fortran string to C string. +!> @details: C string is NULL terminated and, hence, longer by one than the Fortran string. !-------------------------------------------------------------------------------------------------- pure function f_c_string(f_string) result(c_string) character(len=*), intent(in) :: f_string character(kind=C_CHAR), dimension(len_trim(f_string)+1) :: c_string - integer :: i - - do i=1,len_trim(f_string) - c_string(i)=f_string(i:i) - enddo - c_string(len_trim(f_string)+1) = C_NULL_CHAR + c_string = transfer(trim(f_string)//C_NULL_CHAR,c_string,size=size(c_string)) end function f_c_string From a37438ca291dc28d8323b06368f4aba514f07823 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 22:07:13 +0100 Subject: [PATCH 23/26] helpful information about the parallel environment --- src/parallelization.f90 | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/parallelization.f90 b/src/parallelization.f90 index 29deaf724..28ad70d94 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -52,13 +52,13 @@ contains !-------------------------------------------------------------------------------------------------- subroutine parallelization_init - integer(MPI_INTEGER_KIND) :: err_MPI, typeSize + integer(MPI_INTEGER_KIND) :: err_MPI, typeSize, version, subversion, devNull character(len=4) :: rank_str + character(len=MPI_MAX_LIBRARY_VERSION_STRING) :: MPI_library_version !$ integer :: got_env, threadLevel !$ integer(pI32) :: OMP_NUM_THREADS !$ character(len=6) NumThreadsString - PetscErrorCode :: err_PETSc #ifdef _OPENMP ! If openMP is enabled, check if the MPI libary supports it and initialize accordingly. @@ -86,12 +86,22 @@ subroutine parallelization_init if (err_MPI /= 0_MPI_INTEGER_KIND) & error stop 'Could not determine worldrank' - if (worldrank == 0) print'(/,1x,a)', '<<<+- parallelization init -+>>>' + if (worldrank == 0) then + print'(/,1x,a)', '<<<+- parallelization init -+>>>' + + call MPI_Get_library_version(MPI_library_version,devNull,err_MPI) + print'(/,1x,a)', trim(MPI_library_version) + call MPI_Get_version(version,subversion,err_MPI) + print'(1x,a,i0,a,i0)', 'MPI standard: ',version,'.',subversion +#ifdef _OPENMP + print'(1x,a,i0)', 'OpenMP version: ',openmp_version +#endif + end if call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) & error stop 'Could not determine worldsize' - if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize + if (worldrank == 0) print'(/,1x,a,i0)', 'MPI processes: ',worldsize call MPI_Type_size(MPI_INTEGER,typeSize,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) & @@ -128,7 +138,7 @@ subroutine parallelization_init !$ OMP_NUM_THREADS = 4_pI32 !$ endif !$ endif -!$ print'(1x,a,1x,i2)', 'OMP_NUM_THREADS:',OMP_NUM_THREADS +!$ print'(1x,a,i0)', 'OMP_NUM_THREADS: ',OMP_NUM_THREADS !$ call omp_set_num_threads(OMP_NUM_THREADS) end subroutine parallelization_init From 2a59251b1d8a500865738e59f3f1a683ef27d2ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 6 Feb 2022 22:23:21 +0100 Subject: [PATCH 24/26] deprecated reporting crystallite does not exist anymore, information was discretization related --- src/CPFEM.f90 | 1 - src/CPFEM2.f90 | 1 - src/phase.f90 | 24 +++++------------------- 3 files changed, 5 insertions(+), 21 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 0f7590782..c8b7e3cca 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -90,7 +90,6 @@ subroutine CPFEM_initAll call material_init(.false.) call phase_init call homogenization_init - call crystallite_init call CPFEM_init call config_deallocate diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index ed8fb611b..b24ba5480 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -68,7 +68,6 @@ subroutine CPFEM_initAll call material_init(restart=interface_restartInc>0) call phase_init call homogenization_init - call crystallite_init call CPFEM_init call config_deallocate diff --git a/src/phase.f90 b/src/phase.f90 index b7bace81a..e271e6366 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -323,7 +323,6 @@ module phase phase_restore, & plastic_nonlocal_updateCompatibility, & converged, & - crystallite_init, & phase_mechanical_constitutive, & phase_thermal_constitutive, & phase_damage_constitutive, & @@ -401,6 +400,8 @@ subroutine phase_init call damage_init call thermal_init(phases) + call crystallite_init() + end subroutine phase_init @@ -500,22 +501,13 @@ subroutine crystallite_init() ce, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop - el, & !< counter in element loop - cMax, & !< maximum number of integration point components - iMax, & !< maximum number of integration points - eMax !< maximum number of elements + el !< counter in element loop class(tNode), pointer :: & num_crystallite, & phases - print'(/,1x,a)', '<<<+- crystallite init -+>>>' - - cMax = homogenization_maxNconstituents - iMax = discretization_nIPs - eMax = discretization_Nelems - num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal) @@ -549,15 +541,9 @@ subroutine crystallite_init() phases => config_material%get('phase') - print'(/,a42,1x,i10)', ' # of elements: ', eMax - print'( a42,1x,i10)', ' # of integration points/element: ', iMax - print'( a42,1x,i10)', 'max # of constituents/integration point: ', cMax - flush(IO_STDOUT) - - !$OMP PARALLEL DO PRIVATE(ce) - do el = 1, eMax - do ip = 1, iMax + do el = 1, discretization_Nelems + do ip = 1, discretization_nIPs ce = (el-1)*discretization_nIPs + ip do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) call crystallite_orientations(co,ip,el) From e3a9adc7222a145c9b2fd1affa3f4bcb3e7f9e66 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 7 Feb 2022 07:20:15 +0100 Subject: [PATCH 25/26] not used --- src/IO.f90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index ec7cf9704..441a0bf3d 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -584,8 +584,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) character(len=pStringLen) :: formatString select case (warning_ID) - case (42) - msg = 'parameter has no effect' case (47) msg = 'no valid parameter for FFTW, using FFTW_PATIENT' case (207) @@ -594,8 +592,6 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) msg = 'crystallite responds elastically' case (601) msg = 'stiffness close to zero' - case (700) - msg = 'unknown crystal symmetry' case (709) msg = 'read only the first document' case default From bd657e6c62145f2fa47ee7aa83ae292f5bc511db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 9 Feb 2022 12:29:55 +0100 Subject: [PATCH 26/26] fixed style and corrected copy and paste errors --- src/grid/DAMASK_grid.f90 | 2 +- src/grid/grid_mech_FEM.f90 | 2 +- src/math.f90 | 2 +- src/phase_mechanical.f90 | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index f1834d16b..918f477b1 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -32,7 +32,7 @@ program DAMASK_grid implicit none type :: tLoadCase - type(tRotation) :: rot !< rotation of BC + type(tRotation) :: rot !< rotation of BC type(tBoundaryCondition) :: stress, & !< stress BC deformation !< deformation BC (dot_F, F, or L) real(pReal) :: t, & !< length of increment diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 6ad035fc2..3c42c8c23 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -339,7 +339,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(tRotation), intent(in) :: & + type(tRotation), intent(in) :: & rotation_BC PetscErrorCode :: err_PETSc PetscScalar, pointer, dimension(:,:,:,:) :: & diff --git a/src/math.f90 b/src/math.f90 index 5c6b10fdd..dd4690672 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -898,7 +898,7 @@ end function math_33toVoigt6_stress !-------------------------------------------------------------------------------------------------- -!> @brief Convert 3x3 stress tensor into 6 Voigt vector. +!> @brief Convert 3x3 strain tensor into 6 Voigt vector. !-------------------------------------------------------------------------------------------------- pure function math_33toVoigt6_strain(epsilon) result(epsilon_tilde) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 04cc4b946..9a7d1ba62 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -879,10 +879,10 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D plasticState(ph)%state(1:sizeDotState,en) = subState0 & + dotState * Delta_t - broken = integrateStress(F_0 + (F - F_0) * Delta_t * C(stage),subFp0,subFi0,Delta_t * C(stage),co,ip,el) + broken = integrateStress(F_0 + (F-F_0) * Delta_t*C(stage),subFp0,subFi0,Delta_t*C(stage),co,ip,el) if(broken) exit - dotState = plastic_dotState(Delta_t, co,ip,el,ph,en) + dotState = plastic_dotState(Delta_t*C(stage), co,ip,el,ph,en) if (any(IEEE_is_NaN(dotState))) exit enddo