diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 9556e2ba1..27c0d04ce 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -109,7 +109,8 @@ end subroutine !********************************************************* subroutine CPFEM_init() - + + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt use debug, only: debug_verbosity use IO, only: IO_read_jobBinaryFile @@ -128,9 +129,7 @@ subroutine CPFEM_init() crystallite_dPdF0, & crystallite_Tstar0_v use homogenization, only: homogenization_sizeState, & - homogenization_state0, & - materialpoint_F, & - materialpoint_F0 + homogenization_state0 implicit none @@ -231,12 +230,10 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, !*** variables and functions from other modules ***! use prec, only: pReal, & pInt - use numerics, only: relevantStrain, & - defgradTolerance, & + use numerics, only: defgradTolerance, & iJacoStiffness use debug, only: debug_e, & debug_i, & - debug_g, & debug_selectiveDebugger, & debug_verbosity, & debug_stressMaxLocation, & @@ -282,8 +279,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, microstructure_elemhomo, & material_phase use constitutive, only: constitutive_state0,constitutive_state - use crystallite, only: crystallite_F0, & - crystallite_partionedF, & + use crystallite, only: crystallite_partionedF,& + crystallite_F0, & crystallite_Fp0, & crystallite_Fp, & crystallite_Lp0, & @@ -291,8 +288,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP, crystallite_dPdF0, & crystallite_dPdF, & crystallite_Tstar0_v, & - crystallite_Tstar_v, & - crystallite_localConstitution + crystallite_Tstar_v use homogenization, only: homogenization_sizeState, & homogenization_state, & homogenization_state0, & diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 4956f349a..8a71ee54e 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -103,14 +103,12 @@ CONTAINS subroutine crystallite_init(Temperature) !*** variables and functions from other modules ***! -use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) +use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt, & pReal use debug, only: debug_info, & debug_reset, & debug_verbosity -use numerics, only: subStepSizeCryst, & - stepIncreaseCryst use math, only: math_I3, & math_EulerToR, & math_inv33, & @@ -125,14 +123,11 @@ use mesh, only: mesh_element, & mesh_maxNipNeighbors use IO use material -use lattice, only: lattice_symmetryType, & - lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem +use lattice, only: lattice_symmetryType use constitutive, only: constitutive_microstructure use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, & - constitutive_phenopowerlaw_structure, & - constitutive_phenopowerlaw_Nslip + constitutive_phenopowerlaw_structure use constitutive_titanmod, only: constitutive_titanmod_label, & constitutive_titanmod_structure use constitutive_dislotwin, only: constitutive_dislotwin_label, & @@ -141,8 +136,8 @@ use constitutive_nonlocal, only: constitutive_nonlocal_label, & constitutive_nonlocal_structure implicit none -integer(pInt), parameter :: file = 200, & - maxNchunks = 2 +integer(pInt), parameter :: myFile = 200_pInt, & + maxNchunks = 2_pInt !*** input variables ***! real(pReal) Temperature @@ -234,18 +229,18 @@ allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt -if (.not. IO_open_jobFile_stat(file,material_localFileExt)) then ! no local material configuration present... - call IO_open_file(file,material_configFile) ! ...open material.config file +if (.not. IO_open_jobFile_stat(myFile,material_localFileExt)) then ! no local material configuration present... + call IO_open_file(myFile,material_configFile) ! ...open material.config file endif line = '' section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do ! read thru sections of phase part - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section @@ -263,7 +258,7 @@ do ! read thru sections of endif enddo -100 close(file) +100 close(myFile) do i = 1_pInt,material_Ncrystallite ! sanity checks enddo @@ -299,18 +294,18 @@ enddo ! write description file for crystallite output -call IO_write_jobFile(file,'outputCrystallite') +call IO_write_jobFile(myFile,'outputCrystallite') do p = 1_pInt,material_Ncrystallite - write(file,*) - write(file,'(a)') '['//trim(crystallite_name(p))//']' - write(file,*) - do e = 1,crystallite_Noutput(p) - write(file,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p) + write(myFile,*) + write(myFile,'(a)') '['//trim(crystallite_name(p))//']' + write(myFile,*) + do e = 1_pInt,crystallite_Noutput(p) + write(myFile,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p) enddo enddo -close(file) +close(myFile) !$OMP PARALLEL PRIVATE(myNgrains,myPhase,myMat,myStructure) @@ -318,7 +313,7 @@ close(file) do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element - do g = 1,myNgrains + do g = 1_pInt,myNgrains crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation crystallite_F0(1:3,1:3,g,i,e) = math_I3 crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) @@ -342,7 +337,7 @@ crystallite_requested = .true. do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,myNgrains + do g = 1_pInt,myNgrains myPhase = material_phase(g,i,e) myMat = phase_constitutionInstance(myPhase) select case (phase_constitution(myPhase)) @@ -374,7 +369,7 @@ crystallite_orientation0 = crystallite_orientation ! Store initial o do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,myNgrains + do g = 1_pInt,myNgrains call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states enddo @@ -459,7 +454,6 @@ use numerics, only: subStepMinCryst, & pert_Fg, & pert_method, & nCryst, & - iJacoStiffness, & numerics_integrator, & numerics_integrationMode use debug, only: debug_verbosity, & @@ -557,7 +551,7 @@ crystallite_subStep = 0.0_pReal do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains + do g = 1_pInt,myNgrains if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure @@ -710,7 +704,7 @@ enddo .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then write (6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g write (6,*) - write (6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e)) / 1e6 + write (6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)) write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)) write (6,*) @@ -922,10 +916,10 @@ use prec, only: pInt, & pReal use numerics, only: numerics_integrationMode use debug, only: debug_verbosity, & + debug_selectiveDebugger, & debug_e, & debug_i, & debug_g, & - debug_selectiveDebugger, & debug_StateLoopDistribution use FEsolving, only: FEsolving_execElem, & FEsolving_execIP @@ -1029,7 +1023,7 @@ RK4dotTemperature = 0.0_pReal ! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION --- -do n = 1,4 +do n = 1_pInt,4_pInt ! --- state update --- @@ -1177,19 +1171,16 @@ subroutine crystallite_integrateStateRKCK45(gg,ii,ee) use prec, only: pInt, & pReal use debug, only: debug_verbosity, & + debug_selectiveDebugger, & debug_e, & debug_i, & debug_g, & - debug_selectiveDebugger, & debug_StateLoopDistribution use numerics, only: rTol_crystalliteState, & rTol_crystalliteTemperature, & - subStepSizeCryst, & - stepIncreaseCryst, & numerics_integrationMode use FEsolving, only: FEsolving_execElem, & - FEsolving_execIP, & - theInc + FEsolving_execIP use mesh, only: mesh_element, & mesh_NcpElems, & mesh_maxNips @@ -1292,7 +1283,7 @@ else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) iIter(1:2,e) = FEsolving_execIP(1:2,e) - gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))] enddo singleRun = .false. endif @@ -1345,7 +1336,7 @@ endif ! --- SECOND TO SIXTH RUNGE KUTTA STEP --- -do n = 1,5 +do n = 1_pInt,5_pInt ! --- state update --- @@ -1558,7 +1549,7 @@ relTemperatureResiduum = 0.0_pReal do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) - forall (s = 1:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & + forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) if (crystallite_Temperature(g,i,e) > 0) & relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) @@ -1571,7 +1562,7 @@ relTemperatureResiduum = 0.0_pReal .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) #ifndef _OPENMP - if (debug_verbosity > 5 & + if (debug_verbosity > 5_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then write(6,'(a,i8,1x,i3,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g write(6,*) @@ -1614,7 +1605,8 @@ relTemperatureResiduum = 0.0_pReal crystallite_todo(g,i,e) = .false. ! ... integration done if (debug_verbosity > 4) then !$OMP CRITICAL (distributionState) - debug_StateLoopDistribution(6,numerics_integrationMode) = debug_StateLoopDistribution(6,numerics_integrationMode) + 1 + debug_StateLoopDistribution(6,numerics_integrationMode) =& + debug_StateLoopDistribution(6,numerics_integrationMode) + 1_pInt !$OMP END CRITICAL (distributionState) endif else @@ -1666,8 +1658,6 @@ use debug, only: debug_verbosity, & debug_StateLoopDistribution use numerics, only: rTol_crystalliteState, & rTol_crystalliteTemperature, & - subStepSizeCryst, & - stepIncreaseCryst, & numerics_integrationMode use FEsolving, only: FEsolving_execElem, & FEsolving_execIP @@ -1876,9 +1866,9 @@ relTemperatureResiduum = 0.0_pReal ! --- relative residui --- - forall (s = 1:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & + forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) - if (crystallite_Temperature(g,i,e) > 0) & + if (crystallite_Temperature(g,i,e) > 0_pInt) & relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) @@ -2136,10 +2126,6 @@ subroutine crystallite_integrateStateFPI(gg,ii,ee) use prec, only: pInt, & pReal use debug, only: debug_verbosity, & - debug_selectiveDebugger, & - debug_e, & - debug_i, & - debug_g, & debug_StateLoopDistribution use numerics, only: nState, & numerics_integrationMode @@ -2148,9 +2134,7 @@ use FEsolving, only: FEsolving_execElem, & use mesh, only: mesh_element, & mesh_NcpElems use material, only: homogenization_Ngrains -use constitutive, only: constitutive_sizeDotState, & - constitutive_state, & - constitutive_dotState, & +use constitutive, only: constitutive_dotState, & constitutive_collectDotState, & constitutive_dotTemperature, & constitutive_microstructure, & @@ -2192,7 +2176,7 @@ else eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2) iIter(1:2,e) = FEsolving_execIP(1:2,e) - gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/) + gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))] enddo singleRun = .false. endif @@ -2346,7 +2330,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState ) if (debug_verbosity > 4) then !$OMP CRITICAL (distributionState) debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = & - debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1 + debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1_pInt !$OMP END CRITICAL (distributionState) endif endif @@ -2411,11 +2395,8 @@ endsubroutine subroutine crystallite_updateState(done, converged, g, i, e) !*** variables and functions from other modules ***! -use prec, only: pReal, & - pInt, & - pLongInt -use numerics, only: rTol_crystalliteState, & - numerics_integrationMode +use prec, only: pInt +use numerics, only: rTol_crystalliteState use constitutive, only: constitutive_dotState, & constitutive_previousDotState, & constitutive_sizeDotState, & @@ -2424,10 +2405,10 @@ use constitutive, only: constitutive_dotState, & constitutive_aTolState, & constitutive_microstructure use debug, only: debug_verbosity, & - debug_g, & - debug_i, & + debug_selectiveDebugger, & debug_e, & - debug_selectiveDebugger + debug_i, & + debug_g !*** input variables ***! integer(pInt), intent(in):: e, & ! element index @@ -2513,13 +2494,10 @@ endsubroutine subroutine crystallite_updateTemperature(done, converged, g, i, e) !*** variables and functions from other modules ***! -use prec, only: pReal, & - pInt, & - pLongInt +use prec, only: pInt use numerics, only: rTol_crystalliteTemperature use constitutive, only: constitutive_dotTemperature -use debug, only: debug_verbosity - +use debug, only: debug_verbosity !*** input variables ***! integer(pInt), intent(in):: e, & ! element index i, & ! integration point index @@ -2591,17 +2569,16 @@ use numerics, only: nStress, & relevantStrain, & numerics_integrationMode use debug, only: debug_verbosity, & - debug_g, & - debug_i, & - debug_e, & debug_selectiveDebugger, & + debug_e, & + debug_i, & + debug_g, & debug_cumLpCalls, & debug_cumLpTicks, & debug_StressLoopDistribution, & debug_LeapfrogBreakDistribution -use constitutive, only: constitutive_homogenizedC, & - constitutive_LpAndItsTangent, & - constitutive_state +use constitutive, only: constitutive_LpAndItsTangent, & + constitutive_homogenizedC use math, only: math_mul33x33, & math_mul33xx33, & math_mul66x6, & @@ -2739,14 +2716,14 @@ steplength_max = 16.0_pReal jacoCounter = 0_pInt LpLoop: do - NiterationStress = NiterationStress + 1 + NiterationStress = NiterationStress + 1_pInt !* too many loops required ? if (NiterationStress > nStress) then #ifndef _OPENMP - if (debug_verbosity > 4) then + if (debug_verbosity > 4_pInt) then write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g write(6,*) endif @@ -2764,7 +2741,7 @@ LpLoop: do Tstar_v = 0.5_pReal * math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB) - math_I3)) p_hydro = sum(Tstar_v(1:3)) / 3.0_pReal - forall(n=1:3) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor + forall(n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor !* calculate plastic velocity gradient and its tangent according to constitutive law @@ -2875,7 +2852,7 @@ LpLoop: do if (debug_verbosity > 4) then !$OMP CRITICAL (distributionLeapfrogBreak) debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) = & - debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1 + debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1_pInt !$OMP END CRITICAL (distributionLeapfrogBreak) endif endif @@ -2887,7 +2864,7 @@ LpLoop: do if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then dT_dLp = 0.0_pReal - do h=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3 + do h=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do m=1_pInt,3_pInt dT_dLp(3*(h-1)+j,3*(k-1)+l) = dT_dLp(3*(h-1)+j,3*(k-1)+l) + C(h,j,l,m) * AB(k,m) + C(h,j,m,l) * BTA(m,k) enddo; enddo; enddo; enddo; enddo dT_dLp = -0.5_pReal * dt * dT_dLp @@ -2914,13 +2891,13 @@ LpLoop: do return endif deltaLp = 0.0_pReal - do k=1,3; do l=1,3; do m=1,3; do n=1,3 - deltaLp(k,l) = deltaLp(k,l) - inv_dR_dLp(3*(k-1)+l,3*(m-1)+n) * residuum(m,n) + do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt + deltaLp(k,l) = deltaLp(k,l) - inv_dR_dLp(3_pInt*(k-1_pInt)+l,3_pInt*(m-1_pInt)+n) * residuum(m,n) enddo; enddo; enddo; enddo gradientR = 0.0_pReal - do k=1,3; do l=1,3; do m=1,3; do n=1,3 - gradientR(k,l) = gradientR(k,l) + dR_dLp(3*(k-1)+l,3*(m-1)+n) * residuum(m,n) + do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt + gradientR(k,l) = gradientR(k,l) + dR_dLp(3*(k-1)+l,3_pInt*(m-1_pInt)+n) * residuum(m,n) enddo; enddo; enddo; enddo gradientR = gradientR / math_norm33(gradientR) expectedImprovement = math_mul33xx33(deltaLp, gradientR) @@ -2942,8 +2919,8 @@ call math_invert33(invFp_new,Fp_new,det,error) if (error) then #ifndef _OPENMP if (debug_verbosity > 4) then - write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, & - ' ; iteration ', NiterationStress + write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',& + e,i,g, ' ; iteration ', NiterationStress if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then write(6,*) write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new) @@ -2957,7 +2934,7 @@ Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resu !* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress -forall (n=1:3) Tstar_v(n) = Tstar_v(n) + p_hydro +forall (n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) + p_hydro crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose33(invFp_new))) @@ -2976,11 +2953,11 @@ crystallite_integrateStress = .true. #ifndef _OPENMP if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) & .and. numerics_integrationMode == 1_pInt) then - write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1e6 + write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', & - math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1e6 / math_det33(Fg_new) + math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1.0e6_pReal / math_det33(Fg_new) write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', & - math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order + math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)) endif #endif @@ -2988,7 +2965,7 @@ endif if (debug_verbosity > 4) then !$OMP CRITICAL (distributionStress) debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = & - debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1 + debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1_pInt !$OMP END CRITICAL (distributionStress) endif @@ -3007,22 +2984,17 @@ use prec, only: pInt, & use math, only: math_pDecomposition, & math_RtoQuaternion, & math_QuaternionDisorientation, & - inDeg, & math_qConj use FEsolving, only: FEsolving_execElem, & FEsolving_execIP use IO, only: IO_warning use material, only: material_phase, & homogenization_Ngrains, & - phase_constitution, & phase_localConstitution, & phase_constitutionInstance use mesh, only: mesh_element, & mesh_ipNeighborhood, & FE_NipNeighbors -use debug, only: debug_verbosity, & - debug_selectiveDebugger, & - debug_e, debug_i, debug_g use constitutive_nonlocal, only: constitutive_nonlocal_structure, & constitutive_nonlocal_updateCompatibility @@ -3054,7 +3026,7 @@ logical error !$OMP PARALLEL DO PRIVATE(error,U,R,orientation) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,homogenization_Ngrains(mesh_element(3,e)) + do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error) ! polar decomposition of Fe if (error) then @@ -3088,7 +3060,7 @@ logical error ! --- calculate disorientation between me and my neighbor --- - do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors + do n = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors neighboring_e = mesh_ipNeighborhood(1,n,i,e) neighboring_i = mesh_ipNeighborhood(2,n,i,e) if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists @@ -3181,7 +3153,7 @@ function crystallite_postResults(& crystallite_postResults(c+1) = real(crystallite_sizePostResults(crystID),pReal) ! size of results from cryst c = c + 1_pInt - do o = 1,crystallite_Noutput(crystID) + do o = 1_pInt,crystallite_Noutput(crystID) mySize = 0_pInt select case(crystallite_output(o,crystID)) case ('phase') diff --git a/code/homogenization.f90 b/code/homogenization.f90 index d4ec99bad..b2c5bbca1 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -72,6 +72,7 @@ CONTAINS !* Module initialization * !************************************** subroutine homogenization_init(Temperature) +use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pReal,pInt use math, only: math_I3 use debug, only: debug_verbosity @@ -288,10 +289,10 @@ subroutine materialpoint_stressAndItsTangent(& crystallite_converged, & crystallite_stressAndItsTangent, & crystallite_orientations - use debug, only: debug_verbosity, & - debug_selectiveDebugger, & +use debug, only: debug_verbosity, & debug_e, & debug_i, & + debug_selectiveDebugger, & debug_MaterialpointLoopDistribution, & debug_MaterialpointStateLoopDistribution use math, only: math_pDecomposition @@ -358,7 +359,7 @@ subroutine materialpoint_stressAndItsTangent(& if ( materialpoint_converged(i,e) ) then #ifndef _OPENMP if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then - write(6,'(a,1x,f10.8,1x,a,1x,f10.8,1x,a,/)') '<< HOMOG >> winding forward from', & + write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,/)') '<< HOMOG >> winding forward from', & materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', & materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent' endif @@ -409,7 +410,7 @@ subroutine materialpoint_stressAndItsTangent(& #ifndef _OPENMP if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then - write(6,'(a,1x,f10.8,/)') & + write(6,'(a,1x,f12.8,/)') & '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',& materialpoint_subStep(i,e) endif @@ -593,7 +594,7 @@ subroutine homogenization_partitionDeformation(& el & ! element ) - use prec, only: pReal,pInt + use prec, only: pInt use mesh, only: mesh_element use material, only: homogenization_type, homogenization_maxNgrains use crystallite, only: crystallite_partionedF0,crystallite_partionedF @@ -634,7 +635,7 @@ function homogenization_updateState(& ip, & ! integration point el & ! element ) - use prec, only: pReal,pInt + use prec, only: pInt use mesh, only: mesh_element use material, only: homogenization_type, homogenization_maxNgrains use crystallite, only: crystallite_P,crystallite_dPdF,crystallite_partionedF,crystallite_partionedF0 ! modified <<>> @@ -682,7 +683,7 @@ subroutine homogenization_averageStressAndItsTangent(& ip, & ! integration point el & ! element ) - use prec, only: pReal,pInt + use prec, only: pInt use mesh, only: mesh_element use material, only: homogenization_type, homogenization_maxNgrains use crystallite, only: crystallite_P,crystallite_dPdF @@ -724,7 +725,7 @@ subroutine homogenization_averageTemperature(& ip, & ! integration point el & ! element ) - use prec, only: pReal,pInt + use prec, only: pInt use mesh, only: mesh_element use material, only: homogenization_type, homogenization_maxNgrains use crystallite, only: crystallite_Temperature diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 087e4b347..2fd104318 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -63,18 +63,19 @@ CONTAINS !* Module initialization * !************************************** subroutine homogenization_RGC_init(& - file & ! file pointer to material configuration + myFile & ! file pointer to material configuration ) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt, pReal use debug, only: debug_verbosity use math, only: math_Mandel3333to66, math_Voigt66to3333,math_I3,math_sampleRandomOri,math_EulerToR,inRad use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use IO use material - integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 4 - integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt), intent(in) :: myFile + integer(pInt), parameter :: maxNchunks = 4_pInt + integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions integer(pInt) section, maxNinstance, i,j,e, output, mySize, myInstance character(len=64) tag character(len=1024) line @@ -86,7 +87,7 @@ subroutine homogenization_RGC_init(& #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - maxNinstance = count(homogenization_type == homogenization_RGC_label) + maxNinstance = int(count(homogenization_type == homogenization_RGC_label),pInt) if (maxNinstance == 0) return allocate(homogenization_RGC_sizeState(maxNinstance)); homogenization_RGC_sizeState = 0_pInt @@ -100,20 +101,20 @@ subroutine homogenization_RGC_init(& allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),maxNinstance)) homogenization_RGC_sizePostResult = 0_pInt allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems)) - forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems) + forall (i = 1_pInt:mesh_maxNips,e = 1_pInt:mesh_NcpElems) homogenization_RGC_orientation(:,:,i,e) = math_I3 end forall - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do ! read thru sections of phase part - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section @@ -212,9 +213,9 @@ subroutine homogenization_RGC_init(& homogenization_RGC_sizeState(i) & - = 3*(homogenization_RGC_Ngrains(1,i)-1_pInt)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & - + 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1_pInt)*homogenization_RGC_Ngrains(3,i) & - + 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1_pInt) & + = 3_pInt*(homogenization_RGC_Ngrains(1,i)-1_pInt)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) & + + 3_pInt*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1_pInt)*homogenization_RGC_Ngrains(3,i) & + + 3_pInt*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1_pInt) & + 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy, ! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component enddo @@ -254,9 +255,9 @@ subroutine homogenization_RGC_partitionDeformation(& ) use prec, only: pReal,pInt,p_vec use debug, only: debug_verbosity - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance - use FEsolving, only: theInc,cycleCounter,theTime + use FEsolving, only: theInc,cycleCounter implicit none @@ -291,7 +292,7 @@ subroutine homogenization_RGC_partitionDeformation(& !* Compute the deformation gradient of individual grains due to relaxations homID = homogenization_typeInstance(mesh_element(3,el)) F = 0.0_pReal - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) do iFace = 1_pInt,nFace intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain @@ -339,12 +340,11 @@ function homogenization_RGC_updateState(& use prec, only: pReal,pInt,p_vec use debug, only: debug_verbosity, debug_e, debug_i use math, only: math_invert - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains,homogenization_typeInstance, & homogenization_Ngrains use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC, & maxdRelax_RGC,viscPower_RGC,viscModus_RGC,refRelaxRate_RGC - use FEsolving, only: theInc,cycleCounter,theTime implicit none @@ -379,8 +379,8 @@ function homogenization_RGC_updateState(& homID = homogenization_typeInstance(mesh_element(3,el)) nGDim = homogenization_RGC_Ngrains(:,homID) nGrain = homogenization_Ngrains(mesh_element(3,el)) - nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1)*nGDim(3) & - + nGDim(1)*nGDim(2)*(nGDim(3)-1) + nIntFaceTot = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) & + + nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) !* Allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster allocate(resid(3_pInt*nIntFaceTot)); resid = 0.0_pReal @@ -390,10 +390,10 @@ function homogenization_RGC_updateState(& drelax = state%p(1:3_pInt*nIntFaceTot) - state0%p(1:3_pInt*nIntFaceTot) !* Debugging the obtained state - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Obtained state: ' - do i = 1,3*nIntFaceTot + do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,2(e15.8,1x))')state%p(i) enddo write(6,*)' ' @@ -407,16 +407,16 @@ function homogenization_RGC_updateState(& call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID) !* Debugging the mismatch, stress and penalties of grains - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) - do iGrain = 1,nGrain + do iGrain = 1_pInt,nGrain write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain) write(6,*)' ' write(6,'(1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain - do i = 1,3 - write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1,3), & - (R(i,j,iGrain), j = 1,3), & - (D(i,j,iGrain), j = 1,3) + do i = 1_pInt,3_pInt + write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1_pInt,3_pInt), & + (R(i,j,iGrain), j = 1_pInt,3_pInt), & + (D(i,j,iGrain), j = 1_pInt,3_pInt) enddo write(6,*)' ' enddo @@ -426,7 +426,7 @@ function homogenization_RGC_updateState(& !* ------------------------------------------------------------------------------------------------------------- !*** Computing the residual stress from the balance of traction at all (interior) interfaces - do iNum = 1,nIntFaceTot + do iNum = 1_pInt,nIntFaceTot faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index) !* Identify the left/bottom/back grain (-|N) @@ -443,23 +443,23 @@ function homogenization_RGC_updateState(& normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal !* Compute the residual of traction at the interface (in local system, 4-dimensional index) - do i = 1,3 + do i = 1_pInt,3_pInt tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1_pInt)))/(refRelaxRate_RGC*dt))**viscPower_RGC, & - drelax(i+3*(iNum-1))) ! contribution from the relaxation viscosity - do j = 1,3 + drelax(i+3*(iNum-1_pInt))) ! contribution from the relaxation viscosity + do j = 1_pInt,3_pInt tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & + (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j) ! contribution from material stress P, mismatch penalty R, and volume penalty D ! projected into the interface - resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array + resid(i+3_pInt*(iNum-1_pInt)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array enddo enddo !* Debugging the residual stress - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum - write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3) + write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt) write(6,*)' ' !$OMP END CRITICAL (write2out) endif @@ -468,13 +468,13 @@ function homogenization_RGC_updateState(& !* ------------------------------------------------------------------------------------------------------------- !*** Convergence check for stress residual - stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress - stresLoc = maxloc(abs(P)) ! get the location of the maximum stress - residMax = maxval(abs(tract)) ! get the maximum of the residual - residLoc = maxloc(abs(tract)) ! get the position of the maximum residual + stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress + stresLoc = int(maxloc(abs(P)),pInt) ! get the location of the maximum stress + residMax = maxval(abs(tract)) ! get the maximum of the residual + residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual !* Debugging the convergent criteria - if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then + if (debug_verbosity == 4_pInt .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) write(6,'(1x,a)')' ' write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el @@ -504,34 +504,34 @@ function homogenization_RGC_updateState(& !* ... all energy densities computed by time-integration constitutiveWork = state%p(3*nIntFaceTot+1) penaltyEnergy = state%p(3*nIntFaceTot+5) - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy - do i = 1,3 - do j = 1,3 - constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/dble(nGrain) - penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/dble(nGrain) + do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy + do i = 1_pInt,3_pInt + do j = 1_pInt,3_pInt + constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) + penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal) enddo enddo enddo state%p(3*nIntFaceTot+1) = constitutiveWork ! the bulk mechanical/constitutive work - state%p(3*nIntFaceTot+2) = sum(NN(1,:))/dble(nGrain) ! the overall mismatch of all interface normal to e1-direction - state%p(3*nIntFaceTot+3) = sum(NN(2,:))/dble(nGrain) ! the overall mismatch of all interface normal to e2-direction - state%p(3*nIntFaceTot+4) = sum(NN(3,:))/dble(nGrain) ! the overall mismatch of all interface normal to e3-direction + state%p(3*nIntFaceTot+2) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction + state%p(3*nIntFaceTot+3) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction + state%p(3*nIntFaceTot+4) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction state%p(3*nIntFaceTot+5) = penaltyEnergy ! the overall penalty energy state%p(3*nIntFaceTot+6) = volDiscrep ! the overall volume discrepancy - state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/dble(3*nIntFaceTot) ! the average rate of relaxation vectors + state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors - if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then + if (debug_verbosity == 4_pInt .and. debug_e == el .and. debug_i == ip) then !$OMP CRITICAL (write2out) write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork - write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), & - sum(NN(2,:))/dble(nGrain), & - sum(NN(3,:))/dble(nGrain) + write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), & + sum(NN(2,:))/real(nGrain,pReal), & + sum(NN(3,:))/real(nGrain,pReal) write(6,'(1x,a30,1x,e15.8)')'Penalty energy: ',penaltyEnergy write(6,'(1x,a30,1x,e15.8)')'Volume discrepancy: ',volDiscrep write(6,*)'' write(6,'(1x,a30,1x,e15.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt - write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot) + write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) write(6,*)'' call flush(6) !$OMP END CRITICAL (write2out) @@ -575,20 +575,20 @@ function homogenization_RGC_updateState(& !* ... of the constitutive stress tangent, !* assembled from dPdF or material constitutive model "smatrix" allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot)); smatrix = 0.0_pReal - do iNum = 1,nIntFaceTot + do iNum = 1_pInt,nIntFaceTot faceID = homogenization_RGC_interface1to4(iNum,homID) ! assembling of local dPdF into global Jacobian matrix !* Identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate into global grain ID - intFaceN = homogenization_RGC_getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system + intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the interface normal - do iFace = 1,nFace + do iFace = 1_pInt,nFace intFaceN = homogenization_RGC_getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface mornN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get normal of the interfaces iMun = homogenization_RGC_interface4to1(intFaceN,homID) ! translate the interfaces ID into local 4-dimensional index if (iMun .gt. 0) then ! get the corresponding tangent - do i=1,3; do j=1,3; do k=1,3; do l=1,3 + do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l) enddo;enddo;enddo;enddo ! projecting the material tangent dPdF into the interface @@ -598,16 +598,16 @@ function homogenization_RGC_updateState(& !* Identify the right/up/front grain (+|P) iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem + iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate sytem iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate into global grain ID - intFaceP = homogenization_RGC_getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system + intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the connecting interface in local coordinate system normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal - do iFace = 1,nFace + do iFace = 1_pInt,nFace intFaceP = homogenization_RGC_getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface mornP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get normal of the interfaces iMun = homogenization_RGC_interface4to1(intFaceP,homID) ! translate the interfaces ID into local 4-dimensional index if (iMun .gt. 0) then ! get the corresponding tangent - do i=1,3; do j=1,3; do k=1,3; do l=1,3 + do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l) enddo;enddo;enddo;enddo endif @@ -615,11 +615,11 @@ function homogenization_RGC_updateState(& enddo !* Debugging the global Jacobian matrix of stress tangent - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of stress' - do i = 1,3*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot) + do i = 1_pInt,3_pInt*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -631,7 +631,7 @@ function homogenization_RGC_updateState(& allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal - do ipert = 1,3*nIntFaceTot + do ipert = 1_pInt,3_pInt*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector state%p(1:3*nIntFaceTot) = p_relax @@ -641,25 +641,25 @@ function homogenization_RGC_updateState(& !* Computing the global stress residual array from the perturbed state p_resid = 0.0_pReal - do iNum = 1,nIntFaceTot + do iNum = 1_pInt,nIntFaceTot faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index) !* Identify the left/bottom/back grain (-|N) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceN = homogenization_RGC_getInterface(2*faceID(1),iGr3N) ! identifying the interface ID of the grain + intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the interface ID of the grain normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the corresponding interface normal !* Identify the right/up/front grain (+|P) iGr3P = iGr3N - iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) + iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index) iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate the local grain ID into global coordinate system (1-dimensional index) - intFaceP = homogenization_RGC_getInterface(2*faceID(1)-1,iGr3P) ! identifying the interface ID of the grain + intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the interface ID of the grain normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the corresponding normal !* Compute the residual stress (contribution of mismatch and volume penalties) from perturbed state at all interfaces - do i = 1,3 - do j = 1,3 + do i = 1_pInt,3_pInt + do j = 1_pInt,3_pInt p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) & + (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) & + (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) & @@ -674,8 +674,8 @@ function homogenization_RGC_updateState(& if (debug_verbosity == 4) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of penalty' - do i = 1,3*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot) + do i = 1_pInt,3_pInt*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -684,18 +684,18 @@ function homogenization_RGC_updateState(& !* ... of the numerical viscosity traction "rmatrix" allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal - forall (i=1:3*nIntFaceTot) & + forall (i=1_pInt:3_pInt*nIntFaceTot) & rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & (abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! tangent due to numerical viscosity traction appears ! only in the main diagonal term !* Debugging the global Jacobian matrix of numerical viscosity tangent - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix of penalty' - do i = 1,3*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot) + do i = 1_pInt,3_pInt*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -708,8 +708,8 @@ function homogenization_RGC_updateState(& if (debug_verbosity == 4) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian matrix (total)' - do i = 1,3*nIntFaceTot - write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot) + do i = 1_pInt,3_pInt*nIntFaceTot + write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' call flush(6) @@ -727,7 +727,7 @@ function homogenization_RGC_updateState(& if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Jacobian inverse' - do i = 1,3*nIntFaceTot + do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot) enddo write(6,*)' ' @@ -737,8 +737,8 @@ function homogenization_RGC_updateState(& !* Calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration drelax = 0.0_pReal - do i = 1,3*nIntFaceTot - do j = 1,3*nIntFaceTot + do i = 1_pInt,3_pInt*nIntFaceTot + do j = 1_pInt,3_pInt*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable enddo enddo @@ -754,10 +754,10 @@ function homogenization_RGC_updateState(& endif !* Debugging the return state - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) write(6,'(1x,a30)')'Returned state: ' - do i = 1,3*nIntFaceTot + do i = 1_pInt,3_pInt*nIntFaceTot write(6,'(1x,2(e15.8,1x))')state%p(i) enddo write(6,*)' ' @@ -785,7 +785,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(& use prec, only: pReal,pInt,p_vec use debug, only: debug_verbosity - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance use math, only: math_Plain3333to99 implicit none @@ -804,13 +804,13 @@ subroutine homogenization_RGC_averageStressAndItsTangent(& Ngrains = homogenization_Ngrains(mesh_element(3,el)) !* Debugging the grain tangent - if (debug_verbosity == 4) then + if (debug_verbosity == 4_pInt) then !$OMP CRITICAL (write2out) - do iGrain = 1,Ngrains + do iGrain = 1_pInt,Ngrains dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain)) write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain - do i = 1,9 - write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1,9) + do i = 1_pInt,9_pInt + write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1_pInt,9_pInt) enddo write(6,*)' ' enddo @@ -819,8 +819,8 @@ subroutine homogenization_RGC_averageStressAndItsTangent(& endif !* Computing the average first Piola-Kirchhoff stress P and the average tangent dPdF - avgP = sum(P,3)/dble(Ngrains) - dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains) + avgP = sum(P,3)/real(Ngrains,pReal) + dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) endsubroutine @@ -834,7 +834,7 @@ function homogenization_RGC_averageTemperature(& ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains, homogenization_Ngrains implicit none @@ -842,11 +842,11 @@ function homogenization_RGC_averageTemperature(& real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature integer(pInt), intent(in) :: ip,el real(pReal) homogenization_RGC_averageTemperature - integer(pInt) homID, Ngrains + integer(pInt) :: Ngrains !* Computing the average temperature Ngrains = homogenization_Ngrains(mesh_element(3,el)) - homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains) + homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal) endfunction @@ -861,7 +861,7 @@ pure function homogenization_RGC_postResults(& use prec, only: pReal,pInt,p_vec use mesh, only: mesh_element - use material, only: homogenization_typeInstance,homogenization_Noutput,homogenization_Ngrains + use material, only: homogenization_typeInstance,homogenization_Noutput implicit none !* Definition of variables @@ -873,34 +873,34 @@ pure function homogenization_RGC_postResults(& homogenization_RGC_postResults homID = homogenization_typeInstance(mesh_element(3,el)) - nIntFaceTot = (homogenization_RGC_Ngrains(1,homID)-1)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID) + & - homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1)*homogenization_RGC_Ngrains(3,homID) + & - homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1) + nIntFaceTot=(homogenization_RGC_Ngrains(1,homID)-1_pInt)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID)& + + homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1_pInt)*homogenization_RGC_Ngrains(3,homID)& + + homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1_pInt) c = 0_pInt homogenization_RGC_postResults = 0.0_pReal - do o = 1,homogenization_Noutput(mesh_element(3,el)) + do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) select case(homogenization_RGC_output(o,homID)) case('constitutivework') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1) - c = c + 1 + c = c + 1_pInt case('magnitudemismatch') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2) homogenization_RGC_postResults(c+2) = state%p(3*nIntFaceTot+3) homogenization_RGC_postResults(c+3) = state%p(3*nIntFaceTot+4) - c = c + 3 + c = c + 3_pInt case('penaltyenergy') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5) - c = c + 1 + c = c + 1_pInt case('volumediscrepancy') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6) - c = c + 1 + c = c + 1_pInt case('averagerelaxrate') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7) - c = c + 1 + c = c + 1_pInt case('maximumrelaxrate') homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8) - c = c + 1 + c = c + 1_pInt end select enddo @@ -960,24 +960,24 @@ subroutine homogenization_RGC_stressPenalty(& !* ------------------------------------------------------------------------------------------------------------- !*** Computing the mismatch and penalty stress tensor of all grains - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) Gmoduli = homogenization_RGC_equivalentModuli(iGrain,ip,el) muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) ! get the grain ID in local 3-dimensional index (x,y,z)-position !* Looping over all six interfaces of each grain - do iFace = 1,nFace + do iFace = 1_pInt,nFace intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the interface normal iGNghb3 = iGrain3 ! identify the neighboring grain across the interface - iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1)))) + iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt) if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) ! with periodicity along e1 direction - if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1 + if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1_pInt if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) ! with periodicity along e2 direction - if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1 + if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1_pInt if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) ! with periodicity along e3 direction - if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1 + if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1_pInt iGNghb = homogenization_RGC_grain3to1(iGNghb3,homID) ! get the ID of the neighboring grain Gmoduli = homogenization_RGC_equivalentModuli(iGNghb,ip,el) ! collecting the shear modulus and Burgers vector of the neighbor muGNghb = Gmoduli(1) @@ -987,10 +987,10 @@ subroutine homogenization_RGC_stressPenalty(& !* Compute the mismatch tensor of all interfaces nDefNorm = 0.0_pReal nDef = 0.0_pReal - do i = 1,3 - do j = 1,3 - do k = 1,3 - do l = 1,3 + do i = 1_pInt,3_pInt + do j = 1_pInt,3_pInt + do k = 1_pInt,3_pInt + do l = 1_pInt,3_pInt nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l)! compute the interface mismatch tensor from the jump of deformation gradient enddo enddo @@ -1011,10 +1011,10 @@ subroutine homogenization_RGC_stressPenalty(& ! endif !* Compute the stress penalty of all interfaces - do i = 1,3 - do j = 1,3 - do k = 1,3 - do l = 1,3 + do i = 1_pInt,3_pInt + do j = 1_pInt,3_pInt + do k = 1_pInt,3_pInt + do l = 1_pInt,3_pInt rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*homogenization_RGC_xiAlpha(homID) & *surfCorr(abs(intFace(1)))/homogenization_RGC_dAlpha(abs(intFace(1)),homID) & *cosh(homogenization_RGC_ciAlpha(homID)*nDefNorm) & @@ -1073,16 +1073,16 @@ subroutine homogenization_RGC_volumePenalty(& !* Compute the volumes of grains and of cluster vDiscrep = math_det33(fAvg) ! compute the volume of the cluster - do iGrain = 1,nGrain + do iGrain = 1_pInt,nGrain gVol(iGrain) = math_det33(fDef(:,:,iGrain)) ! compute the volume of individual grains - vDiscrep = vDiscrep - gVol(iGrain)/dble(nGrain) ! calculate the difference/dicrepancy between + vDiscrep = vDiscrep - gVol(iGrain)/real(nGrain,pReal) ! calculate the difference/dicrepancy between ! the volume of the cluster and the the total volume of grains enddo !* Calculate the stress and penalty due to volume discrepancy vPen = 0.0_pReal - do iGrain = 1,nGrain - vPen(:,:,iGrain) = -1.0_pReal/dble(nGrain)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & + do iGrain = 1_pInt,nGrain + vPen(:,:,iGrain) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* & sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* & gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain))) @@ -1128,11 +1128,11 @@ function homogenization_RGC_surfaceCorrection(& invC = 0.0_pReal call math_invert33(avgC,invC,detF,error) homogenization_RGC_surfaceCorrection = 0.0_pReal - do iBase = 1,3 + do iBase = 1_pInt,3_pInt intFace = (/iBase,1_pInt,1_pInt,1_pInt/) nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of the interface - do i = 1,3 - do j = 1,3 + do i = 1_pInt,3_pInt + do j = 1_pInt,3_pInt homogenization_RGC_surfaceCorrection(iBase) = & ! compute the component of (the inverse of) the stretch in the direction of the normal homogenization_RGC_surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) enddo @@ -1154,8 +1154,6 @@ function homogenization_RGC_equivalentModuli(& use prec, only: pReal,pInt,p_vec use constitutive, only: constitutive_homogenizedC,constitutive_averageBurgers - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_typeInstance implicit none !* Definition of variables @@ -1188,8 +1186,6 @@ function homogenization_RGC_relaxationVector(& ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_typeInstance implicit none !* Definition of variables @@ -1219,7 +1215,6 @@ function homogenization_RGC_interfaceNormal(& use prec, only: pReal,pInt,p_vec use math, only: math_mul33x3 - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips implicit none !* Definition of variables @@ -1263,7 +1258,7 @@ function homogenization_RGC_getInterface(& integer(pInt) iDir !* Direction of interface normal - iDir = (int(dble(iFace-1)/2.0_pReal)+1)*(-1_pInt)**iFace + iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**iFace homogenization_RGC_getInterface(1) = iDir !* Identify the interface position by the direction of its normal @@ -1281,8 +1276,7 @@ function homogenization_RGC_grain1to3(& homID & ! homogenization ID ) - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use prec, only: pInt,p_vec implicit none !* Definition of variables @@ -1306,8 +1300,7 @@ function homogenization_RGC_grain3to1(& homID & ! homogenization ID ) - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use prec, only: pInt,p_vec implicit none !* Definition of variables @@ -1318,7 +1311,7 @@ function homogenization_RGC_grain3to1(& !* Get the grain ID nGDim = homogenization_RGC_Ngrains(:,homID) - homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1) + nGDim(1)*nGDim(2)*(grain3(3)-1) + homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1_pInt) + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt) endfunction @@ -1330,8 +1323,7 @@ function homogenization_RGC_interface4to1(& homID & ! homogenization ID ) - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use prec, only: pInt,p_vec implicit none !* Definition of variables @@ -1342,22 +1334,22 @@ function homogenization_RGC_interface4to1(& nGDim = homogenization_RGC_Ngrains(:,homID) !* Compute the total number of interfaces, which ... - nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1 - nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2 - nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3 + nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1 + nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2 + nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3 !* Get the corresponding interface ID in 1D global array if (abs(iFace4D(1)) == 1_pInt) then ! ... interface with normal //e1 - homogenization_RGC_interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) & - + nGDim(2)*nGDim(3)*(iFace4D(2)-1) + homogenization_RGC_interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) & + + nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt) if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) homogenization_RGC_interface4to1 = 0_pInt elseif (abs(iFace4D(1)) == 2_pInt) then ! ... interface with normal //e2 - homogenization_RGC_interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) & - + nGDim(3)*nGDim(1)*(iFace4D(3)-1) + nIntFace(1) + homogenization_RGC_interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1_pInt) & + + nGDim(3)*nGDim(1)*(iFace4D(3)-1_pInt) + nIntFace(1) if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) homogenization_RGC_interface4to1 = 0_pInt elseif (abs(iFace4D(1)) == 3_pInt) then ! ... interface with normal //e3 - homogenization_RGC_interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) & - + nGDim(1)*nGDim(2)*(iFace4D(4)-1) + nIntFace(1) + nIntFace(2) + homogenization_RGC_interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1_pInt) & + + nGDim(1)*nGDim(2)*(iFace4D(4)-1_pInt) + nIntFace(1) + nIntFace(2) if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) homogenization_RGC_interface4to1 = 0_pInt endif @@ -1372,7 +1364,6 @@ function homogenization_RGC_interface1to4(& ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips implicit none !* Definition of variables @@ -1383,30 +1374,53 @@ function homogenization_RGC_interface1to4(& nGDim = homogenization_RGC_Ngrains(:,homID) !* Compute the total number of interfaces, which ... - nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1 - nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2 - nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3 + nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1 + nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2 + nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3 !* Get the corresponding interface ID in 4D (normal and local position) - if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! ... interface with normal //e1 - homogenization_RGC_interface1to4(1) = 1 + if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! ... interface with normal //e1 + homogenization_RGC_interface1to4(1) = 1_pInt homogenization_RGC_interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt - homogenization_RGC_interface1to4(4) = mod(int(real(iFace1D-1_pInt,pReal)/real(nGDim(2),pReal),pInt),nGDim(3))+1 - homogenization_RGC_interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal),pInt)+1 - elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then - ! ... interface with normal //e2 - homogenization_RGC_interface1to4(1) = 2 - homogenization_RGC_interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1 - homogenization_RGC_interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1 - homogenization_RGC_interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1 - elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then - ! ... interface with normal //e3 - homogenization_RGC_interface1to4(1) = 3 - homogenization_RGC_interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1 - homogenization_RGC_interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/& - real(nGDim(1),pReal),pInt),nGDim(2))+1 - homogenization_RGC_interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/& - real(nGDim(1),pReal)/real(nGDim(2),pReal),pInt)+1 + homogenization_RGC_interface1to4(4) = mod(& + int(& + real(iFace1D-1_pInt,pReal)/& + real(nGDim(2),pReal)& + ,pInt)& + ,nGDim(3))+1_pInt + homogenization_RGC_interface1to4(2) = int(& + real(iFace1D-1_pInt,pReal)/& + real(nGDim(2),pReal)/& + real(nGDim(3),pReal)& + ,pInt)+1_pInt + elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! ... interface with normal //e2 + homogenization_RGC_interface1to4(1) = 2_pInt + homogenization_RGC_interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt + homogenization_RGC_interface1to4(2) = mod(& + int(& + real(iFace1D-nIntFace(1)-1_pInt,pReal)/& + real(nGDim(3),pReal)& + ,pInt)& + ,nGDim(1))+1_pInt + homogenization_RGC_interface1to4(3) = int(& + real(iFace1D-nIntFace(1)-1_pInt,pReal)/& + real(nGDim(3),pReal)/& + real(nGDim(1),pReal)& + ,pInt)+1_pInt + elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! ... interface with normal //e3 + homogenization_RGC_interface1to4(1) = 3_pInt + homogenization_RGC_interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt + homogenization_RGC_interface1to4(3) = mod(& + int(& + real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/& + real(nGDim(1),pReal)& + ,pInt)& + ,nGDim(2))+1_pInt + homogenization_RGC_interface1to4(4) = int(& + real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/& + real(nGDim(1),pReal)/& + real(nGDim(2),pReal)& + ,pInt)+1_pInt endif endfunction @@ -1442,18 +1456,18 @@ subroutine homogenization_RGC_grainDeformation(& integer(pInt), dimension (3) :: iGrain3 integer(pInt) homID, iGrain,iFace,i,j ! - integer(pInt), parameter :: nFace = 6 + integer(pInt), parameter :: nFace = 6_pInt !* Compute the deformation gradient of individual grains due to relaxations homID = homogenization_typeInstance(mesh_element(3,el)) F = 0.0_pReal - do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) + do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) - do iFace = 1,nFace + do iFace = 1_pInt,nFace intFace = homogenization_RGC_getInterface(iFace,iGrain3) aVect = homogenization_RGC_relaxationVector(intFace,state,homID) nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) - forall (i=1:3,j=1:3) & + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient diff --git a/code/homogenization_isostrain.f90 b/code/homogenization_isostrain.f90 index e3c82e389..b5df2e9ce 100644 --- a/code/homogenization_isostrain.f90 +++ b/code/homogenization_isostrain.f90 @@ -59,17 +59,18 @@ CONTAINS !* Module initialization * !************************************** subroutine homogenization_isostrain_init(& - file & ! file pointer to material configuration + myFile & ! file pointer to material configuration ) - use, intrinsic :: iso_fortran_env - use prec, only: pInt, pReal + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: pInt use math, only: math_Mandel3333to66, math_Voigt66to3333 use IO use material - integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 2 - integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) section, maxNinstance, i,j, output, mySize + integer(pInt), intent(in) :: myFile + integer(pInt), parameter :: maxNchunks = 2_pInt + integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions + integer(pInt) section, i, j, output, mySize + integer :: maxNinstance, k !no pInt (stores a system dependen value from 'count' character(len=64) tag character(len=1024) line @@ -81,7 +82,7 @@ subroutine homogenization_isostrain_init(& !$OMP END CRITICAL (write2out) maxNinstance = count(homogenization_type == homogenization_isostrain_label) - if (maxNinstance == 0_pInt) return + if (maxNinstance == 0) return allocate(homogenization_isostrain_sizeState(maxNinstance)) ; homogenization_isostrain_sizeState = 0_pInt allocate(homogenization_isostrain_sizePostResults(maxNinstance)); homogenization_isostrain_sizePostResults = 0_pInt @@ -91,16 +92,16 @@ subroutine homogenization_isostrain_init(& allocate(homogenization_isostrain_output(maxval(homogenization_Noutput), & maxNinstance)) ; homogenization_isostrain_output = '' - rewind(file) + rewind(myFile) line = '' section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do ! read thru sections of phase part - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section @@ -121,18 +122,18 @@ subroutine homogenization_isostrain_init(& endif enddo -100 do i = 1_pInt,maxNinstance ! sanity checks +100 do k = 1,maxNinstance ! sanity checks enddo - do i = 1,maxNinstance + do k = 1,maxNinstance homogenization_isostrain_sizeState(i) = 0_pInt - do j = 1,maxval(homogenization_Noutput) + do j = 1_pInt,maxval(homogenization_Noutput) select case(homogenization_isostrain_output(j,i)) case('ngrains') - mySize = 1 + mySize = 1_pInt case default - mySize = 0 + mySize = 0_pInt end select if (mySize > 0_pInt) then ! any meaningful output found @@ -180,7 +181,7 @@ subroutine homogenization_isostrain_partitionDeformation(& el & ! my element ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains,homogenization_Ngrains implicit none @@ -193,7 +194,7 @@ subroutine homogenization_isostrain_partitionDeformation(& integer(pInt) i ! homID = homogenization_typeInstance(mesh_element(3,el)) - forall (i = 1:homogenization_Ngrains(mesh_element(3,el))) & + forall (i = 1_pInt:homogenization_Ngrains(mesh_element(3,el))) & F(1:3,1:3,i) = avgF return @@ -215,7 +216,6 @@ function homogenization_isostrain_updateState(& ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains implicit none @@ -249,7 +249,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(& ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains, homogenization_Ngrains implicit none @@ -263,8 +263,8 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(& ! homID = homogenization_typeInstance(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el)) - avgP = sum(P,3)/dble(Ngrains) - dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains) + avgP = sum(P,3)/real(Ngrains,pReal) + dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal) return @@ -281,7 +281,7 @@ function homogenization_isostrain_averageTemperature(& ) use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips + use mesh, only: mesh_element use material, only: homogenization_maxNgrains, homogenization_Ngrains implicit none @@ -293,7 +293,7 @@ function homogenization_isostrain_averageTemperature(& ! homID = homogenization_typeInstance(mesh_element(3,el)) Ngrains = homogenization_Ngrains(mesh_element(3,el)) - homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains) + homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal) return @@ -325,11 +325,11 @@ pure function homogenization_isostrain_postResults(& c = 0_pInt homogenization_isostrain_postResults = 0.0_pReal - do o = 1,homogenization_Noutput(mesh_element(3,el)) + do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) select case(homogenization_isostrain_output(o,homID)) case ('ngrains') - homogenization_isostrain_postResults(c+1) = real(homogenization_isostrain_Ngrains(homID),pReal) - c = c + 1 + homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal) + c = c + 1_pInt end select enddo