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

started suggestions for F77 --> F90 style 

xkdyad useful..?

looking for error calculation (Ricardo to rescue here)
This commit is contained in:
Philip Eisenlohr 2010-07-02 17:15:53 +00:00
parent d6ba9d54b6
commit 4f76eada31
1 changed files with 55 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
code/mpie_spectral.f90
View File

@ -282,7 +282,7 @@ program mpie_spectral
real(pReal), dimension(3,3) :: disgradmacro,disgradmacroactual real(pReal), dimension(3,3) :: disgradmacro,disgradmacroactual
real(pReal), dimension(3,3) :: ddisgradmacro,ddisgradmacroacum,ddisgrad,ddisgradim real(pReal), dimension(3,3) :: ddisgradmacro,ddisgradmacroacum,ddisgrad,ddisgradim
real(pReal), dimension(3,3) :: defgrad0,defgrad real(pReal), dimension(3,3) :: defgrad0,defgrad
real(pReal), dimension(3,3) :: udot,scauchy,scauav,aux33 real(pReal), dimension(3,3) :: udot,scauchy,scauav,aux33,xkdyad
integer(pInt), dimension(3) :: nn integer(pInt), dimension(3) :: nn
integer(pInt), dimension(2) :: nn2 integer(pInt), dimension(2) :: nn2
@ -495,6 +495,10 @@ program mpie_spectral
do ielem = 1, prodnn !call each element with identity (math_i3) to initialize with high stress do ielem = 1, prodnn !call each element with identity (math_i3) to initialize with high stress
call CPFEM_general(3,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) call CPFEM_general(3,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
enddo
do ielem = 1, prodnn !call each element with identity (math_i3) to initialize with high stress
call CPFEM_general(2,math_i3,math_i3,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
c066 = c066+dsde c066 = c066+dsde
c0 = c0+math_Mandel66to3333(dsde) c0 = c0+math_Mandel66to3333(dsde)
enddo enddo
@ -513,15 +517,8 @@ program mpie_spectral
! INITIALIZATION BEFORE STARTING WITH LOADINGS ! INITIALIZATION BEFORE STARTING WITH LOADINGS
disgradmacro = 0. disgrad = 0.0_pReal
disgradmacro = 0.0_pReal
do k = 1, c
do j = 1, b
do i = 1, a
defgradold(:,:,i,j,k) = math_i3(:,:)
enddo
enddo
enddo
do jload = 1,N !Loop over loadcases defined in the loadcase file do jload = 1,N !Loop over loadcases defined in the loadcase file
udot(:,:) = bc_velocityGrad(:,:,jload) udot(:,:) = bc_velocityGrad(:,:,jload)
@ -551,41 +548,39 @@ program mpie_spectral
do j = 1, b do j = 1, b
do i = 1, a do i = 1, a
ielem = ielem+1 ielem = ielem+1
disgrad(:,:,i,j,k)=disgradmacro(:,:) defgradold(:,:,i,j,k) = math_I3(:,:) + disgrad(:,:,i,j,k) ! wind forward
disgrad(:,:,i,j,k) = disgradmacro(:,:) ! no fluctuations as guess
defgrad0(:,:)=defgradold(:,:,i,j,k) call CPFEM_general(3,defgradold(:,:,i,j,k),math_i3(:,:)+disgrad(:,:,i,j,k),&
defgrad(:,:)=math_i3(:,:)+disgrad(:,:,i,j,k) ! calculate new defgrad temperature,0.0_pReal,ielem,1_pInt,&
call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde) stress,dsde)
! sg(:,:,i,j,k)=math_Mandel6to33(stress) ? ! sg(:,:,i,j,k)=math_Mandel6to33(stress) ?
enddo enddo
enddo enddo
enddo enddo
ielem = 0 ielem = 0
do k = 1, c !loop over FPs do k = 1, c !loop over FPs
do j = 1, b do j = 1, b
do i = 1, a do i = 1, a
ielem = ielem+1 ielem = ielem+1
disgrad(:,:,i,j,k) = disgradmacro(:,:) call CPFEM_general(1,defgradold(:,:,i,j,k),math_i3(:,:)+disgrad(:,:,i,j,k),&
temperature,0.0_pReal,ielem,1_pInt,&
defgrad0(:,:) = defgradold(:,:,i,j,k) stress,dsde)
defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k)
call CPFEM_general(1,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
sg(:,:,i,j,k) = math_Mandel6to33(stress) sg(:,:,i,j,k) = math_Mandel6to33(stress)
enddo enddo
enddo enddo
enddo enddo
ddisgradmacroacum = 0. ddisgradmacroacum = 0.0_pReal
iter = 0 iter = 0_pInt
erre = 2.*error erre = 2.*error
errs = 2.*error errs = 2.*error
do while(iter <= itmax.and.(errs.gt.error.or.erre.gt.error)) do while(iter <= itmax.and.(errs > error .or. erre > error))
iter = iter+1 iter = iter+1
write(*,*)'ITER = ',iter write(*,*) 'ITER = ',iter
write(*,*) 'DIRECT FFT OF STRESS FIELD' write(*,*) 'DIRECT FFT OF STRESS FIELD'
do ii = 1,3 do ii = 1,3
do jj = 1,3 do jj = 1,3
@ -602,7 +597,7 @@ program mpie_spectral
enddo enddo
enddo enddo
enddo enddo
if(c.gt.1) then if(c > 1) then
call fourn(datafft,nn,3,1) call fourn(datafft,nn,3,1)
else else
call fourn(datafft,nn2,2,1) call fourn(datafft,nn2,2,1)
@ -642,18 +637,25 @@ program mpie_spectral
xk(3) = 0. xk(3) = 0.
endif endif
! if (any(xk /= 0.0_pReal) then
! xknorm = 1.0_pReal/(xk(1)**2+xk(2)**2+xk(3)**2)
! forall (i=1:3,j=1:3) xkdyad(i,j) = xknorm * xk(i)*xk(j) ! the dyad is always used anfd could speed up things by using element-wise multiplication plus summation of array
! endif
xknorm = sqrt(xk(1)**2+xk(2)**2+xk(3)**2) xknorm = sqrt(xk(1)**2+xk(2)**2+xk(3)**2)
if (xknorm.ne.0.) then if (xknorm /= 0.0_pReal) then
do i = 1,3 do i = 1,3
xk2(i) = xk(i)/(xknorm*xknorm*2.*pi) xk2(i) = xk(i)/(xknorm*xknorm*2.*pi) ! seems to be not needed...?
xk(i) = xk(i)/xknorm xk(i) = xk(i)/xknorm
enddo enddo
endif endif
! forall(i=1:3,k=1:3) aux33(i,k) = sum(c0(i,:,k,:)*xkdyad(:,:)
! this is probably equiv with below quad looping
aux33 = 0.0_pReal
do i = 1,3 do i = 1,3
do k = 1,3 do k = 1,3
aux33(i,k) = 0.
do j = 1,3 do j = 1,3
do l = 1,3 do l = 1,3
aux33(i,k) = aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l) aux33(i,k) = aux33(i,k)+c0(i,j,k,l)*xk(j)*xk(l)
@ -670,6 +672,8 @@ program mpie_spectral
! pause ! pause
! endif ! endif
! forall (p=1:3,q=1:3,i=1:3,j=1:3) g1(p,q,i,j) = -aux33(p,i)*xk(q)*xk(j)
! could substitute below quad loop
do p = 1,3 do p = 1,3
do q = 1,3 do q = 1,3
do i = 1,3 do i = 1,3
@ -685,7 +689,7 @@ program mpie_spectral
ddisgrad(i,j) = 0. ddisgrad(i,j) = 0.
ddisgradim(i,j) = 0. ddisgradim(i,j) = 0.
! if(kx.eq.0.and.ky.eq.0.and.kz.eq.0.) goto 4 ! if(kx.eq.0.and.ky.eq.0.and.kz.eq.0.) goto 4
if(.not.(kx.eq.0.and.ky.eq.0.and.kz.eq.0.)) then if(.not.(kx == 0. .and. ky == 0. .and. kz == 0.)) then
do k = 1,3 do k = 1,3
do l = 1,3 do l = 1,3
ddisgrad(i,j) = ddisgrad(i,j)+g1(i,j,k,l)*workfft(k,l,kxx,kyy,kzz) ddisgrad(i,j) = ddisgrad(i,j)+g1(i,j,k,l)*workfft(k,l,kxx,kyy,kzz)
@ -723,7 +727,7 @@ program mpie_spectral
enddo enddo
enddo enddo
if(c.gt.1) then if(c > 1) then
call fourn(datafft,nn,3,-1) call fourn(datafft,nn,3,-1)
else else
call fourn(datafft,nn2,2,-1) call fourn(datafft,nn2,2,-1)
@ -753,53 +757,41 @@ program mpie_spectral
do i = 1,a do i = 1,a
ielem = ielem+1 ielem = ielem+1
call CPFEM_general(3,defgradold(:,:,i,j,k),math_i3(:,:)+disgrad(:,:,i,j,k),&
temperature,0.0_pReal,ielem,1_pInt,&
stress,dsde)
defgrad0(:,:) = defgradold(:,:,i,j,k)
defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k)
call CPFEM_general(3,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
sg(:,:,i,j,k) = math_Mandel6to33(stress)
enddo enddo
enddo enddo
enddo enddo
ielem = 0 ielem = 0
scauav = 0.
do k = 1,c do k = 1,c
do j = 1,b do j = 1,b
do i = 1,a do i = 1,a
ielem = ielem+1 ielem = ielem+1
call CPFEM_general(2,defgradold(:,:,i,j,k),math_i3(:,:)+disgrad(:,:,i,j,k),&
defgrad0(:,:) = defgradold(:,:,i,j,k) temperature,0.0_pReal,ielem,1_pInt,&
defgrad(:,:) = math_i3(:,:)+disgrad(:,:,i,j,k) stress,dsde)
call CPFEM_general(2,defgrad0,defgrad,temperature,0.0_pReal,ielem,1_pInt,stress,dsde)
sg(:,:,i,j,k) = math_Mandel6to33(stress) sg(:,:,i,j,k) = math_Mandel6to33(stress)
scauav(:,:) = scauav(:,:)+sg(:,:,i,j,k) ! average stress
enddo enddo
enddo enddo
enddo enddo
!!! call get_smacro scauav = scauav*wgt ! final weighting
! AVERAGE STRESS
scauav = 0.
do k = 1,c
do j = 1,b
do i = 1,a
scauav(:,:) = scauav(:,:)+sg(:,:,i,j,k)*wgt
enddo
enddo
enddo
! MIXED BC ! MIXED BC
do i = 1,3 do i = 1,3
do j = 1,3 do j = 1,3
ddisgradmacro(i,j) = 0. ddisgradmacro(i,j) = 0.
if(iudot(i,j).eq.0) then if(iudot(i,j) == 0) then
! ddisgradmacro(i,j) = ddisgradmacro(i,j)+sum(s0(i,j,:,:)*iscau(:,:)*(scauchy(:,:)-scauav(:,:)))=
! could replace the k,l loop
do k = 1,3 do k = 1,3
do l = 1,3 do l = 1,3
ddisgradmacro(i,j) = ddisgradmacro(i,j)+s0(i,j,k,l)*iscau(k,l)*(scauchy(k,l)-scauav(k,l)) ddisgradmacro(i,j) = ddisgradmacro(i,j)+s0(i,j,k,l)*iscau(k,l)*(scauchy(k,l)-scauav(k,l))
@ -829,13 +821,6 @@ program mpie_spectral
!! write(*,*) 'scauav(1,1),scauav(2,2),scauav(3,3)' !! write(*,*) 'scauav(1,1),scauav(2,2),scauav(3,3)'
!! write(*,*) scauav(1,1),scauav(2,2),scauav(3,3) !! write(*,*) scauav(1,1),scauav(2,2),scauav(3,3)
do k = 1,c
do j = 1,b
do i = 1,a
defgradold(:,:,i,j,k) = math_i3(:,:)+disgrad(:,:,i,j,k)
enddo
enddo
enddo
enddo ! IMICRO ENDDO enddo ! IMICRO ENDDO
enddo ! JLOAD ENDDO Ende looping over loadcases enddo ! JLOAD ENDDO Ende looping over loadcases