223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

following philips suggestion for 8c1fe4fff0 

and removing unneeded if statements for worldrank
This commit is contained in:
Martin Diehl 2016-10-25 22:54:32 +02:00
parent e49081e56e
commit 2b50494951
14 changed files with 179 additions and 242 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
code/DAMASK_spectral.f90
View File

@ -192,7 +192,7 @@ program DAMASK_spectral
if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1_pInt) & ! sanity check if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1_pInt) & ! sanity check
call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase
allocate (loadCases(N_n)) ! array of load cases allocate (loadCases(N_n)) ! array of load cases
loadCases%P%myType='p' loadCases%stress%myType='stress'
do i = 1, size(loadCases) do i = 1, size(loadCases)
allocate(loadCases(i)%ID(nActiveFields)) allocate(loadCases(i)%ID(nActiveFields))
@ -244,10 +244,10 @@ program DAMASK_spectral
temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk
if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable
enddo enddo
loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) loadCases(currentLoadCase)%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3]))
loadCases(currentLoadCase)%P%maskFloat = merge(ones,zeros,& loadCases(currentLoadCase)%stress%maskFloat = merge(ones,zeros,&
loadCases(currentLoadCase)%P%maskLogical) loadCases(currentLoadCase)%stress%maskLogical)
loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector) loadCases(currentLoadCase)%stress%values = math_plain9to33(temp_valueVector)
case('t','time','delta') ! increment time case('t','time','delta') ! increment time
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt) loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt)
case('n','incs','increments','steps') ! number of increments case('n','incs','increments','steps') ! number of increments
@ -318,16 +318,16 @@ program DAMASK_spectral
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. & if (any(loadCases(currentLoadCase)%stress%maskLogical .eqv. &
loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only
if (any(loadCases(currentLoadCase)%P%maskLogical .and. & if (any(loadCases(currentLoadCase)%stress%maskLogical .and. &
transpose(loadCases(currentLoadCase)%P%maskLogical) .and. & transpose(loadCases(currentLoadCase)%stress%maskLogical) .and. &
reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
errorID = 838_pInt ! no rotation is allowed by stress BC errorID = 838_pInt ! no rotation is allowed by stress BC
write(6,'(2x,a)') 'stress / GPa:' write(6,'(2x,a)') 'stress / GPa:'
do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt
if(loadCases(currentLoadCase)%P%maskLogical(i,j)) then if(loadCases(currentLoadCase)%stress%maskLogical(i,j)) then
write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%P%values(i,j)*1e-9_pReal write(6,'(2x,f12.7)',advance='no') loadCases(currentLoadCase)%stress%values(i,j)*1e-9_pReal
else else
write(6,'(2x,12a)',advance='no') ' * ' write(6,'(2x,12a)',advance='no') ' * '
endif endif
@ -524,20 +524,20 @@ program DAMASK_spectral
case (DAMASK_spectral_SolverBasicPETSc_label) case (DAMASK_spectral_SolverBasicPETSc_label)
call BasicPETSc_forward (& call BasicPETSc_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime, & guess,timeinc,timeIncOld,remainingLoadCaseTime, &
F_BC = loadCases(currentLoadCase)%deformation, & deformation_BC = loadCases(currentLoadCase)%deformation, &
P_BC = loadCases(currentLoadCase)%P, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverAL_label) case (DAMASK_spectral_SolverAL_label)
call AL_forward (& call AL_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime, & guess,timeinc,timeIncOld,remainingLoadCaseTime, &
F_BC = loadCases(currentLoadCase)%deformation, & deformation_BC = loadCases(currentLoadCase)%deformation, &
P_BC = loadCases(currentLoadCase)%P, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverPolarisation_label) case (DAMASK_spectral_SolverPolarisation_label)
call Polarisation_forward (& call Polarisation_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime, & guess,timeinc,timeIncOld,remainingLoadCaseTime, &
F_BC = loadCases(currentLoadCase)%deformation, & deformation_BC = loadCases(currentLoadCase)%deformation, &
P_BC = loadCases(currentLoadCase)%P, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
end select end select
@ -558,19 +558,19 @@ program DAMASK_spectral
case (DAMASK_spectral_SolverBasicPETSc_label) case (DAMASK_spectral_SolverBasicPETSc_label)
solres(field) = BasicPETSC_solution (& solres(field) = BasicPETSC_solution (&
incInfo,timeinc,timeIncOld, & incInfo,timeinc,timeIncOld, &
P_BC = loadCases(currentLoadCase)%P, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverAL_label) case (DAMASK_spectral_SolverAL_label)
solres(field) = AL_solution (& solres(field) = AL_solution (&
incInfo,timeinc,timeIncOld, & incInfo,timeinc,timeIncOld, &
P_BC = loadCases(currentLoadCase)%P, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
case (DAMASK_spectral_SolverPolarisation_label) case (DAMASK_spectral_SolverPolarisation_label)
solres(field) = Polarisation_solution (& solres(field) = Polarisation_solution (&
incInfo,timeinc,timeIncOld, & incInfo,timeinc,timeIncOld, &
P_BC = loadCases(currentLoadCase)%P, & stress_BC = loadCases(currentLoadCase)%stress, &
rotation_BC = loadCases(currentLoadCase)%rotation) rotation_BC = loadCases(currentLoadCase)%rotation)
end select end select

8
code/kinematics_hydrogen_strain.f90
View File

@ -67,8 +67,6 @@ subroutine kinematics_hydrogen_strain_init(fileUnit)
KINEMATICS_hydrogen_strain_ID, & KINEMATICS_hydrogen_strain_ID, &
material_Nphase, & material_Nphase, &
MATERIAL_partPhase MATERIAL_partPhase
use numerics,only: &
worldrank
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -79,11 +77,9 @@ subroutine kinematics_hydrogen_strain_init(fileUnit)
tag = '', & tag = '', &
line = '' line = ''
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init -+>>>'
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt) maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return

10
code/porosity_none.f90
View File

@ -23,21 +23,17 @@ subroutine porosity_none_init()
use IO, only: & use IO, only: &
IO_timeStamp IO_timeStamp
use material use material
use numerics, only: &
worldrank
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
homog, & homog, &
NofMyHomog NofMyHomog
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_none_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_none_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myhomog: if (porosity_type(homog) == POROSITY_none_ID) then myhomog: if (porosity_type(homog) == POROSITY_none_ID) then
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)

6
code/spectral_damage.f90
View File

@ -98,11 +98,9 @@ subroutine spectral_damage_init()
DMRestoreGlobalVector, & DMRestoreGlobalVector, &
SNESVISetVariableBounds SNESVISetVariableBounds
mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc

104
code/spectral_mech_AL.f90
View File

@ -49,7 +49,7 @@ module spectral_mech_AL
F_av = 0.0_pReal, & !< average incompatible def grad field F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general
character(len=1024), private :: incInfo !< time and increment information character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: & real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
@ -132,11 +132,9 @@ subroutine AL_init
SNESSetConvergenceTest, & SNESSetConvergenceTest, &
SNESSetFromOptions SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
@ -244,7 +242,7 @@ end subroutine AL_init
!> @brief solution for the AL scheme with internal iterations !> @brief solution for the AL scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function &
AL_solution(incInfoIn,timeinc,timeinc_old,P_BC,rotation_BC) AL_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -267,7 +265,7 @@ type(tSolutionState) function &
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment timeinc_old !< increment in time of last increment
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
P_BC stress_BC
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
@ -285,7 +283,7 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) then if (update_gamma) then
call Utilities_updateGamma(C_minMaxAvg,restartWrite) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
@ -294,11 +292,11 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
mask_stress = P_BC%maskFloat mask_stress = stress_BC%maskFloat
params%P_BC = P_BC%values params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
@ -326,8 +324,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
itmax, & itmax, &
itmin, & itmin, &
polarAlpha, & polarAlpha, &
polarBeta, & polarBeta
worldrank
use mesh, only: & use mesh, only: &
grid3, & grid3, &
grid grid
@ -406,16 +403,14 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
totalIter = totalIter + 1_pInt totalIter = totalIter + 1_pInt
if (worldrank == 0_pInt) then write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim)
math_transpose33(F_aim) flush(6)
flush(6)
endif
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -490,8 +485,7 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
err_curl_tolRel, & err_curl_tolRel, &
err_curl_tolAbs, & err_curl_tolAbs, &
err_stress_tolAbs, & err_stress_tolAbs, &
err_stress_tolRel, & err_stress_tolRel
worldrank
use math, only: & use math, only: &
math_mul3333xx33 math_mul3333xx33
use FEsolving, only: & use FEsolving, only: &
@ -514,9 +508,9 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + &
mask_stress *(P_av - params%P_BC))) ! mask = 0.0 for no bc mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! error calculation ! error calculation
@ -538,24 +532,22 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
if (worldrank == 0_pInt) then write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')'
err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')'
err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')'
err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' write(6,'(/,a)') ' ==========================================================================='
write(6,'(/,a)') ' ===========================================================================' flush(6)
flush(6)
endif
end subroutine AL_converged end subroutine AL_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
@ -583,8 +575,8 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
timeinc, & timeinc, &
loadCaseTime !< remaining time of current load case loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
P_BC, & stress_BC, &
F_BC deformation_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
logical, intent(in) :: & logical, intent(in) :: &
guess guess
@ -600,21 +592,19 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
F => xx_psc(0:8,:,:,:) F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:) F_lambda => xx_psc(9:17,:,:,:)
if (restartWrite) then if (restartWrite) then
if (worldrank == 0_pInt) then write(6,'(/,a)') ' writing converged results for restart'
write(6,'(/,a)') ' writing converged results for restart' flush(6)
flush(6)
endif
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F write (777,rec=1) F
close (777) close (777)
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc write (777,rec=1) F_lastInc
close (777) close (777)
call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file
write (777,rec=1) F_lambda write (777,rec=1) F_lambda
close (777) close (777)
call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lambda_lastInc write (777,rec=1) F_lambda_lastInc
close (777) close (777)
if (worldrank == 0_pInt) then if (worldrank == 0_pInt) then
@ -648,14 +638,14 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values f_aimDot = deformation_BC%maskFloat * deformation_BC%values
elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed
f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime
endif endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

95
code/spectral_mech_Basic.f90
View File

@ -122,11 +122,9 @@ subroutine basicPETSc_init
SNESSetConvergenceTest, & SNESSetConvergenceTest, &
SNESSetFromOptions SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
@ -164,7 +162,7 @@ subroutine basicPETSc_init
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
restart: if (restartInc > 1_pInt) then restart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file' 'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6) flush(6)
@ -219,7 +217,7 @@ end subroutine basicPETSc_init
!> @brief solution for the Basic PETSC scheme with internal iterations !> @brief solution for the Basic PETSC scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function &
basicPETSc_solution(incInfoIn,timeinc,timeinc_old,P_BC,rotation_BC) basicPETSc_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -240,7 +238,7 @@ type(tSolutionState) function &
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment timeinc_old !< increment in time of last increment
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
P_BC stress_BC
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
@ -258,17 +256,17 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
mask_stress = P_BC%maskFloat mask_stress = stress_BC%maskFloat
params%P_BC = P_BC%values params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
@ -296,8 +294,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin itmin
use numerics, only: &
worldrank
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
@ -349,16 +345,14 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
totalIter = totalIter + 1_pInt totalIter = totalIter + 1_pInt
if (worldrank == 0_pInt) then write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim)
math_transpose33(F_aim) flush(6)
flush(6)
endif
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -371,8 +365,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim_lastIter = F_aim F_aim_lastIter = F_aim
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc err_stress = maxval(abs(mask_stress * (P_av - params%stress_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme
@ -400,8 +394,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
err_div_tolRel, & err_div_tolRel, &
err_div_tolAbs, & err_div_tolAbs, &
err_stress_tolRel, & err_stress_tolRel, &
err_stress_tolAbs, & err_stress_tolAbs
worldrank
use FEsolving, only: & use FEsolving, only: &
terminallyIll terminallyIll
@ -435,22 +428,20 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
if (worldrank == 0_pInt) then write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & err_div/divTol, ' (',err_div,' / m, tol =',divTol,')'
err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')'
err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' write(6,'(/,a)') ' ==========================================================================='
write(6,'(/,a)') ' ===========================================================================' flush(6)
flush(6)
endif
end subroutine BasicPETSc_converged end subroutine BasicPETSc_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33 ,& math_mul33x33 ,&
math_rotate_backward33 math_rotate_backward33
@ -476,8 +467,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
timeinc, & timeinc, &
loadCaseTime !< remaining time of current load case loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
P_BC, & stress_BC, &
F_BC deformation_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
logical, intent(in) :: & logical, intent(in) :: &
guess guess
@ -490,10 +481,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restart information for spectral solver ! restart information for spectral solver
if (restartWrite) then if (restartWrite) then
if (worldrank == 0_pInt) then write(6,'(/,a)') ' writing converged results for restart'
write(6,'(/,a)') ' writing converged results for restart' flush(6)
flush(6)
endif
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F write (777,rec=1) F
@ -525,14 +514,14 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values f_aimDot = deformation_BC%maskFloat * deformation_BC%values
elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed
f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime
endif endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -547,8 +536,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update local deformation gradient ! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine BasicPETSc_forward end subroutine BasicPETSc_forward

100
code/spectral_mech_Polarisation.f90
View File

@ -132,11 +132,9 @@ subroutine Polarisation_init
SNESSetConvergenceTest, & SNESSetConvergenceTest, &
SNESSetFromOptions SNESSetFromOptions
mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
@ -176,7 +174,7 @@ subroutine Polarisation_init
F => xx_psc(0:8,:,:,:) F => xx_psc(0:8,:,:,:)
F_tau => xx_psc(9:17,:,:,:) F_tau => xx_psc(9:17,:,:,:)
restart: if (restartInc > 1_pInt) then restart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file' 'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6) flush(6)
@ -244,7 +242,7 @@ end subroutine Polarisation_init
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function & type(tSolutionState) function &
Polarisation_solution(incInfoIn,timeinc,timeinc_old,P_BC,rotation_BC) Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC)
use IO, only: & use IO, only: &
IO_error IO_error
use numerics, only: & use numerics, only: &
@ -267,7 +265,7 @@ type(tSolutionState) function &
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment timeinc_old !< increment in time of last increment
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
P_BC stress_BC
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
@ -285,7 +283,7 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) then if (update_gamma) then
call Utilities_updateGamma(C_minMaxAvg,restartWrite) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
@ -294,11 +292,11 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
mask_stress = P_BC%maskFloat mask_stress = stress_BC%maskFloat
params%P_BC = P_BC%values params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
@ -326,8 +324,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
itmax, & itmax, &
itmin, & itmin, &
polarAlpha, & polarAlpha, &
polarBeta, & polarBeta
worldrank
use mesh, only: & use mesh, only: &
grid3, & grid3, &
grid grid
@ -406,16 +403,14 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new iteration ! report begin of new iteration
totalIter = totalIter + 1_pInt totalIter = totalIter + 1_pInt
if (worldrank == 0_pInt) then write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim)
math_transpose33(F_aim) flush(6)
flush(6)
endif
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -489,8 +484,7 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
err_curl_tolRel, & err_curl_tolRel, &
err_curl_tolAbs, & err_curl_tolAbs, &
err_stress_tolAbs, & err_stress_tolAbs, &
err_stress_tolRel, & err_stress_tolRel
worldrank
use math, only: & use math, only: &
math_mul3333xx33 math_mul3333xx33
use FEsolving, only: & use FEsolving, only: &
@ -513,9 +507,9 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + & err_BC = maxval(abs((-mask_stress+1.0_pReal)*math_mul3333xx33(C_scale,F_aim-F_av) + &
mask_stress *(P_av - params%P_BC))) ! mask = 0.0 for no bc mask_stress *(P_av - params%stress_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! error calculation ! error calculation
@ -537,24 +531,22 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
if (worldrank == 0_pInt) then write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')'
err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')' write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')'
err_div/divTol, ' (',err_div, ' / m, tol =',divTol,')' write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
write(6,' (a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')'
err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')' write(6,'(/,a)') ' ==========================================================================='
write(6,'(/,a)') ' ===========================================================================' flush(6)
flush(6)
endif
end subroutine Polarisation_converged end subroutine Polarisation_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
@ -582,8 +574,8 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
timeinc, & timeinc, &
loadCaseTime !< remaining time of current load case loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
P_BC, & stress_BC, &
F_BC deformation_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
logical, intent(in) :: & logical, intent(in) :: &
guess guess
@ -599,19 +591,19 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
F => xx_psc(0:8,:,:,:) F => xx_psc(0:8,:,:,:)
F_tau => xx_psc(9:17,:,:,:) F_tau => xx_psc(9:17,:,:,:)
if (restartWrite) then if (restartWrite) then
if (worldrank == 0_pInt) write(6,'(/,a)') ' writing converged results for restart' write(6,'(/,a)') ' writing converged results for restart'
flush(6) flush(6)
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F write (777,rec=1) F
close (777) close (777)
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc write (777,rec=1) F_lastInc
close (777) close (777)
call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file
write (777,rec=1) F_tau write (777,rec=1) F_tau
close (777) close (777)
call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_tau_lastInc write (777,rec=1) F_tau_lastInc
close (777) close (777)
if (worldrank == 0_pInt) then if (worldrank == 0_pInt) then
@ -645,14 +637,14 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) f_aimDot = deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values f_aimDot = deformation_BC%maskFloat * deformation_BC%values
elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed elseif(deformation_BC%myType=='f') then ! aim at end of load case is prescribed
f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime f_aimDot = deformation_BC%maskFloat * (deformation_BC%values -F_aim)/loadCaseTime
endif endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old if (guess) f_aimDot = f_aimDot + stress_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

4
code/spectral_utilities.f90
View File

@ -85,7 +85,7 @@ module spectral_utilities
type, public :: tLoadCase type, public :: tLoadCase
real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC
type(tBoundaryCondition) :: P, & !< stress BC type(tBoundaryCondition) :: stress, & !< stress BC
deformation !< deformation BC (Fdot or L) deformation !< deformation BC (Fdot or L)
real(pReal) :: time = 0.0_pReal !< length of increment real(pReal) :: time = 0.0_pReal !< length of increment
integer(pInt) :: incs = 0_pInt, & !< number of increments integer(pInt) :: incs = 0_pInt, & !< number of increments
@ -97,7 +97,7 @@ module spectral_utilities
end type tLoadCase end type tLoadCase
type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase including mask type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase including mask
real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal), dimension(3,3) :: stress_BC, rotation_BC
real(pReal) :: timeinc real(pReal) :: timeinc
real(pReal) :: timeincOld real(pReal) :: timeincOld
real(pReal) :: density real(pReal) :: density

8
code/thermal_adiabatic.f90
View File

@ -74,8 +74,6 @@ subroutine thermal_adiabatic_init(fileUnit)
temperature, & temperature, &
temperatureRate, & temperatureRate, &
material_partHomogenization material_partHomogenization
use numerics,only: &
worldrank
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -88,11 +86,9 @@ subroutine thermal_adiabatic_init(fileUnit)
tag = '', & tag = '', &
line = '' line = ''
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_ADIABATIC_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt) maxNinstance = int(count(thermal_type == THERMAL_adiabatic_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return

8
code/thermal_conduction.f90
View File

@ -75,8 +75,6 @@ subroutine thermal_conduction_init(fileUnit)
temperature, & temperature, &
temperatureRate, & temperatureRate, &
material_partHomogenization material_partHomogenization
use numerics,only: &
worldrank
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -89,11 +87,9 @@ subroutine thermal_conduction_init(fileUnit)
tag = '', & tag = '', &
line = '' line = ''
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_CONDUCTION_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(thermal_type == THERMAL_conduction_ID),pInt) maxNinstance = int(count(thermal_type == THERMAL_conduction_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return

10
code/thermal_isothermal.f90
View File

@ -23,8 +23,6 @@ subroutine thermal_isothermal_init()
use IO, only: & use IO, only: &
IO_timeStamp IO_timeStamp
use material use material
use numerics, only: &
worldrank
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
@ -32,13 +30,11 @@ subroutine thermal_isothermal_init()
NofMyHomog, & NofMyHomog, &
sizeState sizeState
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- thermal_'//THERMAL_isothermal_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myhomog: if (thermal_type(homog) == THERMAL_isothermal_ID) then myhomog: if (thermal_type(homog) == THERMAL_isothermal_ID) then
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)

8
code/vacancyflux_cahnhilliard.f90
View File

@ -90,8 +90,6 @@ subroutine vacancyflux_cahnhilliard_init(fileUnit)
vacancyflux_initialCv, & vacancyflux_initialCv, &
material_partHomogenization, & material_partHomogenization, &
material_partPhase material_partPhase
use numerics,only: &
worldrank
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -104,11 +102,9 @@ subroutine vacancyflux_cahnhilliard_init(fileUnit)
tag = '', & tag = '', &
line = '' line = ''
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID),pInt) maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return

8
code/vacancyflux_isochempot.f90
View File

@ -72,8 +72,6 @@ subroutine vacancyflux_isochempot_init(fileUnit)
vacancyConcRate, & vacancyConcRate, &
vacancyflux_initialCv, & vacancyflux_initialCv, &
material_partHomogenization material_partHomogenization
use numerics,only: &
worldrank
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
@ -86,11 +84,9 @@ subroutine vacancyflux_isochempot_init(fileUnit)
tag = '', & tag = '', &
line = '' line = ''
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isochempot_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isochempot_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_isochempot_ID),pInt) maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_isochempot_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return

10
code/vacancyflux_isoconc.f90
View File

@ -23,21 +23,17 @@ subroutine vacancyflux_isoconc_init()
use IO, only: & use IO, only: &
IO_timeStamp IO_timeStamp
use material use material
use numerics, only: &
worldrank
implicit none implicit none
integer(pInt) :: & integer(pInt) :: &
homog, & homog, &
NofMyHomog NofMyHomog
mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isoconc_label//' init -+>>>'
write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isoconc_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
endif mainProcess
initializeInstances: do homog = 1_pInt, material_Nhomogenization initializeInstances: do homog = 1_pInt, material_Nhomogenization
myhomog: if (vacancyflux_type(homog) == VACANCYFLUX_isoconc_ID) then myhomog: if (vacancyflux_type(homog) == VACANCYFLUX_isoconc_ID) then
NofMyHomog = count(material_homog == homog) NofMyHomog = count(material_homog == homog)