2011-08-26 19:36:37 +05:30
! Copyright 2011 Max-Planck-Institut fuer Eisenforschung GmbH
2011-04-04 19:39:54 +05:30
!
! This file is part of DAMASK,
2011-11-04 01:02:11 +05:30
! the Duesseldorf Advanced Material Simulation Kit.
2011-04-04 19:39:54 +05:30
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
2010-06-08 15:40:57 +05:30
!* $Id$
2010-06-08 15:38:15 +05:30
!********************************************************************
! Material subroutine for BVP solution using spectral method
!
! written by P. Eisenlohr,
! F. Roters,
! L. Hantcherli,
2011-02-07 20:05:42 +05:30
! W.A. Counts,
! D.D. Tjahjanto,
! C. Kords,
! M. Diehl,
2010-06-08 15:38:15 +05:30
! R. Lebensohn
!
! MPI fuer Eisenforschung, Duesseldorf
!
!********************************************************************
! Usage:
2011-10-18 20:15:32 +05:30
! - start program with DAMASK_spectral
! -g (--geom, --geometry) PathToGeomFile/NameOfGeom.geom
! -l (--load, --loadcase) PathToLoadFile/NameOfLoadFile.load
2011-02-07 20:05:42 +05:30
! - PathToGeomFile will be the working directory
2010-06-08 15:38:15 +05:30
! - make sure the file "material.config" exists in the working
2011-11-15 23:24:18 +05:30
! directory. For further configuration use "numerics.config" and
! "numerics.config"
2010-06-10 20:21:10 +05:30
!********************************************************************
2011-05-11 22:08:45 +05:30
program DAMASK_spectral
2010-06-10 20:21:10 +05:30
!********************************************************************
2010-06-08 15:38:15 +05:30
2011-05-11 22:31:03 +05:30
use DAMASK_interface
2011-12-01 17:31:13 +05:30
use prec , only : pInt , pReal , DAMASK_NaN
2010-06-10 20:21:10 +05:30
use IO
2011-11-21 23:42:40 +05:30
use debug , only : spectral_debug_verbosity
2010-07-05 21:31:36 +05:30
use math
2011-05-26 14:53:13 +05:30
use mesh , only : mesh_ipCenterOfGravity
2011-11-11 19:47:43 +05:30
use CPFEM , only : CPFEM_general , CPFEM_initAll
use FEsolving , only : restartWrite , restartReadSpectral , restartReadStep
2011-11-18 03:41:05 +05:30
use numerics , only : err_div_tol , err_stress_tolrel , rotation_tol , &
2011-11-15 23:24:18 +05:30
itmax , memory_efficient , DAMASK_NumThreadsInt , divergence_correction , &
2011-10-25 19:08:24 +05:30
fftw_planner_flag , fftw_timelimit
2010-10-27 22:45:49 +05:30
use homogenization , only : materialpoint_sizeResults , materialpoint_results
2011-11-21 23:42:40 +05:30
!$ use OMP_LIB ! the openMP function library
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
2010-06-08 15:38:15 +05:30
implicit none
2011-01-07 18:26:45 +05:30
! variables to read from loadcase and geom file
2011-11-21 23:42:40 +05:30
real ( pReal ) , dimension ( 9 ) :: temp_valueVector ! stores information temporarily from loadcase file
2011-11-15 23:24:18 +05:30
logical , dimension ( 9 ) :: temp_maskVector
2011-10-18 20:15:32 +05:30
integer ( pInt ) , parameter :: maxNchunksLoadcase = &
2011-10-24 23:56:34 +05:30
( 1_pInt + 9_pInt ) * 3_pInt + & ! deformation, rotation, and stress
2011-11-07 16:34:57 +05:30
( 1_pInt + 1_pInt ) * 5_pInt + & ! time, (log)incs, temp, restartfrequency, and outputfrequency
2011-10-24 23:56:34 +05:30
1_pInt ! dropguessing
2011-11-21 23:42:40 +05:30
integer ( pInt ) , dimension ( 1_pInt + maxNchunksLoadcase * 2_pInt ) :: posLoadcase
2011-10-25 19:08:24 +05:30
integer ( pInt ) , parameter :: maxNchunksGeom = 7_pInt ! 4 identifiers, 3 values
2011-11-21 23:42:40 +05:30
integer ( pInt ) , dimension ( 1_pInt + maxNchunksGeom * 2_pInt ) :: posGeom
2011-11-15 23:24:18 +05:30
integer ( pInt ) :: headerLength , N_l = 0_pInt , N_t = 0_pInt , N_n = 0_pInt , N_Fdot = 0_pInt
2011-11-07 23:55:10 +05:30
integer ( pInt ) , parameter :: myUnit = 234_pInt
2011-10-18 20:15:32 +05:30
character ( len = 1024 ) :: path , line , keyword
2011-11-07 23:55:10 +05:30
logical :: gotResolution = . false . , gotDimension = . false . , gotHomogenization = . false .
2011-11-15 23:24:18 +05:30
type bc_type
real ( pReal ) , dimension ( 3 , 3 ) :: deformation , & ! applied velocity gradient or time derivative of deformation gradient
stress , & ! stress BC (if applicable)
rotation ! rotation of BC (if applicable)
real ( pReal ) :: timeIncrement , & ! length of increment
temperature ! isothermal starting conditions
integer ( pInt ) :: steps , & ! number of steps
outputfrequency , & ! frequency of result writes
restartfrequency , & ! frequency of result writes
logscale ! linear/logaritmic time step flag
logical :: followFormerTrajectory , & ! follow trajectory of former loadcase
velGradApplied ! decide wether velocity gradient or fdot is given
logical , dimension ( 3 , 3 ) :: maskDeformation , & ! mask of boundary conditions
maskStress
logical , dimension ( 9 ) :: maskStressVector ! linear mask of boundary conditions
end type
2011-12-04 15:31:32 +05:30
type ( bc_type ) , allocatable , dimension ( : ) :: bc
type ( bc_type ) :: bc_init
character ( len = 3 ) :: loadcase_string
2010-07-01 20:50:06 +05:30
2011-01-07 18:26:45 +05:30
! variables storing information from geom file
2011-12-04 15:31:32 +05:30
real ( pReal ) :: wgt
real ( pReal ) , dimension ( 3 ) :: geomdimension = 0.0_pReal ! physical dimension of volume element per direction
integer ( pInt ) :: Npoints , & ! number of Fourier points
homog ! homogenization scheme used
integer ( pInt ) , dimension ( 3 ) :: res = 1_pInt ! resolution (number of Fourier points) in each direction
logical :: spectralPictureMode = . false . ! indicating 1 to 1 mapping of FP to microstructure
2010-07-01 20:50:06 +05:30
2011-11-07 23:55:10 +05:30
! stress, stiffness and compliance average etc.
2011-12-04 15:31:32 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: pstress , pstress_av , defgrad_av , &
defgradAim = math_I3 , defgradAimOld = math_I3 , defgradAimCorr = math_I3 , &
mask_stress , mask_defgrad , fDot , &
pstress_av_load , defgradAim_lab ! quantities rotated to other coordinate system
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dPdF , c0_reference , c_current = 0.0_pReal , s_prev , c_prev ! stiffness and compliance
real ( pReal ) , dimension ( 6 ) :: cstress ! cauchy stress
real ( pReal ) , dimension ( 6 , 6 ) :: dsde ! small strain stiffness
real ( pReal ) , dimension ( 9 , 9 ) :: s_prev99 , c_prev99 ! compliance and stiffness in matrix notation
real ( pReal ) , dimension ( : , : ) , allocatable :: s_reduced , c_reduced ! reduced compliance and stiffness (only for stress BC)
integer ( pInt ) :: size_reduced = 0.0_pReal ! number of stress BCs
2011-10-18 20:15:32 +05:30
! pointwise data
2011-11-15 23:24:18 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: workfft , defgrad , defgradold
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: coordinates
real ( pReal ) , dimension ( : , : , : ) , allocatable :: temperature
2011-10-18 20:15:32 +05:30
2011-10-25 19:08:24 +05:30
! variables storing information for spectral method and FFTW
2011-08-26 19:36:37 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: xiDyad ! product of wave vectors
real ( pReal ) , dimension ( : , : , : , : , : , : , : ) , allocatable :: gamma_hat ! gamma operator (field) for spectral method
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: xi ! wave vector field
integer ( pInt ) , dimension ( 3 ) :: k_s
2011-11-15 23:24:18 +05:30
integer * 8 , dimension ( 3 ) :: fftw_plan ! plans for fftw (forward and backward)
2011-02-09 23:17:28 +05:30
2011-02-07 20:05:42 +05:30
! loop variables, convergence etc.
2011-11-11 19:47:43 +05:30
real ( pReal ) :: time = 0.0_pReal , time0 = 0.0_pReal , timeinc ! elapsed time, begin of interval, time interval
2011-11-18 03:41:05 +05:30
real ( pReal ) :: guessmode , err_div , err_stress , err_stress_tol , p_hat_avg
2011-11-15 23:24:18 +05:30
complex ( pReal ) , parameter :: img = cmplx ( 0.0_pReal , 1.0_pReal )
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: ones = 1.0_pReal , zeroes = 0.0_pReal
complex ( pReal ) , dimension ( 3 , 3 ) :: temp33_Complex
real ( pReal ) , dimension ( 3 , 3 ) :: temp33_Real
2011-12-04 15:31:32 +05:30
integer ( pInt ) :: i , j , k , l , m , n , p , errorID
2011-11-18 03:41:05 +05:30
integer ( pInt ) :: N_Loadcases , loadcase , step , iter , ielem , CPFEM_mode , &
ierr , notConvergedCounter = 0_pInt , totalStepsCounter = 0_pInt
2011-11-15 23:24:18 +05:30
logical :: errmatinv , regrid = . false .
2011-11-11 19:47:43 +05:30
real ( pReal ) :: defgradDet , defgradDetMax , defgradDetMin
2011-11-15 23:24:18 +05:30
real ( pReal ) :: correctionFactor
2011-12-04 15:31:32 +05:30
! --- debugging variables
2011-11-21 23:42:40 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: divergence
2011-11-15 23:24:18 +05:30
real ( pReal ) :: p_real_avg , err_div_max , err_real_div_avg , err_real_div_max
logical :: debugGeneral = . false . , debugDivergence = . false . , debugRestart = . false .
2011-07-25 22:00:21 +05:30
2011-12-04 15:31:32 +05:30
! --- initialize default value for loadcase
2011-11-21 23:42:40 +05:30
bc_init % deformation = zeroes ; bc_init % stress = zeroes ; bc_init % rotation = zeroes
2011-12-04 15:31:32 +05:30
bc_init % timeIncrement = 0.0_pReal ; bc_init % temperature = 30 0.0_pReal
bc_init % steps = 0_pInt ; bc_init % logscale = 0_pInt
bc_init % outputfrequency = 1_pInt ; bc_init % restartfrequency = 1_pInt
bc_init % maskDeformation = . false . ; bc_init % maskStress = . false .
bc_init % maskStressVector = . false . ; bc_init % velGradApplied = . false .
2011-11-21 23:42:40 +05:30
bc_init % followFormerTrajectory = . true .
bc_init % rotation = math_I3 ! assume no rotation
2011-12-04 15:31:32 +05:30
! --- initializing model size independed parameters
2011-11-18 03:41:05 +05:30
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
if ( . not . ( command_argument_count ( ) == 4 . or . command_argument_count ( ) == 6 ) ) & ! check for correct number of given arguments
call IO_error ( error_ID = 102_pInt )
2011-08-26 19:36:37 +05:30
2011-11-04 01:02:11 +05:30
call DAMASK_interface_init ( )
2011-08-02 19:28:28 +05:30
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-11-15 23:24:18 +05:30
print '(a,a)' , ' <<<+- DAMASK_spectral init -+>>>'
print '(a,a)' , ' $Id$'
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-11-15 23:24:18 +05:30
print '(a,a)' , ' Working Directory: ' , trim ( getSolverWorkingDirectoryName ( ) )
print '(a,a)' , ' Solver Job Name: ' , trim ( getSolverJobName ( ) )
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-11-18 03:41:05 +05:30
! Reading the loadcase file and allocate variables for loadcases
2011-11-04 01:02:11 +05:30
path = getLoadcaseName ( )
2011-12-04 15:31:32 +05:30
if ( . not . IO_open_file ( myUnit , path ) ) call IO_error ( error_ID = 30_pInt , ext_msg = trim ( path ) )
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
2010-06-10 20:21:10 +05:30
do
2011-10-18 20:15:32 +05:30
read ( myUnit , '(a1024)' , END = 100 ) line
2011-11-04 01:02:11 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2011-10-18 20:15:32 +05:30
posLoadcase = IO_stringPos ( line , maxNchunksLoadcase )
2011-11-15 23:24:18 +05:30
do i = 1_pInt , maxNchunksLoadcase , 1_pInt ! reading compulsory parameters for loadcase
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posLoadcase , i ) ) )
2011-10-25 19:08:24 +05:30
case ( 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' )
2011-11-15 23:24:18 +05:30
N_l = N_l + 1_pInt
2011-07-13 22:03:12 +05:30
case ( 'fdot' )
2011-11-15 23:24:18 +05:30
N_Fdot = N_Fdot + 1_pInt
2011-07-13 22:03:12 +05:30
case ( 't' , 'time' , 'delta' )
2011-11-15 23:24:18 +05:30
N_t = N_t + 1_pInt
2011-07-13 22:03:12 +05:30
case ( 'n' , 'incs' , 'increments' , 'steps' , 'logincs' , 'logsteps' )
2011-11-15 23:24:18 +05:30
N_n = N_n + 1_pInt
2010-06-10 20:21:10 +05:30
end select
2011-11-04 01:02:11 +05:30
enddo ! count all identifiers to allocate memory and do sanity check
2010-06-10 20:21:10 +05:30
enddo
2010-06-10 14:20:04 +05:30
2011-10-18 20:15:32 +05:30
100 N_Loadcases = N_n
2011-11-04 01:02:11 +05:30
if ( ( N_l + N_Fdot / = N_n ) . or . ( N_n / = N_t ) ) & ! sanity check
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 37_pInt , ext_msg = trim ( path ) ) ! error message for incomplete loadcase
allocate ( bc ( N_Loadcases ) )
2010-06-10 20:21:10 +05:30
2011-12-04 15:31:32 +05:30
! --- reading the loadcase and assign values to the allocated data structure
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
2011-07-14 15:07:31 +05:30
loadcase = 0_pInt
2010-06-10 20:21:10 +05:30
do
2011-10-18 20:15:32 +05:30
read ( myUnit , '(a1024)' , END = 101 ) line
2011-11-15 23:24:18 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
loadcase = loadcase + 1_pInt
2011-11-21 23:42:40 +05:30
bc ( loadcase ) = bc_init
2011-10-18 20:15:32 +05:30
posLoadcase = IO_stringPos ( line , maxNchunksLoadcase )
2011-11-15 23:24:18 +05:30
do j = 1_pInt , maxNchunksLoadcase
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) )
2011-11-21 23:42:40 +05:30
case ( 'fdot' , 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' ) ! assign values for the deformation BC matrix
bc ( loadcase ) % velGradApplied = ( IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'l' . or . & ! in case of given L, set flag to true
2011-11-15 23:24:18 +05:30
IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'velocitygrad' . or . &
IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'velgrad' . or . &
IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'velocitygradient' )
temp_valueVector = 0.0_pReal
temp_maskVector = . false .
forall ( k = 1_pInt : 9_pInt ) temp_maskVector ( k ) = IO_stringValue ( line , posLoadcase , j + k ) / = '*'
do k = 1_pInt , 9_pInt
if ( temp_maskVector ( k ) ) temp_valueVector ( k ) = IO_floatValue ( line , posLoadcase , j + k )
2010-06-10 20:21:10 +05:30
enddo
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % maskDeformation = transpose ( reshape ( temp_maskVector , ( / 3 , 3 / ) ) )
bc ( loadcase ) % deformation = math_plain9to33 ( temp_valueVector )
2011-11-04 01:02:11 +05:30
case ( 'p' , 'pk1' , 'piolakirchhoff' , 'stress' )
2011-11-15 23:24:18 +05:30
temp_valueVector = 0.0_pReal
forall ( k = 1_pInt : 9_pInt ) bc ( loadcase ) % maskStressVector ( k ) = IO_stringValue ( line , posLoadcase , j + k ) / = '*'
do k = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( k ) ) temp_valueVector ( k ) = IO_floatValue ( line , posLoadcase , j + k ) ! assign values for the bc(loadcase)%stress matrix
2010-06-10 20:21:10 +05:30
enddo
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % maskStress = transpose ( reshape ( bc ( loadcase ) % maskStressVector , ( / 3 , 3 / ) ) )
bc ( loadcase ) % stress = math_plain9to33 ( temp_valueVector )
case ( 't' , 'time' , 'delta' ) ! increment time
bc ( loadcase ) % timeIncrement = IO_floatValue ( line , posLoadcase , j + 1_pInt )
case ( 'temp' , 'temperature' ) ! starting temperature
bc ( loadcase ) % temperature = IO_floatValue ( line , posLoadcase , j + 1_pInt )
case ( 'n' , 'incs' , 'increments' , 'steps' ) ! steps
bc ( loadcase ) % steps = IO_intValue ( line , posLoadcase , j + 1_pInt )
case ( 'logincs' , 'logsteps' ) ! true, if log scale
bc ( loadcase ) % steps = IO_intValue ( line , posLoadcase , j + 1_pInt )
bc ( loadcase ) % logscale = 1_pInt
case ( 'f' , 'freq' , 'frequency' , 'outputfreq' ) ! frequency of result writings
bc ( loadcase ) % outputfrequency = IO_intValue ( line , posLoadcase , j + 1_pInt )
case ( 'r' , 'restart' , 'restartwrite' ) ! frequency of writing restart information
bc ( loadcase ) % restartfrequency = IO_intValue ( line , posLoadcase , j + 1_pInt )
2011-07-14 15:07:31 +05:30
case ( 'guessreset' , 'dropguessing' )
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % followFormerTrajectory = . false . ! do not continue to predict deformation along former trajectory
case ( 'euler' ) ! rotation of loadcase given in euler angles
p = 0_pInt ! assuming values given in radians
l = 1_pInt ! assuming keyword indicating degree/radians
select case ( IO_lc ( IO_stringValue ( line , posLoadcase , j + 1_pInt ) ) )
2011-10-18 20:15:32 +05:30
case ( 'deg' , 'degree' )
2011-11-15 23:24:18 +05:30
p = 1_pInt ! for conversion from degree to radian
2011-10-18 20:15:32 +05:30
case ( 'rad' , 'radian' )
case default
2011-11-15 23:24:18 +05:30
l = 0_pInt ! imediately reading in angles, assuming radians
2011-10-18 20:15:32 +05:30
end select
2011-11-15 23:24:18 +05:30
forall ( k = 1_pInt : 3_pInt ) temp33_Real ( k , 1 ) = IO_floatValue ( line , posLoadcase , j + l + k ) * real ( p , pReal ) * inRad
bc ( loadcase ) % rotation = math_EulerToR ( temp33_Real ( : , 1 ) )
case ( 'rotation' , 'rot' ) ! assign values for the rotation of loadcase matrix
temp_valueVector = 0.0_pReal
forall ( k = 1_pInt : 9_pInt ) temp_valueVector ( k ) = IO_floatValue ( line , posLoadcase , j + k )
bc ( loadcase ) % rotation = math_plain9to33 ( temp_valueVector )
2010-06-10 20:21:10 +05:30
end select
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
enddo ; enddo
2011-10-18 20:15:32 +05:30
101 close ( myUnit )
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
2011-12-04 15:31:32 +05:30
! --- read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90
2011-02-21 20:07:38 +05:30
path = getModelName ( )
2011-11-18 03:41:05 +05:30
2011-10-18 20:15:32 +05:30
if ( . not . IO_open_file ( myUnit , trim ( path ) / / InputFileExtension ) ) &
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 101_pInt , ext_msg = trim ( path ) / / InputFileExtension )
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
read ( myUnit , '(a1024)' ) line
2011-11-15 23:24:18 +05:30
posGeom = IO_stringPos ( line , 2_pInt )
keyword = IO_lc ( IO_StringValue ( line , posGeom , 2_pInt ) )
2011-10-18 20:15:32 +05:30
if ( keyword ( 1 : 4 ) == 'head' ) then
2011-11-15 23:24:18 +05:30
headerLength = IO_intValue ( line , posGeom , 1_pInt ) + 1_pInt
2011-10-18 20:15:32 +05:30
else
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 42_pInt )
2011-10-18 20:15:32 +05:30
endif
rewind ( myUnit )
2011-11-15 23:24:18 +05:30
do i = 1_pInt , headerLength
2011-10-18 20:15:32 +05:30
read ( myUnit , '(a1024)' ) line
2011-01-07 18:26:45 +05:30
posGeom = IO_stringPos ( line , maxNchunksGeom )
select case ( IO_lc ( IO_StringValue ( line , posGeom , 1 ) ) )
2010-06-25 17:01:05 +05:30
case ( 'dimension' )
2010-09-22 17:34:43 +05:30
gotDimension = . true .
2011-11-15 23:24:18 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posGeom , j ) ) )
2010-09-22 17:34:43 +05:30
case ( 'x' )
2011-11-15 23:24:18 +05:30
geomdimension ( 1 ) = IO_floatValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'y' )
2011-11-15 23:24:18 +05:30
geomdimension ( 2 ) = IO_floatValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'z' )
2011-11-15 23:24:18 +05:30
geomdimension ( 3 ) = IO_floatValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
end select
enddo
2010-06-25 17:01:05 +05:30
case ( 'homogenization' )
2010-09-22 17:34:43 +05:30
gotHomogenization = . true .
2011-11-15 23:24:18 +05:30
homog = IO_intValue ( line , posGeom , 2_pInt )
2010-06-25 17:01:05 +05:30
case ( 'resolution' )
2010-09-22 17:34:43 +05:30
gotResolution = . true .
2011-11-15 23:24:18 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posGeom , j ) ) )
2010-09-22 17:34:43 +05:30
case ( 'a' )
2011-11-21 23:42:40 +05:30
res ( 1 ) = IO_intValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'b' )
2011-11-21 23:42:40 +05:30
res ( 2 ) = IO_intValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'c' )
2011-11-21 23:42:40 +05:30
res ( 3 ) = IO_intValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
end select
enddo
2011-10-18 20:15:32 +05:30
case ( 'picture' )
spectralPictureMode = . true .
2010-06-25 17:01:05 +05:30
end select
enddo
2011-10-18 20:15:32 +05:30
close ( myUnit )
2011-12-04 15:31:32 +05:30
if ( . not . ( gotDimension . and . gotHomogenization . and . gotResolution ) ) call IO_error ( error_ID = 45_pInt )
2011-02-14 22:51:31 +05:30
2011-11-21 23:42:40 +05:30
if ( mod ( res ( 1 ) , 2_pInt ) / = 0_pInt . or . &
mod ( res ( 2 ) , 2_pInt ) / = 0_pInt . or . &
2011-12-04 15:31:32 +05:30
( mod ( res ( 3 ) , 2_pInt ) / = 0_pInt . and . res ( 3 ) / = 1_pInt ) ) call IO_error ( error_ID = 103_pInt )
2011-08-26 19:36:37 +05:30
2011-12-04 15:31:32 +05:30
Npoints = res ( 1 ) * res ( 2 ) * res ( 3 )
! --- initialization of CPFEM_general (= constitutive law)
2011-11-15 23:24:18 +05:30
call CPFEM_initAll ( bc ( 1 ) % temperature , 1_pInt , 1_pInt )
2011-11-11 19:47:43 +05:30
2011-12-04 15:31:32 +05:30
! --- debugging parameters
2011-11-18 03:41:05 +05:30
if ( iand ( spectral_debug_verbosity , 1_pInt ) == 1_pInt ) debugGeneral = . true .
if ( iand ( spectral_debug_verbosity , 2_pInt ) == 2_pInt ) debugDivergence = . true .
if ( iand ( spectral_debug_verbosity , 4_pInt ) == 4_pInt ) debugRestart = . true .
2011-11-21 23:42:40 +05:30
2011-12-04 15:31:32 +05:30
! --- output of geometry
print '(a)' , ''
print '(a)' , '#############################################################'
print '(a)' , 'DAMASK spectral:'
print '(a)' , 'The spectral method boundary value problem solver for'
print '(a)' , 'the Duesseldorf Advanced Material Simulation Kit'
print '(a)' , '#############################################################'
print '(a,a)' , 'geometry file: ' , trim ( path ) / / '.geom'
print '(a)' , '============================================================='
2011-11-21 23:42:40 +05:30
print '(a,i12,i12,i12)' , 'resolution a b c:' , res
2011-11-15 23:24:18 +05:30
print '(a,f12.5,f12.5,f12.5)' , 'dimension x y z:' , geomdimension
2011-12-04 15:31:32 +05:30
print '(a,i5)' , 'homogenization: ' , homog
print '(a,L)' , 'spectralPictureMode: ' , spectralPictureMode
print '(a)' , '#############################################################'
print '(a,a)' , 'loadcase file: ' , trim ( getLoadcaseName ( ) )
2011-11-15 23:24:18 +05:30
if ( bc ( 1 ) % followFormerTrajectory ) then
2011-12-04 15:31:32 +05:30
call IO_warning ( warning_ID = 33_pInt ) ! cannot guess along trajectory for first step of first loadcase
2011-11-15 23:24:18 +05:30
bc ( 1 ) % followFormerTrajectory = . false .
2011-10-18 20:15:32 +05:30
endif
2011-11-18 03:41:05 +05:30
2011-12-04 15:31:32 +05:30
! --- consistency checks and output of loadcase
errorID = 0_pInt
2011-11-11 19:47:43 +05:30
do loadcase = 1_pInt , N_Loadcases
2011-11-02 20:08:42 +05:30
write ( loadcase_string , '(i3)' ) loadcase
2011-12-04 15:31:32 +05:30
print '(a)' , '============================================================='
print '(a,i5)' , 'loadcase: ' , loadcase
if ( . not . bc ( loadcase ) % followFormerTrajectory ) print '(a)' , 'drop guessing along trajectory'
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % velGradApplied ) then
do j = 1_pInt , 3_pInt
2011-12-04 15:31:32 +05:30
if ( any ( bc ( loadcase ) % maskDeformation ( j , 1 : 3 ) == . true . ) . and . &
any ( bc ( loadcase ) % maskDeformation ( j , 1 : 3 ) == . false . ) ) errorID = 32_pInt ! each row should be either fully or not at all defined
2011-10-18 20:15:32 +05:30
enddo
2011-12-04 15:31:32 +05:30
print '(a)' , 'velocity gradient:'
2011-10-18 20:15:32 +05:30
else
2011-12-04 15:31:32 +05:30
print '(a)' , 'deformation gradient rate:'
2011-10-18 20:15:32 +05:30
endif
2011-11-15 23:24:18 +05:30
print '(3(3(f12.6,x)/)\)' , merge ( math_transpose3x3 ( bc ( loadcase ) % deformation ) , &
2011-12-04 15:31:32 +05:30
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , ( / 3 , 3 / ) ) , transpose ( bc ( loadcase ) % maskDeformation ) )
print '(a,/,3(3(f12.6,x)/)\)' , 'stress / GPa:' , 1e-9 * merge ( math_transpose3x3 ( bc ( loadcase ) % stress ) , &
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , ( / 3 , 3 / ) ) , &
transpose ( bc ( loadcase ) % maskStress ) )
if ( any ( bc ( loadcase ) % rotation / = math_I3 ) ) &
print '(a,3(3(f12.6,x)/)\)' , 'rotation of loadframe:' , math_transpose3x3 ( bc ( loadcase ) % rotation )
print '(a,f12.6)' , 'temperature:' , bc ( loadcase ) % temperature
print '(a,f12.6)' , 'time: ' , bc ( loadcase ) % timeIncrement
print '(a,i5)' , 'steps: ' , bc ( loadcase ) % steps
print '(a,i5)' , 'output frequency: ' , bc ( loadcase ) % outputfrequency
print '(a,i5)' , 'restart frequency: ' , bc ( loadcase ) % restartfrequency
if ( any ( bc ( loadcase ) % maskStress == bc ( loadcase ) % maskDeformation ) ) errorID = 31 ! exclusive or masking only
if ( any ( bc ( loadcase ) % maskStress . and . transpose ( bc ( loadcase ) % maskStress ) . and . &
reshape ( ( / . false . , . true . , . true . , . true . , . false . , . true . , . true . , . true . , . false . / ) , ( / 3 , 3 / ) ) ) ) &
errorID = 38_pInt ! no rotation is allowed by stress BC
if ( any ( abs ( math_mul33x33 ( bc ( loadcase ) % rotation , math_transpose3x3 ( bc ( loadcase ) % rotation ) ) - math_I3 ) &
> reshape ( spread ( rotation_tol , 1 , 9 ) , ( / 3 , 3 / ) ) ) &
. or . abs ( math_det3x3 ( bc ( loadcase ) % rotation ) ) > 1.0_pReal + rotation_tol ) &
errorID = 46_pInt ! given rotation matrix contains strain
if ( bc ( loadcase ) % timeIncrement < 0.0_pReal ) errorID = 34_pInt ! negative time increment
if ( bc ( loadcase ) % steps < 1_pInt ) errorID = 35_pInt ! non-positive increment count
if ( bc ( loadcase ) % outputfrequency < 1_pInt ) errorID = 36_pInt ! non-positive result frequency
if ( bc ( loadcase ) % restartfrequency < 1_pInt ) errorID = 39_pInt ! non-positive restart frequency
if ( errorID > 0_pInt ) call IO_error ( error_ID = errorID , ext_msg = loadcase_string )
2011-10-18 20:15:32 +05:30
enddo
2011-07-25 22:00:21 +05:30
! Initialization of fftw (see manual on fftw.org for more details)
2011-10-24 23:56:34 +05:30
#ifdef _OPENMP
2011-12-04 15:31:32 +05:30
if ( DAMASK_NumThreadsInt > 0_pInt ) then
2011-10-24 23:56:34 +05:30
call dfftw_init_threads ( ierr )
2011-12-04 15:31:32 +05:30
if ( ierr == 0_pInt ) call IO_error ( error_ID = 104_pInt )
2011-10-24 23:56:34 +05:30
call dfftw_plan_with_nthreads ( DAMASK_NumThreadsInt )
endif
#endif
2011-12-01 17:31:13 +05:30
call dfftw_set_timelimit ( fftw_timelimit )
2011-12-04 15:31:32 +05:30
!*************************************************************
2011-04-06 15:28:17 +05:30
! Loop over loadcases defined in the loadcase file
2011-11-18 03:41:05 +05:30
do loadcase = 1_pInt , N_Loadcases
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
2011-07-13 22:03:12 +05:30
time0 = time ! loadcase start time
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % followFormerTrajectory ) then ! continue to guess along former trajectory where applicable
2011-07-13 22:03:12 +05:30
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step
endif
2011-07-06 18:40:18 +05:30
2011-11-15 23:24:18 +05:30
mask_defgrad = merge ( ones , zeroes , bc ( loadcase ) % maskDeformation )
mask_stress = merge ( ones , zeroes , bc ( loadcase ) % maskStress )
size_reduced = count ( bc ( loadcase ) % maskStressVector )
2011-08-26 19:36:37 +05:30
allocate ( c_reduced ( size_reduced , size_reduced ) ) ; c_reduced = 0.0_pReal
allocate ( s_reduced ( size_reduced , size_reduced ) ) ; s_reduced = 0.0_pReal
2011-08-10 21:32:13 +05:30
2011-11-15 23:24:18 +05:30
timeinc = bc ( loadcase ) % timeIncrement / bc ( loadcase ) % steps ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
fDot = bc ( loadcase ) % deformation ! only valid for given fDot. will be overwritten later in case L is given
2011-11-18 03:41:05 +05:30
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
! loop oper steps defined in input file for current loadcase
2011-11-18 03:41:05 +05:30
do step = 1_pInt , bc ( loadcase ) % steps
2011-11-15 23:24:18 +05:30
!*************************************************************
2011-11-18 03:41:05 +05:30
! forwarding time
if ( bc ( loadcase ) % logscale == 1_pInt ) then ! loglinear scale
if ( loadcase == 1_pInt ) then ! 1st loadcase of loglinear scale
if ( step == 1_pInt ) then ! 1st step of 1st loadcase of loglinear scale
timeinc = bc ( 1 ) % timeIncrement * ( 2.0_pReal ** real ( 1_pInt - bc ( 1 ) % steps , pReal ) ) ! assume 1st step is equal to 2nd
else ! not-1st step of 1st loadcase of loglinear scale
timeinc = bc ( 1 ) % timeIncrement * ( 2.0_pReal ** real ( step - 1_pInt - bc ( 1 ) % steps , pReal ) )
endif
else ! not-1st loadcase of loglinear scale
timeinc = time0 * ( ( 1.0_pReal + bc ( loadcase ) % timeIncrement / time0 ) ** real ( step / bc ( loadcase ) % steps , pReal ) &
- ( 1.0_pReal + bc ( loadcase ) % timeIncrement / time0 ) ** real ( ( step - 1_pInt ) / bc ( loadcase ) % steps , pReal ) )
2011-11-11 19:47:43 +05:30
endif
2011-11-07 16:34:57 +05:30
endif
2011-11-18 03:41:05 +05:30
time = time + timeinc
totalStepsCounter = totalStepsCounter + 1_pInt
2011-11-15 23:24:18 +05:30
!*************************************************************
2011-11-18 03:41:05 +05:30
! Initialization Start
2011-11-15 23:24:18 +05:30
!*************************************************************
2011-11-18 03:41:05 +05:30
if ( totalStepsCounter > = restartReadStep ) then ! Do calculations (otherwise just forwarding)
2011-11-21 23:42:40 +05:30
2011-11-18 03:41:05 +05:30
if ( regrid == . true . ) then ! 'DeInitialize' the values changing in case of regridding
regrid = . false .
call dfftw_destroy_plan ( fftw_plan ( 1 ) ) ; call dfftw_destroy_plan ( fftw_plan ( 2 ) )
if ( debugDivergence ) call dfftw_destroy_plan ( fftw_plan ( 3 ) )
deallocate ( defgrad )
deallocate ( defgradold )
deallocate ( coordinates )
deallocate ( temperature )
deallocate ( xi )
deallocate ( workfft )
!ToDo: here we have to create the new geometry and assign the values from the previous step
endif
2011-11-21 23:42:40 +05:30
if ( totalStepsCounter == restartReadStep ) then ! Initialize values
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step
allocate ( defgrad ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 , 3 ) ) ; defgrad = 0.0_pReal
allocate ( defgradold ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 , 3 ) ) ; defgradold = 0.0_pReal
allocate ( coordinates ( 3 , res ( 1 ) , res ( 2 ) , res ( 3 ) ) ) ; coordinates = 0.0_pReal
allocate ( temperature ( res ( 1 ) , res ( 2 ) , res ( 3 ) ) ) ; temperature = bc ( 1 ) % temperature ! start out isothermally
allocate ( xi ( 3 , res ( 1 ) / 2 + 1 , res ( 2 ) , res ( 3 ) ) ) ; xi = 0.0_pReal
allocate ( workfft ( res ( 1 ) + 2 , res ( 2 ) , res ( 3 ) , 3 , 3 ) ) ; workfft = 0.0_pReal
if ( debugDivergence ) allocate ( divergence ( res ( 1 ) + 2 , res ( 2 ) , res ( 3 ) , 3 ) ) ; divergence = 0.0_pReal
2011-12-04 15:31:32 +05:30
wgt = 1.0_pReal / real ( Npoints , pReal )
2011-11-21 23:42:40 +05:30
call dfftw_plan_many_dft_r2c ( fftw_plan ( 1 ) , 3 , ( / res ( 1 ) , res ( 2 ) , res ( 3 ) / ) , 9 , &
workfft , ( / res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) / ) , 1 , ( res ( 1 ) + 2_pInt ) * res ( 2 ) * res ( 3 ) , &
2011-12-01 17:31:13 +05:30
workfft , ( / res ( 1 ) / 2_pInt + 1_pInt , res ( 2 ) , res ( 3 ) / ) , 1 , ( res ( 1 ) / 2_pInt + 1_pInt ) * res ( 2 ) * res ( 3 ) , fftw_planner_flag )
2011-11-21 23:42:40 +05:30
call dfftw_plan_many_dft_c2r ( fftw_plan ( 2 ) , 3 , ( / res ( 1 ) , res ( 2 ) , res ( 3 ) / ) , 9 , &
workfft , ( / res ( 1 ) / 2_pInt + 1_pInt , res ( 2 ) , res ( 3 ) / ) , 1 , ( res ( 1 ) / 2_pInt + 1_pInt ) * res ( 2 ) * res ( 3 ) , &
2011-12-01 17:31:13 +05:30
workfft , ( / res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) / ) , 1 , ( res ( 1 ) + 2_pInt ) * res ( 2 ) * res ( 3 ) , fftw_planner_flag )
2011-11-21 23:42:40 +05:30
if ( debugDivergence ) &
call dfftw_plan_many_dft_c2r ( fftw_plan ( 3 ) , 3 , ( / res ( 1 ) , res ( 2 ) , res ( 3 ) / ) , 3 , &
divergence , ( / res ( 1 ) / 2_pInt + 1_pInt , res ( 2 ) , res ( 3 ) / ) , 1 , ( res ( 1 ) / 2_pInt + 1_pInt ) * res ( 2 ) * res ( 3 ) , &
2011-12-01 17:31:13 +05:30
divergence , ( / res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) / ) , 1 , ( res ( 1 ) + 2_pInt ) * res ( 2 ) * res ( 3 ) , fftw_planner_flag )
2011-11-18 03:41:05 +05:30
if ( debugGeneral ) then
!$OMP CRITICAL (write2out)
write ( 6 , * ) 'FFTW initialized'
!$OMP END CRITICAL (write2out)
endif
2011-11-21 23:42:40 +05:30
if ( restartReadStep == 1_pInt ) then ! not restarting, no deformation at the beginning
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
enddo ; enddo ; enddo
else ! using old values
if ( IO_read_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , trim ( getSolverJobName ( ) ) , size ( defgrad ) ) ) then
read ( 777 , rec = 1 ) defgrad
close ( 777 )
endif
defgradold = defgrad
defgradAim = 0.0_pReal
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
defgradAim = defgradAim + defgrad ( i , j , k , 1 : 3 , 1 : 3 ) ! calculating old average deformation
enddo ; enddo ; enddo
defgradAim = defgradAim * wgt
defgradAimOld = defgradAim
guessmode = 0.0_pInt
endif
2011-11-21 23:42:40 +05:30
2011-11-18 03:41:05 +05:30
ielem = 0_pInt
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
ielem = ielem + 1_pInt
coordinates ( 1 : 3 , i , j , k ) = mesh_ipCenterOfGravity ( 1 : 3 , 1 , ielem ) ! set to initial coordinates ToDo: SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!! But do we know them? I don't think so. Otherwise we don't need geometry reconstruction
call CPFEM_general ( 2_pInt , coordinates ( 1 : 3 , i , j , k ) , math_I3 , math_I3 , temperature ( i , j , k ) , 0.0_pReal , ielem , 1_pInt , cstress , dsde , pstress , dPdF )
c_current = c_current + dPdF
enddo ; enddo ; enddo
c0_reference = c_current * wgt ! linear reference material stiffness
c_prev = math_rotate_forward3x3x3x3 ( c0_reference , bc ( loadcase ) % rotation ) ! rotate_forward: lab -> load system
if ( debugGeneral ) then
!$OMP CRITICAL (write2out)
2011-12-04 15:31:32 +05:30
write ( 6 , * ) 'first call to CPFEM_general finished'
2011-11-18 03:41:05 +05:30
!$OMP END CRITICAL (write2out)
endif
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
2011-11-18 03:41:05 +05:30
k_s ( 3 ) = k - 1_pInt
2011-11-21 23:42:40 +05:30
if ( k > res ( 3 ) / 2_pInt + 1_pInt ) k_s ( 3 ) = k_s ( 3 ) - res ( 3 )
do j = 1_pInt , res ( 2 )
2011-11-18 03:41:05 +05:30
k_s ( 2 ) = j - 1_pInt
2011-11-21 23:42:40 +05:30
if ( j > res ( 2 ) / 2_pInt + 1_pInt ) k_s ( 2 ) = k_s ( 2 ) - res ( 2 )
do i = 1 , res ( 1 ) / 2_pInt + 1_pInt
2011-11-18 03:41:05 +05:30
k_s ( 1 ) = i - 1_pInt
xi ( 3 , i , j , k ) = 0.0_pReal ! 2D case
2011-11-21 23:42:40 +05:30
if ( res ( 3 ) > 1_pInt ) xi ( 3 , i , j , k ) = real ( k_s ( 3 ) , pReal ) / geomdimension ( 3 ) ! 3D case
2011-11-18 03:41:05 +05:30
xi ( 2 , i , j , k ) = real ( k_s ( 2 ) , pReal ) / geomdimension ( 2 )
xi ( 1 , i , j , k ) = real ( k_s ( 1 ) , pReal ) / geomdimension ( 1 )
enddo ; enddo ; enddo
2011-12-04 15:31:32 +05:30
! remove highest frequencies for calculation of divergence (CAREFUL, they will be used for pre calculatet gamma operator!)
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 ) / 2_pInt + 1_pInt
if ( k == res ( 3 ) / 2_pInt + 1_pInt ) xi ( 3 , i , j , k ) = 0.0_pReal
if ( j == res ( 2 ) / 2_pInt + 1_pInt ) xi ( 2 , i , j , k ) = 0.0_pReal
if ( i == res ( 1 ) / 2_pInt + 1_pInt ) xi ( 1 , i , j , k ) = 0.0_pReal
2011-11-18 03:41:05 +05:30
enddo ; enddo ; enddo
if ( memory_efficient ) then ! allocate just single fourth order tensor
allocate ( gamma_hat ( 1 , 1 , 1 , 3 , 3 , 3 , 3 ) ) ; gamma_hat = 0.0_pReal
else ! precalculation of gamma_hat field
2011-11-21 23:42:40 +05:30
allocate ( gamma_hat ( res ( 1 ) / 2_pInt + 1_pInt , res ( 2 ) , res ( 3 ) , 3 , 3 , 3 , 3 ) ) ; gamma_hat = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 ) / 2_pInt + 1_pInt
2011-11-18 03:41:05 +05:30
if ( any ( xi ( : , i , j , k ) / = 0.0_pReal ) ) then
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
xiDyad ( l , m ) = xi ( l , i , j , k ) * xi ( m , i , j , k )
enddo ; enddo
temp33_Real = math_inv3x3 ( math_mul3333xx33 ( c0_reference , xiDyad ) )
else
xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt ; do p = 1_pInt , 3_pInt
gamma_hat ( i , j , k , l , m , n , p ) = - 0.25 * ( temp33_Real ( l , n ) + temp33_Real ( n , l ) ) * &
( xiDyad ( m , p ) + xiDyad ( p , m ) )
enddo ; enddo ; enddo ; enddo
enddo ; enddo ; enddo
endif
! write header of output file
!$OMP CRITICAL (write2out)
open ( 538 , file = trim ( getSolverWorkingDirectoryName ( ) ) / / trim ( getSolverJobName ( ) ) &
/ / '.spectralOut' , form = 'UNFORMATTED' , status = 'REPLACE' )
write ( 538 ) , 'load' , trim ( getLoadcaseName ( ) )
write ( 538 ) , 'workingdir' , trim ( getSolverWorkingDirectoryName ( ) )
write ( 538 ) , 'geometry' , trim ( getSolverJobName ( ) ) / / InputFileExtension
2011-11-21 23:42:40 +05:30
write ( 538 ) , 'resolution' , res
2011-11-18 03:41:05 +05:30
write ( 538 ) , 'dimension' , geomdimension
write ( 538 ) , 'materialpoint_sizeResults' , materialpoint_sizeResults
write ( 538 ) , 'loadcases' , N_Loadcases
2011-11-21 23:42:40 +05:30
write ( 538 ) , 'logscale' , bc ( 1 : N_Loadcases ) % logscale ! one entry per loadcase (0: linear, 1: log)
write ( 538 ) , 'frequencies' , bc ( 1 : N_Loadcases ) % outputfrequency ! one entry per loadcase
write ( 538 ) , 'times' , bc ( 1 : N_Loadcases ) % timeIncrement ! one entry per loadcase
2011-11-18 03:41:05 +05:30
bc ( 1 ) % steps = bc ( 1 ) % steps + 1_pInt
2011-11-21 23:42:40 +05:30
write ( 538 ) , 'increments' , bc ( 1 : N_Loadcases ) % steps ! one entry per loadcase ToDo: rename keyword to steps
2011-11-18 03:41:05 +05:30
bc ( 1 ) % steps = bc ( 1 ) % steps - 1_pInt
2011-11-21 23:42:40 +05:30
write ( 538 ) , 'startingIncrement' , restartReadStep - 1_pInt ! start with writing out the previous step
write ( 538 ) , 'eoh' ! end of header
2011-12-04 15:31:32 +05:30
write ( 538 ) , materialpoint_results ( 1_pInt : materialpoint_sizeResults , 1 , 1_pInt : Npoints ) ! initial (non-deformed) results
2011-11-18 03:41:05 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
endif
2011-11-18 03:41:05 +05:30
!*************************************************************
! Initialization End
!*************************************************************
if ( mod ( step - 1_pInt , bc ( loadcase ) % restartFrequency ) == 0_pInt ) then ! at frequency of writing restart information
restartWrite = . true . ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
else
restartWrite = . false .
2011-10-25 19:08:24 +05:30
endif
2011-11-18 03:41:05 +05:30
if ( bc ( loadcase ) % velGradApplied ) & ! calculate fDot from given L and current F
fDot = math_mul33x33 ( bc ( loadcase ) % deformation , defgradAim )
!winding forward of deformation aim in loadcase system
temp33_Real = defgradAim
defgradAim = defgradAim &
+ guessmode * mask_stress * ( defgradAim - defgradAimOld ) &
+ mask_defgrad * fDot * timeinc
defgradAimOld = temp33_Real
! update local deformation gradient
if ( any ( bc ( loadcase ) % rotation / = math_I3 ) ) then ! lab and loadcase coordinate system are NOT the same
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
temp33_Real = defgrad ( i , j , k , 1 : 3 , 1 : 3 )
if ( bc ( loadcase ) % velGradApplied ) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33 ( bc ( loadcase ) % deformation , &
math_rotate_forward3x3 ( defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , bc ( loadcase ) % rotation ) )
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = defgrad ( i , j , k , 1 : 3 , 1 : 3 ) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) - defgradold ( i , j , k , 1 : 3 , 1 : 3 ) ) & ! guessing...
+ math_rotate_backward3x3 ( ( 1.0_pReal - guessmode ) * mask_defgrad * fDot , &
bc ( loadcase ) % rotation ) * timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Real
enddo ; enddo ; enddo
else ! one coordinate system for lab and loadcase, save some multiplications
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
temp33_Real = defgrad ( i , j , k , 1 : 3 , 1 : 3 )
if ( bc ( loadcase ) % velGradApplied ) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33 ( bc ( loadcase ) % deformation , defgradold ( i , j , k , 1 : 3 , 1 : 3 ) )
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = defgrad ( i , j , k , 1 : 3 , 1 : 3 ) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) - defgradold ( i , j , k , 1 : 3 , 1 : 3 ) ) & ! guessing...
+ ( 1.0_pReal - guessmode ) * mask_defgrad * fDot * timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Real
enddo ; enddo ; enddo
endif
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
CPFEM_mode = 1_pInt ! winding forward
iter = 0_pInt
err_div = 2.0_pReal * err_div_tol ! go into loop
if ( size_reduced > 0_pInt ) then ! calculate compliance in case stress BC is applied
c_prev99 = math_Plain3333to99 ( c_prev )
k = 0_pInt ! build reduced stiffness
do n = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( n ) ) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( m ) ) then
2011-08-26 19:36:37 +05:30
j = j + 1_pInt
2011-11-18 03:41:05 +05:30
c_reduced ( k , j ) = c_prev99 ( n , m )
endif ; enddo ; endif ; enddo
call math_invert ( size_reduced , c_reduced , s_reduced , i , errmatinv ) ! invert reduced stiffness
if ( errmatinv ) call IO_error ( error_ID = 800 )
s_prev99 = 0.0_pReal ! build full compliance
k = 0_pInt
do n = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( n ) ) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( m ) ) then
j = j + 1_pInt
s_prev99 ( n , m ) = s_reduced ( k , j )
endif ; enddo ; endif ; enddo
s_prev = ( math_Plain99to3333 ( s_prev99 ) )
endif
2011-11-15 23:24:18 +05:30
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-18 03:41:05 +05:30
print '(a)' , '#############################################################'
print '(A,I5.5,A,es12.6)' , 'Increment ' , totalStepsCounter , ' Time ' , time
if ( restartWrite ) then
2011-12-04 15:31:32 +05:30
print '(A)' , 'writing converged results of previous step for restart'
2011-11-18 03:41:05 +05:30
if ( IO_write_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , size ( defgrad ) ) ) then ! and writing deformation gradient field to file
write ( 777 , rec = 1 ) defgrad
close ( 777 )
2011-08-26 19:36:37 +05:30
endif
2011-11-18 03:41:05 +05:30
endif
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-18 03:41:05 +05:30
!*************************************************************
! convergence loop
do while ( iter < itmax . and . &
( err_div > err_div_tol . or . &
err_stress > err_stress_tol ) )
iter = iter + 1_pInt
!*************************************************************
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2011-11-18 03:41:05 +05:30
print '(a)' , '============================================================='
2011-12-04 15:31:32 +05:30
print '(5(a,i5.5))' , 'Loadcase ' , loadcase , ' Step ' , step , '/' , bc ( loadcase ) % steps , '@Iteration ' , iter , '/' , itmax
2011-11-18 03:41:05 +05:30
do n = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
2011-11-21 23:42:40 +05:30
defgrad_av ( m , n ) = sum ( defgrad ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , m , n ) ) * wgt
2011-11-18 03:41:05 +05:30
enddo ; enddo
!$OMP CRITICAL (write2out)
2011-12-04 15:31:32 +05:30
print '(a,/,3(3(f12.7,x)/)\)' , 'deformation gradient:' , math_transpose3x3 ( defgrad_av )
print '(l)' , restartWrite
print '(a)' , '... update stress field P(F) ................................'
2011-11-18 03:41:05 +05:30
!$OMP END CRITICAL (write2out)
2011-11-21 23:42:40 +05:30
defgradDetMax = - huge ( 1.0_pReal )
defgradDetMin = + huge ( 1.0_pReal )
2011-11-18 03:41:05 +05:30
ielem = 0_pInt
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
defgradDet = math_det3x3 ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) )
defgradDetMax = max ( defgradDetMax , defgradDet )
defgradDetMin = min ( defgradDetMin , defgradDet )
ielem = ielem + 1_pInt
call CPFEM_general ( 3_pInt , & ! collect cycle
coordinates ( 1 : 3 , i , j , k ) , defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) , &
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , &
cstress , dsde , pstress , dPdF )
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
2011-11-18 03:41:05 +05:30
2011-12-04 15:31:32 +05:30
print '(a,x,es10.4)' , 'max determinant of deformation:' , defgradDetMax
print '(a,x,es10.4)' , 'min determinant of deformation:' , defgradDetMin
2011-11-18 03:41:05 +05:30
workfft = 0.0_pReal ! needed because of the padding for FFTW
c_current = 0.0_pReal
ielem = 0_pInt
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2011-11-18 03:41:05 +05:30
ielem = ielem + 1_pInt
call CPFEM_general ( CPFEM_mode , & ! first element in first iteration retains CPFEM_mode 1,
coordinates ( 1 : 3 , i , j , k ) , &
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) , & ! others get 2 (saves winding forward effort)
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , &
cstress , dsde , pstress , dPdF )
CPFEM_mode = 2_pInt
workfft ( i , j , k , 1 : 3 , 1 : 3 ) = pstress ! build up average P-K stress
c_current = c_current + dPdF
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
2011-11-18 03:41:05 +05:30
restartWrite = . false . ! ToDo: don't know if we need it. Depends on how CPFEM_general is writing results
do n = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
2011-11-21 23:42:40 +05:30
pstress_av ( m , n ) = sum ( workfft ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , m , n ) ) * wgt
2011-11-18 03:41:05 +05:30
enddo ; enddo
!$OMP CRITICAL (write2out)
2011-12-04 15:31:32 +05:30
print '(a,/,3(3(f12.7,x)/)\)' , 'Piola-Kirchhoff stress / MPa: ' , math_transpose3x3 ( pstress_av ) / 1.e6
2011-11-18 03:41:05 +05:30
err_stress_tol = 0.0_pReal
pstress_av_load = math_rotate_forward3x3 ( pstress_av , bc ( loadcase ) % rotation )
2011-12-04 15:31:32 +05:30
2011-11-18 03:41:05 +05:30
if ( size_reduced > 0_pInt ) then ! calculate stress BC if applied
err_stress = maxval ( abs ( mask_stress * ( pstress_av_load - bc ( loadcase ) % stress ) ) ) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval ( abs ( mask_defgrad * pstress_av_load ) ) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
2011-12-04 15:31:32 +05:30
print '(a)' , ''
print '(a,es10.4,a,f6.2)' , 'error stress = ' , err_stress , ', ' , err_stress / err_stress_tol
print '(a)' , '... correcting deformation gradient to fulfill BCs ..........'
2011-11-18 03:41:05 +05:30
defgradAimCorr = - math_mul3333xx33 ( s_prev , ( ( pstress_av_load - bc ( loadcase ) % stress ) ) ) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr
2011-12-04 15:31:32 +05:30
print '(a,/,3(3(f12.7,x)/)\)' , 'new deformation aim: ' , &
math_transpose3x3 ( math_rotate_backward3x3 ( defgradAim , bc ( loadcase ) % rotation ) )
print '(a,x,es10.4)' , 'with determinant: ' , math_det3x3 ( defgradAim )
2011-11-18 03:41:05 +05:30
endif
2011-12-04 15:31:32 +05:30
print '(a)' , ''
print '(a)' , '... calculating equilibrium with spectral method ............'
2011-11-18 03:41:05 +05:30
!$OMP END CRITICAL (write2out)
call dfftw_execute_dft_r2c ( fftw_plan ( 1 ) , workfft , workfft ) ! FFT of pstress
p_hat_avg = sqrt ( maxval ( math_eigenvalues3x3 ( math_mul33x33 ( workfft ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) , & ! L_2 norm of average stress in fourier space,
math_transpose3x3 ( workfft ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) ) ) ) ) ) ! ignore imaginary part as it is always zero for real only input))
err_div = 0.0_pReal
err_div_max = 0.0_pReal
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 ) / 2_pInt + 1_pInt
2011-11-18 03:41:05 +05:30
err_div = err_div + sqrt ( sum ( ( & ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain)
math_mul33x3_complex ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) + &
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img , &
xi ( 1 : 3 , i , j , k ) ) &
) ** 2.0_pReal ) )
if ( debugDivergence ) &
2011-11-21 23:42:40 +05:30
err_div_max = max ( err_div_max , abs ( sqrt ( sum ( ( & ! maximum of L two norm of div(stress) in fourier space (Suquet large strain)
2011-11-18 03:41:05 +05:30
math_mul33x3_complex ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) + &
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img , &
xi ( 1 : 3 , i , j , k ) ) &
) ** 2.0_pReal ) ) ) )
2011-11-15 23:24:18 +05:30
enddo ; enddo ; enddo
2011-12-04 15:31:32 +05:30
if ( divergence_correction ) then
if ( res ( 3 ) == 1_pInt ) then
correctionFactor = minval ( geomdimension ( 1 : 2 ) ) * wgt ** ( - 1.0_pReal / 4.0_pReal ) ! 2D case, ToDo: correct? PE: Do we need this in the loop or can be pre-calculated?
else
correctionFactor = minval ( geomdimension ( 1 : 3 ) ) * wgt ** ( - 1.0_pReal / 4.0_pReal ) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency
endif
else
correctionFactor = 1.0_pReal
endif
2011-11-18 03:41:05 +05:30
err_div = err_div * wgt / p_hat_avg * correctionFactor ! weighting by points and average stress and multiplying with correction factor
err_div_max = err_div_max / p_hat_avg * correctionFactor ! weighting by average stress and multiplying with correction factor
if ( memory_efficient ) then ! memory saving version, on-the-fly calculation of gamma_hat
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 ) / 2_pInt + 1_pInt
2011-11-18 03:41:05 +05:30
if ( any ( xi ( : , i , j , k ) / = 0.0_pReal ) ) then
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
xiDyad ( l , m ) = xi ( l , i , j , k ) * xi ( m , i , j , k )
enddo ; enddo
temp33_Real = math_inv3x3 ( math_mul3333xx33 ( c0_reference , xiDyad ) )
else
xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt ; do p = 1_pInt , 3_pInt
gamma_hat ( 1 , 1 , 1 , l , m , n , p ) = - 0.25_pReal * ( temp33_Real ( l , n ) + temp33_Real ( n , l ) ) * &
( xiDyad ( m , p ) + xiDyad ( p , m ) )
enddo ; enddo ; enddo ; enddo
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
temp33_Complex ( m , n ) = sum ( gamma_hat ( 1 , 1 , 1 , m , n , 1 : 3 , 1 : 3 ) * ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) &
+ workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img ) )
enddo ; enddo
workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) = real ( temp33_Complex )
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) = aimag ( temp33_Complex )
enddo ; enddo ; enddo
else ! use precalculated gamma-operator
2011-11-21 23:42:40 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 ) / 2_pInt + 1_pInt
2011-11-18 03:41:05 +05:30
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
temp33_Complex ( m , n ) = sum ( gamma_hat ( i , j , k , m , n , 1 : 3 , 1 : 3 ) * ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) &
+ workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img ) )
enddo ; enddo
workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) = real ( temp33_Complex )
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) = aimag ( temp33_Complex )
enddo ; enddo ; enddo
endif
if ( debugDivergence ) then
2011-11-21 23:42:40 +05:30
divergence = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 ) / 2_pInt + 1_pInt
! real part at i*2-1, imaginary part at i*2 and multiply by i ==> switch and change sign
divergence ( i * 2_pInt - 1_pInt , j , k , 1 : 3 ) = workfft ( i * 2_pInt , j , k , 1 : 3 , 1 ) * xi ( 1 : 3 , i , j , k ) * pi * 2.0_pReal &
+ workfft ( i * 2_pInt , j , k , 1 : 3 , 2 ) * xi ( 1 : 3 , i , j , k ) * pi * 2.0_pReal &
+ workfft ( i * 2_pInt , j , k , 1 : 3 , 3 ) * xi ( 1 : 3 , i , j , k ) * pi * 2.0_pReal
divergence ( i * 2_pInt , j , k , 1 : 3 ) = - workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 ) * xi ( 1 : 3 , i , j , k ) * pi * 2.0_pReal &
- workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 2 ) * xi ( 1 : 3 , i , j , k ) * pi * 2.0_pReal &
- workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 3 ) * xi ( 1 : 3 , i , j , k ) * pi * 2.0_pReal
2011-11-18 03:41:05 +05:30
enddo ; enddo ; enddo
2011-11-21 23:42:40 +05:30
divergence = divergence * correctionFactor
2011-11-18 03:41:05 +05:30
call dfftw_execute_dft_c2r ( fftw_plan ( 3 ) , divergence , divergence )
2011-11-21 23:42:40 +05:30
divergence = divergence * wgt
err_real_div_avg = 0.0_pReal
err_real_div_max = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
err_real_div_avg = err_real_div_avg + sqrt ( sum ( ( divergence ( i , j , k , 1 : 3 ) ) ** 2.0_pReal ) ) ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain)
err_real_div_max = max ( err_real_div_max , abs ( sqrt ( sum ( ( divergence ( i , j , k , 1 : 3 ) ) ** 2.0_pReal ) ) ) ) ! maximum of L two norm of div(stress) in fourier space (Suquet large strain)
enddo ; enddo ; enddo
p_real_avg = sqrt ( maxval ( math_eigenvalues3x3 ( math_mul33x33 ( pstress_av , & ! L_2 norm of average stress in fourier space,
math_transpose3x3 ( pstress_av ) ) ) ) ) ! ignore imaginary part as it is always zero for real only input))
err_real_div_avg = err_real_div_avg * wgt / p_real_avg
err_real_div_max = err_real_div_max / p_real_avg
2011-11-18 03:41:05 +05:30
endif
! average strain
workfft ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) = defgrad_av - math_I3 ! zero frequency (real part)
workfft ( 2 , 1 , 1 , 1 : 3 , 1 : 3 ) = 0.0_pReal ! zero frequency (imaginary part)
call dfftw_execute_dft_c2r ( fftw_plan ( 2 ) , workfft , workfft )
2011-11-21 23:42:40 +05:30
defgrad = defgrad + workfft ( 1 : res ( 1 ) , : , : , : , : ) * wgt
2011-12-04 15:31:32 +05:30
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
2011-11-18 03:41:05 +05:30
defgrad_av ( m , n ) = sum ( defgrad ( : , : , : , m , n ) ) * wgt
enddo ; enddo
defgradAim_lab = math_rotate_backward3x3 ( defgradAim , bc ( loadcase ) % rotation )
2011-12-04 15:31:32 +05:30
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
2011-11-18 03:41:05 +05:30
defgrad ( : , : , : , m , n ) = defgrad ( : , : , : , m , n ) + ( defgradAim_lab ( m , n ) - defgrad_av ( m , n ) ) ! anticipated target minus current state
enddo ; enddo
!$OMP CRITICAL (write2out)
2011-12-04 15:31:32 +05:30
if ( debugDivergence ) then
print '(a,es10.4,a,f6.2)' , 'error divergence FT avg = ' , err_div , ', ' , err_div / err_div_tol
print '(a,es10.4)' , 'error divergence FT max = ' , err_div_max
print '(a,es10.4)' , 'error divergence Real avg = ' , err_real_div_avg
print '(a,es10.4)' , 'error divergence Real max = ' , err_real_div_max
2011-11-21 23:42:40 +05:30
else
2011-12-04 15:31:32 +05:30
print '(a,es10.4,a,f6.2)' , 'error divergence = ' , err_div , ', ' , err_div / err_div_tol
2011-11-21 23:42:40 +05:30
endif
2011-11-18 03:41:05 +05:30
!$OMP END CRITICAL (write2out)
enddo ! end looping when convergency is achieved
2011-08-26 19:36:37 +05:30
2011-11-18 03:41:05 +05:30
c_prev = math_rotate_forward3x3x3x3 ( c_current * wgt , bc ( loadcase ) % rotation ) ! calculate stiffness for next step
!ToDo: Incfluence for next loadcase
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2011-11-18 03:41:05 +05:30
print '(a)' , '============================================================='
2011-12-04 15:31:32 +05:30
if ( err_div > err_div_tol . or . err_stress > err_stress_tol ) then
print '(A,I5.5,A)' , 'increment ' , totalStepsCounter , ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
2011-11-18 03:41:05 +05:30
else
2011-12-04 15:31:32 +05:30
print '(A,I5.5,A)' , 'increment ' , totalStepsCounter , ' converged'
2011-11-18 03:41:05 +05:30
endif
2011-12-04 15:31:32 +05:30
if ( mod ( totalStepsCounter - 1_pInt , bc ( loadcase ) % outputfrequency ) == 0_pInt ) then ! at output frequency
print '(a)' , ''
print '(a)' , '... writing results to file .................................'
write ( 538 ) , materialpoint_results ( 1_pInt : materialpoint_sizeResults , 1 , 1_pInt : Npoints ) ! write result to file
2011-11-18 03:41:05 +05:30
endif
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
endif
2010-09-06 15:30:59 +05:30
enddo ! end looping over steps in current loadcase
2011-08-26 19:36:37 +05:30
deallocate ( c_reduced )
deallocate ( s_reduced )
enddo ! end looping over loadcases
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2011-11-15 23:24:18 +05:30
print '(a)' , '#############################################################'
2011-12-04 15:31:32 +05:30
print '(a,i5.5,a,i5.5,a)' , 'of ' , totalStepsCounter - restartReadStep + 1_pInt , ' calculated steps, ' , notConvergedCounter , ' steps did not converge!'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-01-07 18:26:45 +05:30
close ( 538 )
2011-10-18 20:15:32 +05:30
call dfftw_destroy_plan ( fftw_plan ( 1 ) ) ; call dfftw_destroy_plan ( fftw_plan ( 2 ) )
2011-12-04 15:31:32 +05:30
if ( debugDivergence ) call dfftw_destroy_plan ( fftw_plan ( 3 ) )
2011-05-11 22:08:45 +05:30
end program DAMASK_spectral
2010-06-10 20:21:10 +05:30
!********************************************************************
! quit subroutine to satisfy IO_error
!
!********************************************************************
2010-06-08 15:38:15 +05:30
subroutine quit ( id )
use prec
implicit none
integer ( pInt ) id
stop
2011-05-24 21:27:59 +05:30
end subroutine