From a5c228fd024dd5343dcc88331b3f4e2285c45dd9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Sep 2010 08:05:11 +0000 Subject: [PATCH] changed back transform to complex-to-real, removed redundant variables, reduced size of arrays where possible --- code/makefile | 46 +++++------ code/mpie_spectral.f90 | 179 ++++++++++++++++++++--------------------- 2 files changed, 112 insertions(+), 113 deletions(-) diff --git a/code/makefile b/code/makefile index 4e8b39690..ab0f9697d 100644 --- a/code/makefile +++ b/code/makefile @@ -1,22 +1,22 @@ cpspectral.out: mpie_spectral.o CPFEM.a - ifort -o cpspectral.out mpie_spectral.o CPFEM.a fft.o libfftw3.a constitutive.a advanced.a basics.a + ifort -o cpspectral.out mpie_spectral.o CPFEM.a libfftw3_threads.a libfftw3.a constitutive.a advanced.a basics.a -lpthread mpie_spectral.o: mpie_spectral.f90 CPFEM.o - ifort -c -heap-arrays 500000000 mpie_spectral.f90 + ifort -c -O3 -heap-arrays 500000000 mpie_spectral.f90 CPFEM.a: CPFEM.o ar rc CPFEM.a homogenization.o homogenization_RGC.o homogenization_isostrain.o crystallite.o CPFEM.o constitutive.o CPFEM.o: CPFEM.f90 homogenization.o - ifort -c -heap-arrays 500000000 CPFEM.f90 + ifort -c -O3 -heap-arrays 500000000 CPFEM.f90 homogenization.o: homogenization.f90 homogenization_isostrain.o homogenization_RGC.o crystallite.o - ifort -c -heap-arrays 500000000 homogenization.f90 + ifort -c -O3 -heap-arrays 500000000 homogenization.f90 homogenization_RGC.o: homogenization_RGC.f90 constitutive.a - ifort -c -heap-arrays 500000000 homogenization_RGC.f90 + ifort -c -O3 -heap-arrays 500000000 homogenization_RGC.f90 homogenization_isostrain.o: homogenization_isostrain.f90 basics.a advanced.a - ifort -c -heap-arrays 500000000 homogenization_isostrain.f90 + ifort -c -O3 -heap-arrays 500000000 homogenization_isostrain.f90 crystallite.o: crystallite.f90 constitutive.a - ifort -c -heap-arrays 500000000 crystallite.f90 + ifort -c -O3 -heap-arrays 500000000 crystallite.f90 @@ -24,22 +24,22 @@ constitutive.a: constitutive.o ar rc constitutive.a constitutive.o constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o basics.a advanced.a constitutive.o: constitutive.f90 constitutive_titanmod.o constitutive_nonlocal.o constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o - ifort -c -heap-arrays 500000000 constitutive.f90 + ifort -c -O3 -heap-arrays 500000000 constitutive.f90 constitutive_titanmod.o: constitutive_titanmod.f90 basics.a advanced.a - ifort -c -heap-arrays 500000000 constitutive_titanmod.f90 + ifort -c -O3 -heap-arrays 500000000 constitutive_titanmod.f90 constitutive_nonlocal.o: constitutive_nonlocal.f90 basics.a advanced.a - ifort -c -heap-arrays 500000000 constitutive_nonlocal.f90 + ifort -c -O3 -heap-arrays 500000000 constitutive_nonlocal.f90 constitutive_dislotwin.o: constitutive_dislotwin.f90 basics.a advanced.a - ifort -c -heap-arrays 500000000 constitutive_dislotwin.f90 + ifort -c -O3 -heap-arrays 500000000 constitutive_dislotwin.f90 constitutive_j2.o: constitutive_j2.f90 basics.a advanced.a - ifort -c -heap-arrays 500000000 constitutive_j2.f90 + ifort -c -O3 -heap-arrays 500000000 constitutive_j2.f90 constitutive_phenopowerlaw.o: constitutive_phenopowerlaw.f90 basics.a advanced.a - ifort -c -heap-arrays 500000000 constitutive_phenopowerlaw.f90 + ifort -c -O3 -heap-arrays 500000000 constitutive_phenopowerlaw.f90 @@ -47,13 +47,13 @@ advanced.a: lattice.o ar rc advanced.a FEsolving.o mesh.o material.o lattice.o lattice.o: lattice.f90 material.o - ifort -c -heap-arrays 500000000 lattice.f90 + ifort -c -O3 -heap-arrays 500000000 lattice.f90 material.o: material.f90 mesh.o - ifort -c -heap-arrays 500000000 material.f90 + ifort -c -O3 -heap-arrays 500000000 material.f90 mesh.o: mesh.f90 FEsolving.o - ifort -c -heap-arrays 500000000 mesh.f90 + ifort -c -O3 -heap-arrays 500000000 mesh.f90 FEsolving.o: FEsolving.f90 basics.a - ifort -c -heap-arrays 500000000 FEsolving.f90 + ifort -c -O3 -heap-arrays 500000000 FEsolving.f90 @@ -61,16 +61,16 @@ basics.a: debug.o math.o ar rc basics.a debug.o math.o numerics.o IO.o mpie_spectral_interface.o prec.o debug.o: debug.f90 numerics.o - ifort -c debug.f90 + ifort -c -O3 debug.f90 math.o: math.f90 numerics.o - ifort -c math.f90 + ifort -c -O3 math.f90 numerics.o: numerics.f90 IO.o - ifort -c numerics.f90 + ifort -c -O3 numerics.f90 IO.o: IO.f90 mpie_spectral_interface.o - ifort -c IO.f90 + ifort -c -O3 IO.f90 mpie_spectral_interface.o: mpie_spectral_interface.f90 prec.o - ifort -c mpie_spectral_interface.f90 + ifort -c -O3 mpie_spectral_interface.f90 prec.o: prec.f90 - ifort -c prec.f90 + ifort -c -O3 prec.f90 diff --git a/code/mpie_spectral.f90 b/code/mpie_spectral.f90 index 8e598ad1d..daaf6d8e3 100644 --- a/code/mpie_spectral.f90 +++ b/code/mpie_spectral.f90 @@ -31,7 +31,7 @@ program mpie_spectral use CPFEM, only: CPFEM_general implicit none - include 'fftw3.f' !header file for fftw3 (declaring variables) library is also needed + include 'fftw3.f' !header file for fftw3 (declaring variables) library file is also needed !variables to read in from loadcase and mesh file real(pReal), dimension(9) :: valuevector ! stores information temporarily from loadcase file @@ -61,33 +61,34 @@ program mpie_spectral real(pReal), dimension(3) :: meshdimension ! stress etc. - real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation - real(pReal), dimension (3,3) :: pstress ! Piola-Kirchhoff stress in Matrix notation - real(pReal), dimension (3,3,3,3) :: dPdF,c0,s0 ! ??, reference stiffnes, (reference stiffness)^-1 - real(pReal), dimension(6,6) :: dsde,c066,s066 ! mandel notation - real(pReal), dimension(3,3) :: disgradmacro - real(pReal), dimension(3,3) :: cstress_av, defgrad_av, aux33 - real(pReal), dimension(:,:,:,:,:), allocatable :: cstress_field, defgrad, defgradold, start + real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation + real(pReal), dimension(3,3) :: pstress ! Piola-Kirchhoff stress in Matrix notation + real(pReal), dimension(3,3,3,3) :: dPdF,c0,s0 ! ??, reference stiffnes, (reference stiffness)^-1 + real(pReal), dimension(6,6) :: dsde,c066,s066 ! mandel notation + real(pReal), dimension(3,3) :: disgradmacro + real(pReal), dimension(3,3) :: cstress_av, defgrad_av, temp33_Real + real(pReal), dimension(:,:,:,:,:), allocatable :: cstress_field, defgrad, defgradold, ddefgrad ! variables storing information for spectral method - complex(pReal), dimension(:,:,:,:,:), allocatable :: workfft - complex(pReal), dimension(3,3) :: ddefgrad - real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat - real(pReal), dimension(3) :: xk - real(pReal), dimension(3,3) :: xknormdyad - integer(pInt), dimension(3) :: k_s - integer*8, dimension(2) :: plan + complex(pReal), dimension(:,:,:,:,:), allocatable :: workfft + complex(pReal), dimension(3,3) :: temp33_Complex + real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat + real(pReal), dimension(3) :: xk + real(pReal), dimension(3,3) :: xknormdyad + integer(pInt), dimension(3) :: k_s + integer*8, dimension(2,3,3) :: plan_fft ! convergency etc. logical errmatinv integer(pInt) itmax, ierr real(pReal) error, err_stress_av, err_stress_max, err_strain_av, err_strain_max real(pReal), dimension(3,3) :: strain_err, cstress_err - + ! loop variables etc. integer(pInt) i, j, k, l, m, n, p integer(pInt) loadcase, ielem, ial, iter, calcmode - + real(pReal), dimension(2) :: guessmode + real(pReal) temperature ! not used, but needed !gmsh @@ -123,7 +124,7 @@ program mpie_spectral if (IargC() < 2) call IO_error(102) ! check for correct Nr. of arguments given -! Reading the loadcase file and assingn variables +! Reading the loadcase file and assign variables path = getLoadcaseName() print*,'Loadcase: ',trim(path) print*,'Workingdir: ',trim(getSolverWorkingDirectoryName()) @@ -260,15 +261,21 @@ program mpie_spectral print *,'homogenization',homog print *, '' - allocate (workfft(resolution(1),resolution(2),resolution(3),3,3)); workfft = .0_pReal - allocate (gamma_hat(3,3,3,3,resolution(1),resolution(2),resolution(3))); gamma_hat = .0_pReal - allocate (cstress_field(resolution(1),resolution(2),resolution(3),3,3)); cstress_field = .0_pReal - allocate (defgrad(3,3,resolution(1),resolution(2),resolution(3))); defgrad = .0_pReal - allocate (defgradold(3,3,resolution(1),resolution(2),resolution(3))); defgradold = .0_pReal - allocate (start(3,3,resolution(1),resolution(2),resolution(3))); start = .0_pReal + allocate (workfft(resolution(1)/2+1,resolution(2),resolution(3),3,3)); workfft = .0_pReal + allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = .0_pReal + allocate (cstress_field(resolution(1),resolution(2),resolution(3),3,3)); cstress_field = .0_pReal + allocate (defgrad(resolution(1),resolution(2),resolution(3),3,3)); defgrad = .0_pReal + allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = .0_pReal + allocate (ddefgrad(resolution(1),resolution(2),resolution(3),3,3)); ddefgrad = .0_pReal - call dfftw_plan_dft_r2c_3d(plan(1),resolution(1),resolution(2),resolution(3), cstress_field(:,:,:,:,:),workfft(1:resolution(1)/2+1,:,:,:,:), FFTW_PATIENT) - call dfftw_plan_dft_3d(plan(2), resolution(1),resolution(2),resolution(3), workfft(:,:,:,:,:),workfft(:,:,:,:,:), FFTW_FORWARD, FFTW_PATIENT) + call dfftw_init_threads(ierr) + call dfftw_plan_with_nthreads(4) + do m = 1,3; do n = 1,3 + call dfftw_plan_dft_r2c_3d(plan_fft(1,m,n),resolution(1),resolution(2),resolution(3),& + cstress_field(:,:,:,m,n), workfft(:,:,:,m,n), FFTW_PATIENT, FFTW_DESTROY_INPUT) + call dfftw_plan_dft_c2r_3d(plan_fft(2,m,n),resolution(1),resolution(2),resolution(3),& + workfft(:,:,:,m,n), ddefgrad(:,:,:,m,n), FFTW_PATIENT) + enddo; enddo prodnn = resolution(1)*resolution(2)*resolution(3) wgt = 1._pReal/real(prodnn, pReal) @@ -276,11 +283,11 @@ program mpie_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) - defgradold(:,:,i,j,k) = math_I3 !to fit calculation of first step to calculation of following steps - defgrad(:,:,i,j,k) = math_I3 !to fit calculation of first step to calculation of following steps + defgradold(i,j,k,:,:) = math_I3 !to fit calculation of first step to calculation of following steps + defgrad(i,j,k,:,:) = math_I3 !to fit calculation of first step to calculation of following steps ielem = ielem +1 !loop over FPs and determine elastic constants of reference material - call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde, pstress, dPdF) - c066 = c066+dsde + call CPFEM_general(2,math_I3,math_I3,temperature,0.0_pReal,ielem,1_pInt,cstress,dsde,pstress,dPdF) + c066 = c066 + dsde enddo; enddo; enddo c066 = c066*wgt @@ -291,18 +298,18 @@ program mpie_spectral s0 = math_Mandel66to3333(s066) c0 = math_Mandel66to3333(c066) +!calculation of gamma_hat (only approx half of it) do k = 1, resolution(3) k_s(3) = k-1 - if(k > resolution(3)/2) k_s(3) = k_s(3)-resolution(3) + if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) xk(3) = .0_pReal if(resolution(3) > 1) xk(3) = real(k_s(3), pReal)/meshdimension(3) do j = 1, resolution(2) k_s(2) = j-1 - if(j > resolution(2)/2) k_s(2) = k_s(2)-resolution(2) + if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) xk(2) = real(k_s(2), pReal)/meshdimension(2) - do i = 1, resolution(1) + do i = 1, resolution(1)/2+1 k_s(1) = i-1 - if(i > resolution(1)/2) k_s(1) = k_s(1) -resolution(1) xk(1) = real(k_s(1), pReal)/meshdimension(1) xknormdyad=.0_pReal @@ -314,15 +321,15 @@ program mpie_spectral endif !forall loops don't work for the next 2 loop constructs!!! - aux33 = .0_pReal + temp33_Real = .0_pReal do l = 1,3; do m = 1,3; do n = 1,3; do p = 1,3 - aux33(l,m) = aux33(l,m)+c0(l,n,m,p)*xknormdyad(n,p) + temp33_Real(l,m) = temp33_Real(l,m)+c0(l,n,m,p)*xknormdyad(n,p) enddo; enddo; enddo; enddo - aux33 = math_inv3x3(aux33) + temp33_Real = math_inv3x3(temp33_Real) do l=1,3; do m=1,3; do n=1,3; do p=1,3 - gamma_hat(l,m,n,p,i,j,k) = -aux33(l,n)*xknormdyad(m,p) + gamma_hat(i,j,k,l,m,n,p) = -temp33_Real(l,n)*xknormdyad(m,p) enddo; enddo; enddo; enddo enddo; enddo; enddo @@ -340,9 +347,8 @@ program mpie_spectral timeinc = bc_timeIncrement(loadcase)/bc_steps(loadcase) - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) !no fluctuation as guess for new loadcase (last two summands will disappear) - start(:,:,i,j,k) = bc_velocityGrad(:,:,loadcase)*timeinc -defgrad(:,:,i,j,k) + defgradold(:,:,i,j,k) - enddo; enddo; enddo + guessmode(1) =.0_pReal + guessmode(2) = 1_pReal !************************************************************* ! loop oper steps defined in input file for current loadcase @@ -352,13 +358,15 @@ program mpie_spectral write(*,*) 'STEP = ',steps do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) - aux33 = defgrad(:,:,i,j,k) - defgrad(:,:,i,j,k) = 2 * defgrad(:,:,i,j,k) - defgradold(:,:,i,j,k) + start(:,:,i,j,k) ! old fluctuations as guess - defgradold(:,:,i,j,k) = aux33 ! wind forward + temp33_Real = defgrad(i,j,k,:,:) + defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) + guessmode(1)*(defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& ! old fluctuations as guess for new step + + guessmode(2)*bc_velocityGrad(:,:,loadcase)*timeinc ! homogeneous fluctuations for new loadcase + defgradold(i,j,k,:,:) = temp33_Real ! wind forward enddo; enddo; enddo - disgradmacro = disgradmacro + bc_velocityGrad(:,:,loadcase)*timeinc !update macroscopic displacementgradient (stores the desired BCs of defgrad) - start = .0_pReal + disgradmacro = disgradmacro + bc_velocityGrad(:,:,loadcase)*timeinc !update macroscopic displacementgradient (stores the desired BCs of defgrad) + guessmode(1)= 1_pReal + guessmode(2)=.0_pReal calcmode = 1_pInt iter = 0_pInt err_stress_av = 2.*error; err_strain_av = 2.*error @@ -377,7 +385,7 @@ program mpie_spectral ielem = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1 - call CPFEM_general(3, defgradold(:,:,i,j,k), defgrad(:,:,i,j,k),& + call CPFEM_general(3, defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo @@ -387,22 +395,22 @@ program mpie_spectral do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1 call CPFEM_general(calcmode,& ! first element in first iteration retains calcMode 1, others get 2 (saves winding forward effort) - defgradold(:,:,i,j,k), defgrad(:,:,i,j,k),& + defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& temperature,timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) calcmode = 2 - aux33 = math_Mandel6to33(cstress) + temp33_Real = math_Mandel6to33(cstress) do m = 1,3; do n = 1,3 ! calculate stress error - if(abs(aux33(m,n)) > 0.1 * abs(cstress_err(m,n))) then ! only stress components larger than 10% are taking under consideration - err_stress_av = err_stress_av + abs((cstress_field(i,j,k,m,n)-aux33(m,n))/aux33(m,n)) !any find maxval> leave loop - err_stress_max = max(err_stress_max, abs((cstress_field(i,j,k,m,n)-aux33(m,n))/aux33(m,n))) + if(abs(temp33_Real(m,n)) > 0.1 * abs(cstress_err(m,n))) then ! only stress components larger than 10% are taking under consideration + err_stress_av = err_stress_av + abs((cstress_field(i,j,k,m,n)-temp33_Real(m,n))/temp33_Real(m,n)) !any find maxval> leave loop + err_stress_max = max(err_stress_max, abs((cstress_field(i,j,k,m,n)-temp33_Real(m,n))/temp33_Real(m,n))) l=l+1 endif enddo; enddo - cstress_field(i,j,k,:,:) = aux33 - cstress_av = cstress_av + aux33 ! average stress + cstress_field(i,j,k,:,:) = temp33_Real + cstress_av = cstress_av + temp33_Real ! average stress enddo; enddo; enddo err_stress_av = err_stress_av/l ! do the weighting of the error @@ -411,42 +419,30 @@ program mpie_spectral write(*,*) 'SPECTRAL METHOD TO GET CHANGE OF DEFORMATION GRADIENT FIELD' do m = 1,3; do n = 1,3 - call dfftw_execute_dft_r2c(plan(1), cstress_field(:,:,:,m,n),workfft(1:resolution(1)/2+1,:,:,m,n)) + call dfftw_execute_dft_r2c(plan_fft(1,m,n), cstress_field(:,:,:,m,n),workfft(:,:,:,m,n)) enddo; enddo - workfft=conjg(workfft) - do i = 0, resolution(1)/2-2 !unpack fft data for conj complex symmetric part. can be directly used in calculation of cstress_field - m = 1 - do k = 1, resolution(3) - n = 1 - do j = 1, resolution(2) - workfft(resolution(1)-i,j,k,:,:) = conjg(workfft(2+i,n,m,:,:)) - if(n == 1) n = resolution(2) +1 - n = n-1 - enddo - if(m == 1) m = resolution(3) +1 - m = m -1 - enddo; enddo - - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) - ddefgrad = .0_pReal + + do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 !no better way to do the calculation??? + temp33_Complex = .0_pReal do m = 1,3; do n = 1,3 - ddefgrad(m,n) = ddefgrad(m,n) +sum(gamma_hat(m,n,:,:,i,j,k)*workfft(i,j,k,:,:)) + temp33_Complex(m,n) = temp33_Complex(m,n) +sum(gamma_hat(i,j,k,m,n,:,:) * workfft(i,j,k,:,:)) enddo; enddo - workfft(i,j,k,:,:) = ddefgrad(:,:) + workfft(i,j,k,:,:) = temp33_Complex(:,:) enddo; enddo; enddo do m = 1,3; do n = 1,3 - call dfftw_execute_dft(plan(2), workfft(:,:,:,m,n), workfft(:,:,:,m,n)) - defgrad(m,n,:,:,:) = defgrad(m,n,:,:,:) + real(workfft(:,:,:,m,n), pReal)*wgt + call dfftw_execute_dft_c2r(plan_fft(2,m,n), workfft(:,:,:,m,n),ddefgrad(:,:,:,m,n)) enddo; enddo + + defgrad = defgrad + ddefgrad * wgt l = 0_pInt do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) - defgrad_av(:,:) = defgrad_av(:,:) + defgrad(:,:,i,j,k) ! calculate average strain + defgrad_av(:,:) = defgrad_av(:,:) + defgrad(i,j,k,:,:) ! calculate average strain do m = 1,3; do n = 1,3 ! calculate strain error - if(abs(defgrad(m,n,i,j,k)) > 0.1 * abs(strain_err(m,n))) then - err_strain_av = err_strain_av + abs((real(workfft(i,j,k,m,n), pReal)*wgt)/defgrad(m,n,i,j,k)) - err_strain_max = max(err_strain_max, abs((real(workfft(i,j,k,m,n), pReal)*wgt)/defgrad(m,n,i,j,k))) + if(abs(defgrad(i,j,k,m,n)) > 0.1 * abs(strain_err(m,n))) then + err_strain_av = err_strain_av + abs((real(ddefgrad(i,j,k,m,n), pReal)*wgt)/defgrad(i,j,k,m,n)) + err_strain_max = max(err_strain_max, abs((real(ddefgrad(i,j,k,m,n), pReal)*wgt)/defgrad(i,j,k,m,n))) l=l+1 endif enddo; enddo @@ -458,17 +454,17 @@ program mpie_spectral do m = 1,3; do n = 1,3 if(bc_mask(m,n,1,loadcase)) then !adjust defgrad to achieve displacement BC (disgradmacro) - defgrad(m,n,:,:,:) = defgrad(m,n,:,:,:) + (disgradmacro(m,n)+math_I3(m,n)-defgrad_av(m,n)) + defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (disgradmacro(m,n)+math_I3(m,n)-defgrad_av(m,n)) endif if(bc_mask(m,n,2,loadcase)) then !adjust defgrad to achieve convergency in stress - defgrad(m,n,:,:,:) = defgrad(m,n,:,:,:) + sum(s0(m,n,:,:)*bc_stress_i(:,:)*(bc_stress(:,:,loadcase)-cstress_av(:,:))) + defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + sum(s0(m,n,:,:)*bc_stress_i(:,:)*(bc_stress(:,:,loadcase)-cstress_av(:,:))) endif enddo; enddo - write(*,*) 'STRESS FIELD ERROR AV = ',err_strain_av - write(*,*) 'STRAIN FIELD ERROR AV = ',err_stress_av - write(*,*) 'STRESS FIELD ERROR MAX = ',err_strain_max - write(*,*) 'STRAIN FIELD ERROR MAX = ',err_stress_max + write(*,*) 'STRESS FIELD ERROR AV = ',err_stress_av + write(*,*) 'STRAIN FIELD ERROR AV = ',err_strain_av + write(*,*) 'STRESS FIELD ERROR MAX = ',err_stress_max + write(*,*) 'STRAIN FIELD ERROR MAX = ',err_strain_max enddo ! end looping when convergency is achieved write(539,'(f12.6,a,f12.6)'),defgrad_av(3,3)-1,' ',cstress_av(3,3) @@ -482,7 +478,7 @@ program mpie_spectral !gsmh output write(nriter, *) iter write(nrstep, *) steps - nrstep='defgrad'//trim(adjustl(nrstep))//trim(adjustl(nriter))//'_cpfem.msh' + nrstep='defgrad'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh' open(589,file=nrstep) write(589, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn do i = 1, prodnn @@ -500,14 +496,14 @@ program mpie_spectral do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) ielem = ielem + 1 write(589, '(i10,f16.8,tr2,f16.8,tr2,f16.8,tr2,f16.8,tr2,f16.8,tr2,f16.8,& - tr2,f16.8,tr2,f16.8,tr2,f16.8,tr2)'), ielem, defgrad(:,:,i,j,k) + tr2,f16.8,tr2,f16.8,tr2,f16.8,tr2)'), ielem, defgrad(i,j,k,:,:) enddo; enddo; enddo write(589, *), '$EndNodeData' close(589) write(nriter, *) iter write(nrstep, *) steps - nrstep = 'stress'//trim(adjustl(nrstep))//trim(adjustl(nriter))//'_cpfem.msh' + nrstep = 'stress'//trim(adjustl(nrstep))//'-'//trim(adjustl(nriter))//'_cpfem.msh' open(589,file = nrstep) write(589, '(A, /, A, /, A, /, A, /, I10)'), '$MeshFormat', '2.1 0 8', '$EndMeshFormat', '$Nodes', prodnn do i = 1, prodnn @@ -534,8 +530,11 @@ program mpie_spectral enddo ! end loping over steps in current loadcase enddo ! end looping over loadcases close(539) -call dfftw_destroy_plan(plan(2)) -call dfftw_destroy_plan(plan(1)) + +do i=1,2; do m = 1,3; do n = 1,3 + call dfftw_destroy_plan(plan_fft(i,m,n)) +enddo; enddo; enddo + end program mpie_spectral !********************************************************************