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

corrected syntax errors (long lines, line continuation by \) and logical mistake in mpie_cpfem_marc / abq_std which went unnoticed in ifort. Using SunStudio f90 surfaced those...

This commit is contained in:
Philip Eisenlohr 2010-08-03 23:47:00 +00:00
parent 2c4f5c37dc
commit 298cecbfec
11 changed files with 98 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
code/FEsolving.f90
View File

@ -43,10 +43,10 @@
rewind(fileunit) rewind(fileunit)
do do
read (fileunit,'(a1024)',END=100) line read (fileunit,'(a1024)',END=100) line
positions = IO_stringPos(line,1) positions(1:1+2*1) = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,positions,1)) == 'solver' ) then if( IO_lc(IO_stringValue(line,positions,1)) == 'solver' ) then
read (fileunit,'(a1024)',END=100) line ! Garbage line read (fileunit,'(a1024)',END=100) line ! Garbage line
positions = IO_stringPos(line,2) positions(1:1+2*2) = IO_stringPos(line,2)
symmetricSolver = (IO_intValue(line,positions,2) /= 1_pInt) symmetricSolver = (IO_intValue(line,positions,2) /= 1_pInt)
exit exit
endif endif

2
code/IO.f90
View File

@ -135,7 +135,7 @@ end function
use mpie_interface use mpie_interface
implicit none implicit none
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
character(len=*) relPath character(len=*) relPath
integer(pInt) unit integer(pInt) unit

66
code/constitutive_dislotwin.f90
View File

@ -453,7 +453,7 @@ do i = 1,maxNinstance
'resolved_stress_slip', & 'resolved_stress_slip', &
'threshold_stress_slip', & 'threshold_stress_slip', &
'edge_dipole_distance', & 'edge_dipole_distance', &
'stress_exponent' & 'stress_exponent' &
) )
mySize = constitutive_dislotwin_totalNslip(i) mySize = constitutive_dislotwin_totalNslip(i)
case('twin_fraction', & case('twin_fraction', &
@ -817,7 +817,7 @@ forall (s = 1:ns) &
state(g,ip,el)%p(5*ns+3*nt+s) = & state(g,ip,el)%p(5*ns+3*nt+s) = &
constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*& constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*&
sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),& sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),&
constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance)))
!* threshold stress for growing twin !* threshold stress for growing twin
forall (t = 1:nt) & forall (t = 1:nt) &
@ -920,7 +920,8 @@ do f = 1,lattice_maxNslipFamily ! loop over all
constitutive_dislotwin_v0PerSlipSystem(f,myInstance) constitutive_dislotwin_v0PerSlipSystem(f,myInstance)
!* Shear rates due to slip !* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*&
sign(1.0_pReal,tau_slip(j))
!* Derivatives of shear rates !* Derivatives of shear rates
dgdot_dtauslip(j) = & dgdot_dtauslip(j) = &
@ -953,10 +954,10 @@ do f = 1,lattice_maxNtwinFamily ! loop over all
!* Resolved shear stress on twin system !* Resolved shear stress on twin system
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure))
!* Stress ratios !* Stress ratios
StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance) StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance)
!* Shear rates and their derivatives due to twin !* Shear rates and their derivatives due to twin
if ( tau_twin(j) > 0.0_pReal ) then if ( tau_twin(j) > 0.0_pReal ) then
gdot_twin(j) = & gdot_twin(j) = &
(constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*lattice_shearTwin(index_myFamily+i,myStructure)*& (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*lattice_shearTwin(index_myFamily+i,myStructure)*&
@ -1008,10 +1009,13 @@ function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,g,ip,el)
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use math, only: pi use math, only: pi
use mesh, only: mesh_NcpElems,mesh_maxNips use mesh, only: mesh_NcpElems, mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance use material, only: homogenization_maxNgrains, material_phase, phase_constitutionInstance
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & use lattice, only: lattice_Sslip,lattice_Sslip_v, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin lattice_Stwin,lattice_Stwin_v, &
lattice_maxNslipFamily,lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem, &
lattice_shearTwin
implicit none implicit none
!* Input-Output variables !* Input-Output variables
@ -1029,7 +1033,8 @@ real(pReal), dimension(constitutive_dislotwin_totalNslip(phase_constitutionInsta
gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,& gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,&
ClimbVelocity,DotRhoEdgeDipClimb,DotRhoDipFormation ClimbVelocity,DotRhoEdgeDipClimb,DotRhoDipFormation
real(pReal), dimension(constitutive_dislotwin_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: gdot_twin,tau_twin real(pReal), dimension(constitutive_dislotwin_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
gdot_twin,tau_twin
!* Shortened notation !* Shortened notation
myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myInstance = phase_constitutionInstance(material_phase(g,ip,el))
@ -1064,7 +1069,8 @@ do f = 1,lattice_maxNslipFamily ! loop over all
constitutive_dislotwin_v0PerSlipSystem(f,myInstance) constitutive_dislotwin_v0PerSlipSystem(f,myInstance)
!* Shear rates due to slip !* Shear rates due to slip
gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*&
sign(1.0_pReal,tau_slip(j))
!* Multiplication !* Multiplication
DotRhoMultiplication(j) = abs(gdot_slip(j))/& DotRhoMultiplication(j) = abs(gdot_slip(j))/&
@ -1079,8 +1085,8 @@ do f = 1,lattice_maxNslipFamily ! loop over all
EdgeDipDistance(j) = & EdgeDipDistance(j) = &
(3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/& (3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/&
(16.0_pReal*pi*abs(tau_slip(j))) (16.0_pReal*pi*abs(tau_slip(j)))
if (EdgeDipDistance(j)>state(g,ip,el)%p(4*ns+2*nt+j)) EdgeDipDistance(j)=state(g,ip,el)%p(4*ns+2*nt+j) if (EdgeDipDistance(j)>state(g,ip,el)%p(4*ns+2*nt+j)) EdgeDipDistance(j)=state(g,ip,el)%p(4*ns+2*nt+j)
if (EdgeDipDistance(j)<EdgeDipMinDistance) EdgeDipDistance(j)=EdgeDipMinDistance if (EdgeDipDistance(j)<EdgeDipMinDistance) EdgeDipDistance(j)=EdgeDipMinDistance
DotRhoDipFormation(j) = & DotRhoDipFormation(j) = &
((2.0_pReal*EdgeDipDistance(j))/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& ((2.0_pReal*EdgeDipDistance(j))/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*&
state(g,ip,el)%p(j)*abs(gdot_slip(j)) state(g,ip,el)%p(j)*abs(gdot_slip(j))
@ -1102,14 +1108,14 @@ do f = 1,lattice_maxNslipFamily ! loop over all
VacancyDiffusion = & VacancyDiffusion = &
constitutive_dislotwin_D0(myInstance)*exp(-constitutive_dislotwin_Qsd(myInstance)/(kB*Temperature)) constitutive_dislotwin_D0(myInstance)*exp(-constitutive_dislotwin_Qsd(myInstance)/(kB*Temperature))
if (tau_slip(j) == 0.0_pReal) then if (tau_slip(j) == 0.0_pReal) then
DotRhoEdgeDipClimb(j) = 0.0_pReal DotRhoEdgeDipClimb(j) = 0.0_pReal
else else
ClimbVelocity(j) = & ClimbVelocity(j) = &
((3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*VacancyDiffusion*AtomicVolume)/(2.0_pReal*pi*kB*Temperature))*& ((3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*VacancyDiffusion*AtomicVolume)/(2.0_pReal*pi*kB*Temperature))*&
(1/(EdgeDipDistance(j)+EdgeDipMinDistance)) (1/(EdgeDipDistance(j)+EdgeDipMinDistance))
DotRhoEdgeDipClimb(j) = & DotRhoEdgeDipClimb(j) = &
(4.0_pReal*ClimbVelocity(j)*state(g,ip,el)%p(ns+j))/(EdgeDipDistance(j)-EdgeDipMinDistance) (4.0_pReal*ClimbVelocity(j)*state(g,ip,el)%p(ns+j))/(EdgeDipDistance(j)-EdgeDipMinDistance)
endif endif
!* Edge dislocation density rate of change !* Edge dislocation density rate of change
constitutive_dislotwin_dotState(j) = & constitutive_dislotwin_dotState(j) = &
@ -1280,17 +1286,17 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
case ('threshold_stress_slip') case ('threshold_stress_slip')
constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+3*nt+1):(6*ns+3*nt)) constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+3*nt+1):(6*ns+3*nt))
c = c + ns c = c + ns
case ('edge_dipole_distance') case ('edge_dipole_distance')
j = 0_pInt j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) = & constitutive_dislotwin_postResults(c+j) = &
(3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/& (3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/&
(16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)))) (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))))
constitutive_dislotwin_postResults(c+j) = min(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j)) constitutive_dislotwin_postResults(c+j) = min(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j))
! constitutive_dislotwin_postResults(c+j) = max(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j)) ! constitutive_dislotwin_postResults(c+j) = max(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j))
enddo; enddo enddo; enddo
c = c + ns c = c + ns
case ('twin_fraction') case ('twin_fraction')
@ -1326,8 +1332,8 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
if (nt > 0_pInt) then if (nt > 0_pInt) then
j = 0_pInt j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all slip families do f = 1,lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family
do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure))
enddo; enddo enddo; enddo
@ -1366,11 +1372,11 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal)
!* Stress exponent !* Stress exponent
if (gdot_slip==0.0_pReal) then if (gdot_slip==0.0_pReal) then
constitutive_dislotwin_postResults(c+j) = 0.0_pReal constitutive_dislotwin_postResults(c+j) = 0.0_pReal
else else
constitutive_dislotwin_postResults(c+j) = (tau/gdot_slip)*dgdot_dtauslip constitutive_dislotwin_postResults(c+j) = (tau/gdot_slip)*dgdot_dtauslip
endif endif
enddo ; enddo enddo ; enddo
c = c + ns c = c + ns

8
code/constitutive_nonlocal.f90
View File

@ -583,11 +583,11 @@ do i = 1,maxNinstance
constitutive_nonlocal_forestProjectionEdge(s1, s2, i) & constitutive_nonlocal_forestProjectionEdge(s1, s2, i) &
= abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_st(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane lattice_st(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
constitutive_nonlocal_forestProjectionScrew(s1, s2, i) & constitutive_nonlocal_forestProjectionScrew(s1, s2, i) &
= abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
!*** calculation of interaction matrices !*** calculation of interaction matrices
@ -1128,7 +1128,7 @@ ns = constitutive_nonlocal_totalNslip(myInstance)
forall (t = 1:8) & forall (t = 1:8) &
rhoSgl(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns) rhoSgl(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
@ -1473,7 +1473,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_fluxdensity = fluxdensity ! ... then copy flux density to neighbor to ensure zero gradient in fluxdensity neighboring_fluxdensity = fluxdensity ! ... then copy flux density to neighbor to ensure zero gradient in fluxdensity
endif endif
else ! if no neighbor existent... else ! if no neighbor existent...
if ( all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,n,ip,el))))) ) then ! ... and we want periodic fluxes at surface... if ( all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(:,n,ip,el))))) ) then ! ... and we want periodic fluxes at surface...
forall (t = 1:4) & ! ... then mirror fluxes forall (t = 1:4) & ! ... then mirror fluxes
neighboring_fluxdensity(:,t) = fluxdensity(:,t-1_pInt+2_pInt*mod(t,2_pInt)) neighboring_fluxdensity(:,t) = fluxdensity(:,t-1_pInt+2_pInt*mod(t,2_pInt))
else ! ... and we have a free surface... else ! ... and we have a free surface...

9
code/constitutive_titanmod.f90
View File

@ -15,8 +15,8 @@ character(len=18), dimension(2), parameter:: constitutive_titanmod_listBasicSlip
'rho_screw'/) 'rho_screw'/)
character(len=18), dimension(1), parameter:: constitutive_titanmod_listBasicTwinStates = (/'twinFraction'/) character(len=18), dimension(1), parameter:: constitutive_titanmod_listBasicTwinStates = (/'twinFraction'/)
character(len=18), dimension(5), parameter:: constitutive_titanmod_listDependentSlipStates =(/'invLambdaSlipe', & character(len=18), dimension(5), parameter:: constitutive_titanmod_listDependentSlipStates =(/'invLambdaSlipe', &
'invLambdaSlips', & 'invLambdaSlips', &
'etauSlipThreshold', & 'etauSlipThreshold', &
'stauSlipThreshold', & 'stauSlipThreshold', &
'invLambdaSlipTwin'/) 'invLambdaSlipTwin'/)
@ -1500,7 +1500,8 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
!* Stress ratios !* Stress ratios
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_titanmod_pe_PerSlipSystem(j,myInstance) StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_titanmod_pe_PerSlipSystem(j,myInstance)
StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_titanmod_pe_PerSlipSystem(j,myInstance)-1.0_pReal) StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**&
(constitutive_titanmod_pe_PerSlipSystem(j,myInstance)-1.0_pReal)
!* Boltzmann ratio !* Boltzmann ratio
BoltzmannRatio = constitutive_titanmod_f0_PerSlipSystem(j,myInstance)/(kB*Temperature) BoltzmannRatio = constitutive_titanmod_f0_PerSlipSystem(j,myInstance)/(kB*Temperature)
!* Initial shear rates !* Initial shear rates
@ -1511,7 +1512,7 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
!* Shear rates due to slip !* Shear rates due to slip
constitutive_titanmod_postResults(c+j) = & constitutive_titanmod_postResults(c+j) = &
DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_titanmod_qe_PerSlipSystem(j,myInstance))* & DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_titanmod_qe_PerSlipSystem(j,myInstance))* &
sign(1.0_pReal,tau) sign(1.0_pReal,tau)
enddo ; enddo enddo ; enddo
! invLambdaSlipe', & ! invLambdaSlipe', &

36
code/crystallite.f90
View File

@ -941,10 +941,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do g = 1,myNgrains do g = 1,myNgrains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) ! stress integration crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) ! stress integration
if ( .not. crystallite_localConstitution(g,i,e) & if ( .not. crystallite_localConstitution(g,i,e) &
.and. .not. crystallite_todo(g,i,e)) & ! if broken non-local... .and. .not. crystallite_todo(g,i,e)) & ! if broken non-local...
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -954,7 +954,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_todo(g,i,e)) & if (crystallite_todo(g,i,e)) &
constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState constitutive_dotState(g,i,e)%p = 0.0_pReal ! zero out dotState
enddo; enddo; enddo enddo; enddo; enddo
crystallite_statedamper = 1.0_pReal crystallite_statedamper = 1.0_pReal
@ -991,21 +991,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
do g = 1,myNgrains do g = 1,myNgrains
selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g) selectiveDebugger = (e == debug_e .and. i == debug_i .and. g == debug_g)
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature
crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) & crystallite_converged(g,i,e) = crystallite_stateConverged(g,i,e) &
.and. crystallite_temperatureConverged(g,i,e) .and. crystallite_temperatureConverged(g,i,e)
if ( .not. crystallite_localConstitution(g,i,e) & if ( .not. crystallite_localConstitution(g,i,e) &
.and. .not. crystallite_todo(g,i,e)) & ! if updateState signals broken non-local... .and. .not. crystallite_todo(g,i,e)) & ! if updateState signals broken non-local...
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! ...all non-locals skipped
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged? if (any(.not. crystallite_converged .and. .not. crystallite_localConstitution)) & ! any non-local not yet converged?
crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged crystallite_converged = crystallite_converged .and. crystallite_localConstitution ! all non-local not converged
crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged crystallite_todo = crystallite_todo .and. .not. crystallite_converged ! skip all converged
enddo enddo
@ -1013,10 +1013,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,myNgrains do g = 1,myNgrains
if (crystallite_converged(g,i,e)) then ! if stiffness calculation converged... if (crystallite_converged(g,i,e)) then ! if stiffness calculation converged...
crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! ... use tangent dP_ij/dFg_kl crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! ... use tangent dP_ij/dFg_kl
elseif (.not. storedConvergenceFlag(g,i,e)) then ! if crystallite didnÕt converge before... elseif (.not. storedConvergenceFlag(g,i,e)) then ! if crystallite didnÕt converge before...
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! ... use (elastic) fallback crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! ... use (elastic) fallback
endif endif
enddo; enddo; enddo enddo; enddo; enddo
@ -1436,7 +1436,8 @@ LpLoop: do
if (error) then if (error) then
if (verboseDebugger .and. selectiveDebugger) then if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e,' ; iteration ', NiterationStress write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e, &
' ; iteration ', NiterationStress
write(6,*) write(6,*)
write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',dRdLp write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',dRdLp
write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive
@ -1469,7 +1470,8 @@ LpLoop: do
if (error) then if (error) then
if (verboseDebugger .and. selectiveDebugger) then if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e,' ; iteration ', NiterationStress write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e, &
' ; iteration ', NiterationStress
write(6,*) write(6,*)
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)

4
code/homogenization_RGC.f90
View File

@ -476,7 +476,9 @@ function homogenization_RGC_updateState(&
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt
if (el == 1 .and. ip == 1) then if (el == 1 .and. ip == 1) then
write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork write(6,'(x,a30,x,e14.8)')'Constitutive work: ',constitutiveWork
write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain),sum(NN(2,:))/dble(nGrain),sum(NN(3,:))/dble(nGrain) write(6,'(x,a30,3(x,e14.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), &
sum(NN(2,:))/dble(nGrain), &
sum(NN(3,:))/dble(nGrain)
write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy write(6,'(x,a30,x,e14.8)')'Penalty energy: ',penaltyEnergy
write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep write(6,'(x,a30,x,e14.8)')'Volume discrepancy: ',volDiscrep
write(6,*)'' write(6,*)''

20
code/mesh.f90
View File

@ -2006,7 +2006,7 @@ subroutine mesh_marc_count_cpSizes (unit)
rewind(unit) rewind(unit)
do do
read (unit,610,END=630) line read (unit,610,END=630) line
pos = IO_stringPos(line,1) pos = IO_stringPos(line,maxNchunks)
if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then
read (unit,610,END=630) line ! Garbage line read (unit,610,END=630) line ! Garbage line
do i=1,mesh_Nelems ! read all elements do i=1,mesh_Nelems ! read all elements
@ -2060,7 +2060,7 @@ subroutine mesh_marc_count_cpSizes (unit)
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
pos = IO_stringPos(line,2) pos = IO_stringPos(line,maxNchunks)
if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
@ -2348,7 +2348,7 @@ subroutine mesh_marc_count_cpSizes (unit)
do while (e < mesh_NcpElems) do while (e < mesh_NcpElems)
read(unit,'(a1024)',END=110) line read(unit,'(a1024)',END=110) line
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
pos = IO_stringPos(line,1) pos(1:1+2*1) = IO_stringPos(line,1)
e = e+1 ! valid element entry e = e+1 ! valid element entry
mesh_element ( 1,e) = e ! FE id mesh_element ( 1,e) = e ! FE id
@ -2399,7 +2399,7 @@ subroutine mesh_marc_count_cpSizes (unit)
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
pos = IO_stringPos(line,1) pos(1:1+2*1) = IO_stringPos(line,1)
if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then if( IO_lc(IO_stringValue(line,pos,1)) == 'connectivity' ) then
read (unit,610,END=620) line ! Garbage line read (unit,610,END=620) line ! Garbage line
do i = 1,mesh_Nelems do i = 1,mesh_Nelems
@ -2421,16 +2421,16 @@ subroutine mesh_marc_count_cpSizes (unit)
620 rewind(unit) ! just in case "initial state" apears before "connectivity" 620 rewind(unit) ! just in case "initial state" apears before "connectivity"
read (unit,610,END=620) line read (unit,610,END=620) line
do do
pos = IO_stringPos(line,2) pos(1:1+2*2) = IO_stringPos(line,2)
if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial') .and. & if( (IO_lc(IO_stringValue(line,pos,1)) == 'initial') .and. &
(IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then (IO_lc(IO_stringValue(line,pos,2)) == 'state') ) then
if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style if (initialcondTableStyle == 2) read (unit,610,END=620) line ! read extra line for new style
read (unit,610,END=630) line ! read line with index of state var read (unit,610,END=630) line ! read line with index of state var
pos = IO_stringPos(line,1) pos(1:1+2*1) = IO_stringPos(line,1)
sv = IO_IntValue(line,pos,1) ! figure state variable index sv = IO_IntValue(line,pos,1) ! figure state variable index
if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest if( (sv == 2).or.(sv == 3) ) then ! only state vars 2 and 3 of interest
read (unit,610,END=620) line ! read line with value of state var read (unit,610,END=620) line ! read line with value of state var
pos = IO_stringPos(line,1) pos(1:1+2*1) = IO_stringPos(line,1)
do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value? do while (scan(IO_stringValue(line,pos,1),'+-',back=.true.)>1) ! is noEfloat value?
val = NINT(IO_fixedNoEFloatValue(line,(/0,20/),1)) ! state var's value val = NINT(IO_fixedNoEFloatValue(line,(/0,20/),1)) ! state var's value
mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of homogenization and microstructure index mesh_maxValStateVar(sv-1) = max(val,mesh_maxValStateVar(sv-1)) ! remember max val of homogenization and microstructure index
@ -2445,7 +2445,7 @@ subroutine mesh_marc_count_cpSizes (unit)
enddo enddo
if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style if (initialcondTableStyle == 0) read (unit,610,END=620) line ! ignore IP range for old table style
read (unit,610,END=630) line read (unit,610,END=630) line
pos = IO_stringPos(line,1) pos(1:1+2*1) = IO_stringPos(line,1)
enddo enddo
endif endif
else else
@ -2487,7 +2487,7 @@ subroutine mesh_marc_count_cpSizes (unit)
rewind(unit) rewind(unit)
do do
read (unit,610,END=620) line read (unit,610,END=620) line
pos = IO_stringPos(line,2) pos(1:1+2*2) = IO_stringPos(line,2)
if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true. if ( IO_lc(IO_stringValue(line,pos,1)) == '*part' ) inPart = .true.
if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. & if ( IO_lc(IO_stringValue(line,pos,1)) == '*end' .and. &
IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false. IO_lc(IO_stringValue(line,pos,2)) == 'part' ) inPart = .false.
@ -2539,7 +2539,7 @@ subroutine mesh_marc_count_cpSizes (unit)
if (i <= mesh_Nmaterials) then ! found one? if (i <= mesh_Nmaterials) then ! found one?
elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet elemSetName = mesh_mapMaterial(i) ! take corresponding elemSet
read (unit,610,END=630) line ! read homogenization and microstructure read (unit,610,END=630) line ! read homogenization and microstructure
pos = IO_stringPos(line,2) pos(1:1+2*2) = IO_stringPos(line,2)
homog = NINT(IO_floatValue(line,pos,1)) homog = NINT(IO_floatValue(line,pos,1))
micro = NINT(IO_floatValue(line,pos,2)) micro = NINT(IO_floatValue(line,pos,2))
do i = 1,mesh_NelemSets ! look thru all elemSet definitions do i = 1,mesh_NelemSets ! look thru all elemSet definitions

4
code/mpie_cpfem_abaqus_std.f
View File

@ -174,7 +174,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect) calcMode(npt,cp_en) = .not. calcMode(npt,cp_en) ! ping pong (calc <--> collect)
if ( calcMode(npt,cp_en) ) then ! now calc if ( calcMode(npt,cp_en) ) then ! now calc
if ( lastMode .ne. calcMode(npt,cp_en) ) then ! first after ping pong if ( lastMode .neqv. calcMode(npt,cp_en) ) then ! first after ping pong
call debug_reset() ! resets debugging call debug_reset() ! resets debugging
outdatedFFN1 = .false. outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1 cycleCounter = cycleCounter + 1
@ -186,7 +186,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
computationMode = 2 ! plain calc computationMode = 2 ! plain calc
endif endif
else ! now collect else ! now collect
if ( lastMode .ne. calcMode(npt,cp_en) ) call debug_info() ! first after ping pong reports debugging if ( lastMode .neqv. calcMode(npt,cp_en) ) call debug_info() ! first after ping pong reports debugging
if ( lastIncConverged ) then if ( lastIncConverged ) then
lastIncConverged = .false. lastIncConverged = .false.
computationMode = 4 ! collect and backup Jacobian after convergence computationMode = 4 ! collect and backup Jacobian after convergence

8
code/mpie_cpfem_marc.f90
View File

@ -57,7 +57,7 @@ end subroutine
function getSolverWorkingDirectoryName() function getSolverWorkingDirectoryName()
implicit none implicit none
character(1024) getSolverWorkingDirectoryName, outName character(1024) getSolverWorkingDirectoryName, outName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
getSolverWorkingDirectoryName='' getSolverWorkingDirectoryName=''
outName='' outName=''
@ -71,7 +71,7 @@ function getSolverJobName()
implicit none implicit none
character(1024) getSolverJobName, outName character(1024) getSolverJobName, outName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer(pInt) extPos integer(pInt) extPos
getSolverJobName='' getSolverJobName=''
@ -262,7 +262,7 @@ subroutine hypela2(&
calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect) calcMode(nn,cp_en) = .not. calcMode(nn,cp_en) ! ping pong (calc <--> collect)
if ( calcMode(nn,cp_en) ) then ! now calc if ( calcMode(nn,cp_en) ) then ! now calc
if ( lastMode .ne. calcMode(nn,cp_en) ) then ! first after ping pong if ( lastMode .neqv. calcMode(nn,cp_en) ) then ! first after ping pong
call debug_reset() ! resets debugging call debug_reset() ! resets debugging
outdatedFFN1 = .false. outdatedFFN1 = .false.
cycleCounter = cycleCounter + 1 cycleCounter = cycleCounter + 1
@ -274,7 +274,7 @@ subroutine hypela2(&
computationMode = 2 ! plain calc computationMode = 2 ! plain calc
endif endif
else ! now collect else ! now collect
if ( lastMode /= calcMode(nn,cp_en) ) call debug_info() ! first after ping pong reports debugging if ( lastMode .neqv. calcMode(nn,cp_en) ) call debug_info() ! first after ping pong reports debugging
if ( lastIncConverged ) then if ( lastIncConverged ) then
lastIncConverged = .false. lastIncConverged = .false.
computationMode = 4 ! collect and backup Jacobian after convergence computationMode = 4 ! collect and backup Jacobian after convergence

9
code/mpie_spectral.f90
View File

@ -55,7 +55,7 @@ function getSolverWorkingDirectoryName()
implicit none implicit none
character(len=1024) cwd,outname,getSolverWorkingDirectoryName character(len=1024) cwd,outname,getSolverWorkingDirectoryName
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
call getarg(2,outname) ! path to loadFile call getarg(2,outname) ! path to loadFile
@ -83,7 +83,7 @@ function getSolverJobName()
implicit none implicit none
character(1024) getSolverJobName, outName, cwd character(1024) getSolverJobName, outName, cwd
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer(pInt) posExt,posSep integer(pInt) posExt,posSep
getSolverJobName = '' getSolverJobName = ''
@ -119,9 +119,10 @@ function getLoadcaseName()
implicit none implicit none
character(len=1024) getLoadcaseName, outName, cwd character(len=1024) getLoadcaseName, outName, cwd
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \ character(len=*), parameter :: pathSep = achar(47)//achar(92) ! forward and backward slash
integer(pInt) posExt,posSep integer(pInt) posExt,posSep
posExt = 0 !not sure if its needed
posExt = 0 ! not sure if required
call getarg(2,getLoadcaseName) call getarg(2,getLoadcaseName)
posExt = scan(getLoadcaseName,'.',back=.true.) posExt = scan(getLoadcaseName,'.',back=.true.)