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

small polishing to be compatible with NAG compiler (basis for automatic differentiation)

This commit is contained in:
Philip Eisenlohr 2008-04-29 12:18:17 +00:00
parent 9626c25bfb
commit 51e596c81d
3 changed files with 55 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
trunk/CPFEM_Taylor.f90
View File

@ -145,7 +145,7 @@
!
cp_en = mesh_FEasCP('elem',CPFEM_en)
if (cp_en == 1 .and. CPFEM_in == 1) &
write(6,'(a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') &
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
'mode',CPFEM_mode
!
@ -156,7 +156,7 @@
CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new
endif
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt
!
@ -177,8 +177,10 @@
H_bar = 0.0_pReal
forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
H_bar(i,j,k,l) = H_bar(i,j,k,l) + &
(CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en)*CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en)*CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - &
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + &
CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en) * &
CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en) * &
CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - &
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en) + &
0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
@ -206,6 +208,7 @@
!
END SUBROUTINE
!
!
!**********************************************************
!*** calculate the material point behaviour ***
!**********************************************************
@ -249,13 +252,13 @@
CPFEM_dt,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),&
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en),CPFEM_ffn_bar(:,:,CPFEM_in,cp_en),&
CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en))
if (msg /= 'ok') then ! solution not reached --> exit
write(6,*) 'grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in
call IO_error(600)
return
endif
CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en),constitutive_state_old(:,grain,CPFEM_in,cp_en))
if (msg /= 'ok') then ! solution not reached --> exit
write(6,*) 'grain loop failed to converge @ EL:',cp_en,' IP:',CPFEM_in
call IO_error(600)
return
endif
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en)
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + volfrac*PK1

30
trunk/constitutive_pheno.f90
View File

@ -427,7 +427,7 @@ fileunit=200
!* First reading: number of materials and textures
!-----------------------------
!* determine material_maxN and texture_maxN from last respective parts
if(IO_open_file(fileunit,filename)==.false.) goto 100
if(.not. IO_open_file(fileunit,filename)) call IO_error (200) ! corrupt mattex file
part = '_dummy_'
do while (part/='')
formerPart = part
@ -440,16 +440,16 @@ do while (part/='')
end select
enddo
!* Array allocation
allocate(material_CrystalStructure(material_maxN)) ; material_CrystalStructure=0_pInt
allocate(material_CrystalStructure(material_maxN)) ; material_CrystalStructure=0_pInt
allocate(material_Nslip(material_maxN)) ; material_Nslip=0_pInt
allocate(material_C11(material_maxN)) ; material_C11=0.0_pReal
allocate(material_C12(material_maxN)) ; material_C12=0.0_pReal
allocate(material_C13(material_maxN)) ; material_C13=0.0_pReal
allocate(material_C33(material_maxN)) ; material_C33=0.0_pReal
allocate(material_C44(material_maxN)) ; material_C44=0.0_pReal
allocate(material_Cslip_66(6,6,material_maxN)) ; material_Cslip_66=0.0_pReal
allocate(material_Cslip_66(6,6,material_maxN)) ; material_Cslip_66=0.0_pReal
allocate(material_s0_slip(material_maxN)) ; material_s0_slip=0.0_pReal
allocate(material_gdot0_slip(material_maxN)) ; material_gdot0_slip=0.0_pReal
allocate(material_gdot0_slip(material_maxN)) ; material_gdot0_slip=0.0_pReal
allocate(material_n_slip(material_maxN)) ; material_n_slip=0.0_pReal
allocate(material_h0(material_maxN)) ; material_h0=0.0_pReal
allocate(material_s_sat(material_maxN)) ; material_s_sat=0.0_pReal
@ -515,16 +515,16 @@ do i=1,material_maxN
case(3) ! hcp
material_Cslip_66(1,1,i)=material_C11(i)
material_Cslip_66(2,2,i)=material_C11(i)
material_Cslip_66(3,3,i)=material_C33(i)
material_Cslip_66(1,2,i)=material_C12(i)
material_Cslip_66(2,1,i)=material_C12(i)
material_Cslip_66(1,3,i)=material_C13(i)
material_Cslip_66(3,1,i)=material_C13(i)
material_Cslip_66(2,3,i)=material_C13(i)
material_Cslip_66(3,2,i)=material_C13(i)
material_Cslip_66(4,4,i)=material_C44(i)
material_Cslip_66(5,5,i)=material_C44(i)
material_Cslip_66(6,6,i)=0.5_pReal*(material_C11(i)-material_C12(i))
material_Cslip_66(3,3,i)=material_C33(i)
material_Cslip_66(1,2,i)=material_C12(i)
material_Cslip_66(2,1,i)=material_C12(i)
material_Cslip_66(1,3,i)=material_C13(i)
material_Cslip_66(3,1,i)=material_C13(i)
material_Cslip_66(2,3,i)=material_C13(i)
material_Cslip_66(3,2,i)=material_C13(i)
material_Cslip_66(4,4,i)=material_C44(i)
material_Cslip_66(5,5,i)=material_C44(i)
material_Cslip_66(6,6,i)=0.5_pReal*(material_C11(i)-material_C12(i))
end select
material_Cslip_66(:,:,i) = math_Mandel3333to66(math_Voigt66to3333(material_Cslip_66(:,:,i)))
! Check
@ -534,7 +534,7 @@ enddo
! MISSING some consistency checks may be..?
! if ODFfile present then set NGauss NFiber =0
return
100 call IO_error(200) ! corrupt materials_textures file
end subroutine

49
trunk/crystallite.f90
View File

@ -95,7 +95,8 @@ CONTAINS
cuttedBack = .false. ! no cut back required
guessNew = .false. ! keep the Lp
subFrac = subFrac + subStep
subStep = 1.0_pReal - subFrac ! try one go
subStep = 1.0_pReal - subFrac ! try one go
if (debugger) write (6,*) '--------- one go -----------++##'
else
nCutbacks = nCutbacks + 1 ! record additional cutback
@ -103,7 +104,8 @@ CONTAINS
guessNew = .true. ! redo plastic Lp guess
subStep = subStep / 2.0_pReal ! cut time step in half
state_bestguess = state_current ! current state is then best guess
if (debugger) write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<'
if (debugger) write (6,*) '>>>>>>>>>>>>>>>>>>>> cutback <<<<<<<<<<<<<<<<<<<<<<'
endif
enddo ! potential substepping
!
@ -160,7 +162,8 @@ CONTAINS
use constitutive, only: constitutive_Nstatevars,&
constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,&
constitutive_Nresults,constitutive_Microstructure,constitutive_post_results
use math
use math
use IO
implicit none
!
@ -207,8 +210,8 @@ Outer: do ! outer iteration: State
!
iInner = 0_pInt
leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step
maxleap = 1024.0_pReal ! preassign maximum acceleration level
maxleap = 1024.0_pReal ! preassign maximum acceleration level
!
Lpguess_old = Lpguess ! consider present Lpguess good
!
Inner: do ! inner iteration: Lp
@ -217,8 +220,8 @@ Inner: do ! inner iteration: Lp
msg = 'limit Inner iteration'
debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1
return
endif
endif
!
B = math_i3 - dt*Lpguess
BT = transpose(B)
AB = matmul(A,B)
@ -229,12 +232,10 @@ Inner: do ! inner iteration: Lp
forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure
call constitutive_LpAndItsTangent(Lp,dLp, &
Tstar_v,state,Temperature,grain,ip,cp_en)
Rinner = Lpguess - Lp ! update current residuum
if (not(any(Rinner.NE.Rinner)) .and. & ! exclude any NaN in residuum
!
Rinner = Lpguess - Lp ! update current residuum
!
if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum
( (maxval(abs(Rinner)) < abstol_Inner) .or. & ! below abs tol .or.
( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and.
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner & ! below rel tol
@ -243,20 +244,21 @@ Inner: do ! inner iteration: Lp
) &
exit Inner ! convergence
!
! check for acceleration/deceleration in Newton--Raphson correction
if (any(Rinner.NE.Rinner) .and. & ! NaN occured at regular speed
leapfrog == 1.0) then
Lpguess = Lpguess_old ! restore known good guess
msg = 'NaN present' ! croak for cutback
return
! check for acceleration/deceleration in Newton--Raphson correction
!
if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed
leapfrog == 1.0) then
Lpguess = Lpguess_old ! restore known good guess
msg = 'NaN present' ! croak for cutback
return
elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ?
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
any(Rinner.NE.Rinner) ) then ! NaN
any(Rinner/=Rinner) ) then ! NaN
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
else ! better residuum
dTdLp = 0.0_pReal ! calc dT/dLp
forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
@ -273,7 +275,8 @@ Inner: do ! inner iteration: Lp
!
Rinner_old = Rinner ! remember current residuum
Lpguess_old = Lpguess ! remember current Lp guess
if (iInner > 1 .and. leapfrog < maxleap) &
if (iInner > 1 .and. leapfrog < maxleap) &
leapfrog = 2.0_pReal * leapfrog ! accelerate if ok
endif
!