2012-01-25 19:57:26 +05:30
! Copyright 2012 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/>.
!
2012-01-25 19:57:26 +05:30
!##################################################################################################
2010-06-08 15:40:57 +05:30
!* $Id$
2012-01-25 19:57:26 +05:30
!##################################################################################################
2010-06-08 15:38:15 +05:30
! Material subroutine for BVP solution using spectral method
!
2012-01-04 23:13:26 +05:30
! Run 'DAMASK_spectral.exe --help' to get usage hints
!
2010-06-08 15:38:15 +05:30
! 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
2012-01-25 19:57:26 +05:30
!##################################################################################################
! used modules
!##################################################################################################
2011-05-11 22:08:45 +05:30
program DAMASK_spectral
2010-06-08 15:38:15 +05:30
2011-05-11 22:31:03 +05:30
use DAMASK_interface
2012-01-25 19:57:26 +05:30
use prec , only : pInt , pReal , DAMASK_NaN
use IO
use debug , only : debug_spectral , &
debug_spectralGeneral , &
debug_spectralDivergence , &
debug_spectralRestart , &
debug_spectralFFTW
2010-07-05 21:31:36 +05:30
use math
2012-01-04 23:13:26 +05:30
use kdtree2_module
2012-01-25 19:57:26 +05:30
use CPFEM , only : CPFEM_general , CPFEM_initAll
2012-02-10 17:29:59 +05:30
use FEsolving , only : restartWrite , restartInc
2012-01-25 19:57:26 +05:30
use numerics , only : err_div_tol , err_stress_tolrel , rotation_tol , itmax , &
memory_efficient , update_gamma , &
simplified_algorithm , divergence_correction , &
DAMASK_NumThreadsInt , &
fftw_planner_flag , fftw_timelimit
use homogenization , only : materialpoint_sizeResults , materialpoint_results
!$ use OMP_LIB ! the openMP function library
!##################################################################################################
! variable declaration
!##################################################################################################
2010-06-08 15:38:15 +05:30
implicit none
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! variables to read from load case and geom file
real ( pReal ) , dimension ( 9 ) :: temp_valueVector ! stores information temporarily from loadcase file
logical , dimension ( 9 ) :: temp_maskVector
integer ( pInt ) , parameter :: maxNchunksLoadcase = ( 1_pInt + 9_pInt ) * 3_pInt + & ! deformation, rotation, and stress
( 1_pInt + 1_pInt ) * 5_pInt + & ! time, (log)incs, temp, restartfrequency, and outputfrequency
1_pInt , & ! dropguessing
maxNchunksGeom = 7_pInt , & ! 4 identifiers, 3 values
myUnit = 234_pInt
integer ( pInt ) , dimension ( 1_pInt + maxNchunksLoadcase * 2_pInt ) :: positions ! this is longer than needed for geometry parsing
integer ( pInt ) :: headerLength , &
N_l = 0_pInt , &
N_t = 0_pInt , &
N_n = 0_pInt , &
N_Fdot = 0_pInt
2011-10-18 20:15:32 +05:30
character ( len = 1024 ) :: path , line , keyword
2012-01-25 19:57:26 +05:30
logical :: gotResolution = . false . , &
gotDimension = . false . , &
gotHomogenization = . false .
2011-11-15 23:24:18 +05:30
type bc_type
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: deformation = 0.0_pReal , & ! applied velocity gradient or time derivative of deformation gradient
stress = 0.0_pReal , & ! stress BC (if applicable)
rotation = math_I3 ! rotation of BC (if applicable)
real ( pReal ) :: time = 0.0_pReal , & ! length of increment
2012-02-13 22:45:02 +05:30
temperature = 30 0.0_pReal ! isothermal starting conditions
2012-01-25 19:57:26 +05:30
integer ( pInt ) :: incs = 0_pInt , & ! number of increments
outputfrequency = 1_pInt , & ! frequency of result writes
restartfrequency = 0_pInt , & ! frequency of restart writes
logscale = 0_pInt ! linear/logaritmic time inc flag
logical :: followFormerTrajectory = . true . , & ! follow trajectory of former loadcase
velGradApplied = . false . ! decide wether velocity gradient or fdot is given
logical , dimension ( 3 , 3 ) :: maskDeformation = . false . , & ! mask of deformation boundary conditions
maskStress = . false . ! mask of stress boundary conditions
logical , dimension ( 9 ) :: maskStressVector = . false . ! linear mask of boundary conditions
2011-11-15 23:24:18 +05:30
end type
2012-01-13 21:48:16 +05:30
2011-12-04 15:31:32 +05:30
type ( bc_type ) , allocatable , dimension ( : ) :: bc
2012-01-12 15:53:05 +05:30
character ( len = 6 ) :: loadcase_string
2010-07-01 20:50:06 +05:30
2012-01-25 19:57:26 +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
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( 3 ) :: geomdim = 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
integer ( pInt ) :: res1_red ! to store res(1)/2 +1
2010-07-01 20:50:06 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2011-11-07 23:55:10 +05:30
! stress, stiffness and compliance average etc.
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: pstress , pstress_av , &
defgradAim = math_I3 , defgradAimOld = math_I3 , &
2012-02-10 17:29:59 +05:30
mask_stress , mask_defgrad , deltaF , &
2012-01-25 19:57:26 +05:30
pstress_av_lab , defgradAim_lab , defgrad_av_lab ! quantities rotated to other coordinate system
2012-02-10 17:29:59 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dPdF , c0_reference , c_current = 0.0_pReal , s_prev , c_prev , & ! stiffness and compliance
s0_reference
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( 6 ) :: cstress ! cauchy stress
2012-02-10 17:29:59 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: dsde , c0_66 , s0_66 ! small strain stiffness
2012-02-16 00:28:38 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: s_prev99 , c_prev99 ! compliance and stiffness in matrix notation
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( : , : ) , allocatable :: s_reduced , c_reduced ! reduced compliance and stiffness (only for stress BC)
2012-02-16 00:28:38 +05:30
integer ( pInt ) :: size_reduced = 0_pInt ! number of stress BCs
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2011-10-18 20:15:32 +05:30
! pointwise data
2012-01-25 19:57:26 +05:30
type ( C_PTR ) :: tensorField , tau ! fields in real an fourier space
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: tensorField_real ! fields in real space (pointer)
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: tau_real
2012-02-10 17:29:59 +05:30
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: tensorField_fourier ! fields in fourier space (pointer)
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: tau_fourier
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: defgrad , defgradold
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: coordinates
real ( pReal ) , dimension ( : , : , : ) , allocatable :: temperature
!--------------------------------------------------------------------------------------------------
2011-10-25 19:08:24 +05:30
! variables storing information for spectral method and FFTW
2012-01-25 19:57:26 +05:30
type ( C_PTR ) :: plan_stress , plan_correction , plan_tau ! plans for fftw
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 for divergence and for gamma operator
2012-02-16 00:28:38 +05:30
integer ( pInt ) , dimension ( 3 ) :: k_s
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2011-02-07 20:05:42 +05:30
! loop variables, convergence etc.
2012-02-10 17:29:59 +05:30
real ( pReal ) :: time = 0.0_pReal , time0 = 0.0_pReal , timeinc = 1.0_pReal , timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval
2012-02-01 00:48:55 +05:30
real ( pReal ) :: guessmode , err_div , err_stress , err_stress_tol
2011-11-15 23:24:18 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: ones = 1.0_pReal , zeroes = 0.0_pReal
2012-01-13 21:48:16 +05:30
complex ( pReal ) , dimension ( 3 ) :: temp3_Complex
2011-11-15 23:24:18 +05:30
complex ( pReal ) , dimension ( 3 , 3 ) :: temp33_Complex
2012-01-13 21:48:16 +05:30
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
2012-01-13 21:48:16 +05:30
integer ( pInt ) :: N_Loadcases , loadcase , inc , iter , ielem , CPFEM_mode , &
2012-02-10 17:29:59 +05:30
ierr , totalIncsCounter = 0_pInt , &
notConvergedCounter = 0_pInt , convergedCounter = 0_pInt
2012-01-13 21:48:16 +05:30
logical :: errmatinv
2012-01-25 19:57:26 +05:30
real ( pReal ) :: defgradDet , correctionFactor
!--------------------------------------------------------------------------------------------------
!variables controlling debugging
logical :: debugGeneral , debugDivergence , debugRestart , debugFFTW
!--------------------------------------------------------------------------------------------------
!variables for additional output due to general debugging
2012-02-16 00:28:38 +05:30
real ( pReal ) :: defgradDetMax , defgradDetMin , maxCorrectionSym , maxCorrectionSkew
2011-11-15 23:24:18 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! variables for additional output of divergence calculations
type ( C_PTR ) :: divergence , plan_divergence
2012-01-13 21:48:16 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , pointer :: divergence_real
2012-02-10 17:29:59 +05:30
complex ( pReal ) , dimension ( : , : , : , : ) , pointer :: divergence_fourier
2012-01-25 19:57:26 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: divergence_postProc
2012-02-10 17:29:59 +05:30
real ( pReal ) :: pstress_av_L2 , err_div_RMS , err_real_div_RMS , &
2012-01-25 19:57:26 +05:30
err_div_max , err_real_div_max , &
max_div_error
!--------------------------------------------------------------------------------------------------
! variables for debugging fft using a scalar field
2012-02-10 17:29:59 +05:30
type ( C_PTR ) :: scalarField_realC , scalarField_fourierC , &
2012-01-30 19:22:41 +05:30
plan_scalarField_forth , plan_scalarField_back
2012-02-10 17:29:59 +05:30
complex ( pReal ) , dimension ( : , : , : ) , pointer :: scalarField_real
complex ( pReal ) , dimension ( : , : , : ) , pointer :: scalarField_fourier
2012-01-13 21:48:16 +05:30
integer ( pInt ) :: row , column
2012-01-25 19:57:26 +05:30
!##################################################################################################
! reading of information from load case file and geometry file
!##################################################################################################
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
2011-11-04 01:02:11 +05:30
call DAMASK_interface_init ( )
2011-12-23 18:00:35 +05:30
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2012-01-13 21:48:16 +05:30
print '(a)' , ' <<<+- DAMASK_spectral init -+>>>'
print '(a)' , ' $Id$'
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
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)' , ''
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! reading the load case file and allocate data structure containing load cases
2011-11-04 01:02:11 +05:30
path = getLoadcaseName ( )
2012-02-13 23:11:27 +05:30
call IO_open_file ( myUnit , 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
2012-01-25 19:57:26 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2012-01-12 15:53:05 +05:30
positions = IO_stringPos ( line , maxNchunksLoadcase )
2012-01-25 19:57:26 +05:30
do i = 1_pInt , maxNchunksLoadcase , 1_pInt ! reading compulsory parameters for loadcase
2012-01-12 15:53:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , positions , i ) ) )
2012-01-13 21:48:16 +05:30
case ( 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' )
2011-11-15 23:24:18 +05:30
N_l = N_l + 1_pInt
2012-02-10 23:15:45 +05:30
case ( 'fdot' , 'dotf' )
2011-11-15 23:24:18 +05:30
N_Fdot = N_Fdot + 1_pInt
2012-01-13 21:48:16 +05:30
case ( 't' , 'time' , 'delta' )
2011-11-15 23:24:18 +05:30
N_t = N_t + 1_pInt
2012-01-13 21:48:16 +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
2012-01-25 19:57:26 +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
2012-01-25 19:57:26 +05:30
if ( ( N_l + N_Fdot / = N_n ) . or . ( N_n / = N_t ) ) & ! sanity check
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 837_pInt , ext_msg = trim ( path ) ) ! error message for incomplete loadcase
2011-11-15 23:24:18 +05:30
allocate ( bc ( N_Loadcases ) )
2010-06-10 20:21:10 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! reading the load case 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
2012-01-25 19:57:26 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2011-11-15 23:24:18 +05:30
loadcase = loadcase + 1_pInt
2012-01-12 15:53:05 +05:30
positions = IO_stringPos ( line , maxNchunksLoadcase )
2011-11-15 23:24:18 +05:30
do j = 1_pInt , maxNchunksLoadcase
2012-01-12 15:53:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , positions , j ) ) )
2012-02-13 18:11:44 +05:30
case ( 'fdot' , 'dotf' , 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' ) ! assign values for the deformation BC matrix
2012-01-25 19:57:26 +05:30
bc ( loadcase ) % velGradApplied = &
( IO_lc ( IO_stringValue ( line , positions , j ) ) == 'l' . or . & ! in case of given L, set flag to true
IO_lc ( IO_stringValue ( line , positions , j ) ) == 'velocitygrad' . or . &
IO_lc ( IO_stringValue ( line , positions , j ) ) == 'velgrad' . or . &
IO_lc ( IO_stringValue ( line , positions , j ) ) == 'velocitygradient' )
2011-11-15 23:24:18 +05:30
temp_valueVector = 0.0_pReal
temp_maskVector = . false .
2012-01-12 15:53:05 +05:30
forall ( k = 1_pInt : 9_pInt ) temp_maskVector ( k ) = IO_stringValue ( line , positions , j + k ) / = '*'
2011-11-15 23:24:18 +05:30
do k = 1_pInt , 9_pInt
2012-01-12 15:53:05 +05:30
if ( temp_maskVector ( k ) ) temp_valueVector ( k ) = IO_floatValue ( line , positions , j + k )
2010-06-10 20:21:10 +05:30
enddo
2012-02-10 17:29:59 +05:30
bc ( loadcase ) % maskDeformation = transpose ( reshape ( temp_maskVector , [ 3 , 3 ] ) )
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % deformation = math_plain9to33 ( temp_valueVector )
2012-01-13 21:48:16 +05:30
case ( 'p' , 'pk1' , 'piolakirchhoff' , 'stress' )
2011-11-15 23:24:18 +05:30
temp_valueVector = 0.0_pReal
2012-01-25 19:57:26 +05:30
forall ( k = 1_pInt : 9_pInt ) bc ( loadcase ) % maskStressVector ( k ) = &
IO_stringValue ( line , positions , j + k ) / = '*'
2011-11-15 23:24:18 +05:30
do k = 1_pInt , 9_pInt
2012-01-25 19:57:26 +05:30
if ( bc ( loadcase ) % maskStressVector ( k ) ) temp_valueVector ( k ) = &
IO_floatValue ( line , positions , j + k ) ! assign values for the bc(loadcase)%stress matrix
2010-06-10 20:21:10 +05:30
enddo
2012-02-10 17:29:59 +05:30
bc ( loadcase ) % maskStress = transpose ( reshape ( bc ( loadcase ) % maskStressVector , [ 3 , 3 ] ) )
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % stress = math_plain9to33 ( temp_valueVector )
2012-01-25 19:57:26 +05:30
case ( 't' , 'time' , 'delta' ) ! increment time
2012-01-13 21:48:16 +05:30
bc ( loadcase ) % time = IO_floatValue ( line , positions , j + 1_pInt )
2012-01-25 19:57:26 +05:30
case ( 'temp' , 'temperature' ) ! starting temperature
2012-01-12 15:53:05 +05:30
bc ( loadcase ) % temperature = IO_floatValue ( line , positions , j + 1_pInt )
2012-01-25 19:57:26 +05:30
case ( 'n' , 'incs' , 'increments' , 'steps' ) ! number of increments
2012-01-13 21:48:16 +05:30
bc ( loadcase ) % incs = IO_intValue ( line , positions , j + 1_pInt )
2012-01-25 19:57:26 +05:30
case ( 'logincs' , 'logincrements' , 'logsteps' ) ! number of increments (switch to log time scaling)
2012-01-13 21:48:16 +05:30
bc ( loadcase ) % incs = IO_intValue ( line , positions , j + 1_pInt )
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % logscale = 1_pInt
2012-01-25 19:57:26 +05:30
case ( 'f' , 'freq' , 'frequency' , 'outputfreq' ) ! frequency of result writings
2012-01-12 15:53:05 +05:30
bc ( loadcase ) % outputfrequency = IO_intValue ( line , positions , j + 1_pInt )
2012-01-25 19:57:26 +05:30
case ( 'r' , 'restart' , 'restartwrite' ) ! frequency of writing restart information
2012-01-12 15:53:05 +05:30
bc ( loadcase ) % restartfrequency = max ( 0_pInt , IO_intValue ( line , positions , j + 1_pInt ) )
2011-07-14 15:07:31 +05:30
case ( 'guessreset' , 'dropguessing' )
2012-01-25 19:57:26 +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
2012-01-12 15:53:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , positions , j + 1_pInt ) ) )
2011-10-18 20:15:32 +05:30
case ( 'deg' , 'degree' )
2012-01-25 19:57:26 +05:30
p = 1_pInt ! for conversion from degree to radian
2011-10-18 20:15:32 +05:30
case ( 'rad' , 'radian' )
case default
2012-01-25 19:57:26 +05:30
l = 0_pInt ! immediately reading in angles, assuming radians
2011-10-18 20:15:32 +05:30
end select
2012-01-25 19:57:26 +05:30
forall ( k = 1_pInt : 3_pInt ) temp33_Real ( k , 1 ) = &
IO_floatValue ( line , positions , j + l + k ) * real ( p , pReal ) * inRad
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % rotation = math_EulerToR ( temp33_Real ( : , 1 ) )
2012-01-25 19:57:26 +05:30
case ( 'rotation' , 'rot' ) ! assign values for the rotation of loadcase matrix
2011-11-15 23:24:18 +05:30
temp_valueVector = 0.0_pReal
2012-01-12 15:53:05 +05:30
forall ( k = 1_pInt : 9_pInt ) temp_valueVector ( k ) = IO_floatValue ( line , positions , j + k )
2011-11-15 23:24:18 +05:30
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
2012-01-25 19:57:26 +05:30
!-------------------------------------------------------------------------------------------------- ToDo: if temperature at CPFEM is treated properly, move this up immediately after interface init
! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry)
call CPFEM_initAll ( bc ( 1 ) % temperature , 1_pInt , 1_pInt )
2012-02-13 23:11:27 +05:30
if ( update_gamma . and . . not . memory_efficient ) call IO_error ( error_ID = 847_pInt )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! read header of geom file to get size information. 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
2012-02-13 23:11:27 +05:30
call IO_open_file ( myUnit , trim ( path ) / / InputFileExtension )
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
read ( myUnit , '(a1024)' ) line
2012-01-12 15:53:05 +05:30
positions = IO_stringPos ( line , 2_pInt )
keyword = IO_lc ( IO_StringValue ( line , positions , 2_pInt ) )
2011-10-18 20:15:32 +05:30
if ( keyword ( 1 : 4 ) == 'head' ) then
2012-01-12 15:53:05 +05:30
headerLength = IO_intValue ( line , positions , 1_pInt ) + 1_pInt
2011-10-18 20:15:32 +05:30
else
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 842_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
2012-01-12 15:53:05 +05:30
positions = IO_stringPos ( line , maxNchunksGeom )
select case ( IO_lc ( IO_StringValue ( line , positions , 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
2012-01-12 15:53:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , positions , j ) ) )
2010-09-22 17:34:43 +05:30
case ( 'x' )
2012-01-25 19:57:26 +05:30
geomdim ( 1 ) = IO_floatValue ( line , positions , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'y' )
2012-01-25 19:57:26 +05:30
geomdim ( 2 ) = IO_floatValue ( line , positions , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'z' )
2012-01-25 19:57:26 +05:30
geomdim ( 3 ) = IO_floatValue ( line , positions , 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 .
2012-01-12 15:53:05 +05:30
homog = IO_intValue ( line , positions , 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
2012-01-12 15:53:05 +05:30
select case ( IO_lc ( IO_stringValue ( line , positions , j ) ) )
2010-09-22 17:34:43 +05:30
case ( 'a' )
2012-01-12 15:53:05 +05:30
res ( 1 ) = IO_intValue ( line , positions , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'b' )
2012-01-12 15:53:05 +05:30
res ( 2 ) = IO_intValue ( line , positions , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'c' )
2012-01-12 15:53:05 +05:30
res ( 3 ) = IO_intValue ( line , positions , j + 1_pInt )
2010-09-22 17:34:43 +05:30
end select
enddo
2010-06-25 17:01:05 +05:30
end select
enddo
2011-10-18 20:15:32 +05:30
close ( myUnit )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks of geometry parameters
if ( . not . ( gotDimension . and . gotHomogenization . and . gotResolution ) ) &
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 845_pInt )
if ( any ( geomdim < = 0.0_pReal ) ) call IO_error ( error_ID = 802_pInt )
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 . &
2012-01-25 19:57:26 +05:30
( mod ( res ( 3 ) , 2_pInt ) / = 0_pInt . and . res ( 3 ) / = 1_pInt ) ) &
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 803_pInt )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! variables derived from resolution
res1_red = res ( 1 ) / 2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
2011-12-04 15:31:32 +05:30
Npoints = res ( 1 ) * res ( 2 ) * res ( 3 )
2012-01-13 21:48:16 +05:30
wgt = 1.0_pReal / real ( Npoints , pReal )
2011-12-04 15:31:32 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! output of geometry
2011-12-04 15:31:32 +05:30
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)' , '============================================================='
2012-01-25 19:57:26 +05:30
print '(a,3(i12 ))' , 'resolution a b c:' , res
print '(a,3(f12.5))' , 'dimension x y z:' , geomdim
2011-12-04 15:31:32 +05:30
print '(a,i5)' , 'homogenization: ' , homog
print '(a)' , '#############################################################'
print '(a,a)' , 'loadcase file: ' , trim ( getLoadcaseName ( ) )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! consistency checks and output of load case
bc ( 1 ) % followFormerTrajectory = . false . ! cannot guess along trajectory for first inc of first loadcase
2011-12-04 15:31:32 +05:30
errorID = 0_pInt
2011-11-11 19:47:43 +05:30
do loadcase = 1_pInt , N_Loadcases
2012-01-12 15:53:05 +05:30
write ( loadcase_string , '(i6)' ) loadcase
2011-12-04 15:31:32 +05:30
print '(a)' , '============================================================='
2012-01-12 15:53:05 +05:30
print '(a,i6)' , 'loadcase: ' , loadcase
2011-12-04 15:31:32 +05:30
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-06 22:28:17 +05:30
if ( any ( bc ( loadcase ) % maskDeformation ( j , 1 : 3 ) . eqv . . true . ) . and . &
2012-02-13 23:11:27 +05:30
any ( bc ( loadcase ) % maskDeformation ( j , 1 : 3 ) . eqv . . false . ) ) errorID = 832_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
2012-02-01 00:48:55 +05:30
write ( * , '(3(3(f12.7,1x)/))' , advance = 'no' ) merge ( math_transpose33 ( bc ( loadcase ) % deformation ) , &
2012-02-10 17:29:59 +05:30
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , [ 3 , 3 ] ) , transpose ( bc ( loadcase ) % maskDeformation ) )
2012-02-01 00:48:55 +05:30
write ( * , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) ' stress / GPa:' , &
2012-02-16 00:28:38 +05:30
1e-9_pReal * merge ( math_transpose33 ( bc ( loadcase ) % stress ) , &
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , [ 3 , 3 ] ) , transpose ( bc ( loadcase ) % maskStress ) )
2011-12-04 15:31:32 +05:30
if ( any ( bc ( loadcase ) % rotation / = math_I3 ) ) &
2012-02-01 00:48:55 +05:30
write ( * , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) ' rotation of loadframe:' , &
math_transpose33 ( bc ( loadcase ) % rotation )
2011-12-04 15:31:32 +05:30
print '(a,f12.6)' , 'temperature:' , bc ( loadcase ) % temperature
2012-01-13 21:48:16 +05:30
print '(a,f12.6)' , 'time: ' , bc ( loadcase ) % time
print '(a,i5)' , 'increments: ' , bc ( loadcase ) % incs
2011-12-06 22:28:17 +05:30
print '(a,i5)' , 'output frequency: ' , bc ( loadcase ) % outputfrequency
print '(a,i5)' , 'restart frequency: ' , bc ( loadcase ) % restartfrequency
2011-12-04 15:31:32 +05:30
2012-02-13 23:11:27 +05:30
if ( any ( bc ( loadcase ) % maskStress . eqv . bc ( loadcase ) % maskDeformation ) ) errorID = 831_pInt ! exclusive or masking only
2011-12-04 15:31:32 +05:30
if ( any ( bc ( loadcase ) % maskStress . and . transpose ( bc ( loadcase ) % maskStress ) . and . &
2012-02-10 17:29:59 +05:30
reshape ( [ . false . , . true . , . true . , . true . , . false . , . true . , . true . , . true . , . false . ] , [ 3 , 3 ] ) ) ) &
2012-02-13 23:11:27 +05:30
errorID = 838_pInt ! no rotation is allowed by stress BC
2012-01-26 19:20:00 +05:30
if ( any ( abs ( math_mul33x33 ( bc ( loadcase ) % rotation , math_transpose33 ( bc ( loadcase ) % rotation ) ) &
2012-02-10 17:29:59 +05:30
- math_I3 ) > reshape ( spread ( rotation_tol , 1 , 9 ) , [ 3 , 3 ] ) ) &
2012-01-26 19:20:00 +05:30
. or . abs ( math_det33 ( bc ( loadcase ) % rotation ) ) > 1.0_pReal + rotation_tol ) &
2012-02-13 23:11:27 +05:30
errorID = 846_pInt ! given rotation matrix contains strain
if ( bc ( loadcase ) % time < 0.0_pReal ) errorID = 834_pInt ! negative time increment
if ( bc ( loadcase ) % incs < 1_pInt ) errorID = 835_pInt ! non-positive incs count
if ( bc ( loadcase ) % outputfrequency < 1_pInt ) errorID = 836_pInt ! non-positive result frequency
2011-12-04 15:31:32 +05:30
if ( errorID > 0_pInt ) call IO_error ( error_ID = errorID , ext_msg = loadcase_string )
2011-10-18 20:15:32 +05:30
enddo
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! debugging parameters
debugGeneral = iand ( debug_spectral , debug_spectralGeneral ) > 0_pInt
debugDivergence = iand ( debug_spectral , debug_spectralDivergence ) > 0_pInt
debugRestart = iand ( debug_spectral , debug_spectralRestart ) > 0_pInt
debugFFTW = iand ( debug_spectral , debug_spectralFFTW ) > 0_pInt
2012-02-10 17:29:59 +05:30
2012-01-25 19:57:26 +05:30
!##################################################################################################
2012-02-10 17:29:59 +05:30
! initialization
2012-01-25 19:57:26 +05:30
!##################################################################################################
2012-02-10 17:29:59 +05:30
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 ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 ) ) ; coordinates = 0.0_pReal
allocate ( temperature ( res ( 1 ) , res ( 2 ) , res ( 3 ) ) ) ; temperature = bc ( 1 ) % temperature ! start out isothermally
allocate ( xi ( 3 , res1_red , res ( 2 ) , res ( 3 ) ) ) ; xi = 0.0_pReal
tensorField = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * 9_pInt , C_SIZE_T ) ) ! allocate continous data using a C function, C_SIZE_T is of type integer(8)
call c_f_pointer ( tensorField , tensorField_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , 3 , 3 ] ) ! place a pointer for the real representation
call c_f_pointer ( tensorField , tensorField_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , 3 , 3 ] ) ! place a pointer for the complex representation
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-10 17:29:59 +05:30
! init fields to no deformation
ielem = 0_pInt
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
ielem = ielem + 1_pInt
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
2012-02-16 00:28:38 +05:30
coordinates ( i , j , k , 1 : 3 ) = geomdim / real ( res * [ i , j , k ] , pReal ) - geomdim / real ( 2_pInt * res , pReal )
2012-02-10 17:29:59 +05:30
call CPFEM_general ( 2_pInt , coordinates ( i , j , k , 1 : 3 ) , 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
c0_66 = math_Mandel3333to66 ( c0_reference )
call math_invert ( 6_pInt , c0_66 , s0_66 , i , errmatinv ) ! invert in mandel notation
2012-02-13 23:11:27 +05:30
if ( errmatinv ) call IO_error ( error_ID = 400_pInt )
2012-02-10 17:29:59 +05:30
s0_reference = math_Mandel66to3333 ( s0_66 )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-10 17:29:59 +05:30
! possible restore deformation gradient from saved state
if ( restartInc > 1_pInt ) then ! using old values from file
if ( debugRestart ) print '(a,i6,a)' , 'Reading values of increment ' , &
restartInc - 1_pInt , ' from file'
2012-02-13 23:11:27 +05:30
call IO_read_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , &
trim ( getSolverJobName ( ) ) , size ( defgrad ) )
read ( 777 , rec = 1 ) defgrad
close ( 777 )
2012-02-10 17:29:59 +05:30
defgradold = defgrad
defgradAim = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
defgradAim = defgradAim + defgrad ( i , j , k , 1 : 3 , 1 : 3 ) ! calculating old average deformation
enddo ; enddo ; enddo
defgradAim = defgradAim * wgt
defgradAimOld = defgradAim
endif
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) and remove the given highest frequencies
do k = 1_pInt , res ( 3 )
k_s ( 3 ) = k - 1_pInt
if ( k > res ( 3 ) / 2_pInt + 1_pInt ) k_s ( 3 ) = k_s ( 3 ) - res ( 3 )
do j = 1_pInt , res ( 2 )
k_s ( 2 ) = j - 1_pInt
if ( j > res ( 2 ) / 2_pInt + 1_pInt ) k_s ( 2 ) = k_s ( 2 ) - res ( 2 )
2012-02-16 00:28:38 +05:30
do i = 1_pInt , res1_red
2012-02-10 17:29:59 +05:30
k_s ( 1 ) = i - 1_pInt
xi ( 1 : 3 , i , j , k ) = real ( k_s , pReal ) / geomdim
enddo ; enddo ; enddo
2012-01-25 19:57:26 +05:30
2012-02-10 17:29:59 +05:30
!--------------------------------------------------------------------------------------------------
! calculate the gamma operator
2012-02-16 00:28:38 +05:30
if ( memory_efficient ) then ! allocate just single fourth order tensor
2012-02-10 17:29:59 +05:30
allocate ( gamma_hat ( 1 , 1 , 1 , 3 , 3 , 3 , 3 ) ) ; gamma_hat = 0.0_pReal
2012-02-16 00:28:38 +05:30
else ! precalculation of gamma_hat field
2012-02-10 17:29:59 +05:30
allocate ( gamma_hat ( res1_red , 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 , res1_red
2012-02-13 18:08:46 +05:30
! if(k==res(3)/2 .or. k==res(3)/2+2 .or.&
! j==res(2)/2 .or. j==res(2)/2+2 .or.&
! i==res(1)/2 .or. i==res(1)/2+2) then
! gamma_hat(i,j,k,1:3,1:3,1:3,1:3) = s0_reference
! else
2012-02-16 00:28:38 +05:30
if ( any ( [ i , j , k ] / = 1_pInt ) ) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
xiDyad ( l , m ) = xi ( l , i , j , k ) * xi ( m , i , j , k )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
temp33_Real ( l , m ) = sum ( c0_reference ( l , m , 1 : 3 , 1 : 3 ) * xiDyad )
temp33_Real = math_inv33 ( temp33_Real )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt , p = 1_pInt : 3_pInt ) &
gamma_hat ( i , j , k , l , m , n , p ) = temp33_Real ( l , n ) * xiDyad ( m , p )
endif
2012-02-10 17:29:59 +05:30
enddo ; enddo ; enddo
2012-02-13 18:08:46 +05:30
gamma_hat ( 1 , 1 , 1 , 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 ) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
2012-02-10 17:29:59 +05:30
endif
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! general initialization of fftw (see manual on fftw.org for more details)
2012-02-16 00:28:38 +05:30
if ( pReal / = C_DOUBLE . or . pInt / = C_INT ) call IO_error ( error_ID = 808_pInt ) ! check for correct precision in C
2012-01-25 19:57:26 +05:30
#ifdef _OPENMP
2012-02-10 17:29:59 +05:30
if ( DAMASK_NumThreadsInt > 0_pInt ) then
ierr = fftw_init_threads ( )
2012-02-13 23:11:27 +05:30
if ( ierr == 0_pInt ) call IO_error ( error_ID = 809_pInt )
2012-02-10 17:29:59 +05:30
call fftw_plan_with_nthreads ( DAMASK_NumThreadsInt )
endif
2012-01-25 19:57:26 +05:30
#endif
2012-02-16 00:28:38 +05:30
call fftw_set_timelimit ( fftw_timelimit ) ! set timelimit for plan creation
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! creating plans
2012-02-10 17:29:59 +05:30
plan_stress = fftw_plan_many_dft_r2c ( 3 , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , 9 , & ! dimensions , length in each dimension in reversed order
tensorField_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , & ! input data , physical length in each dimension in reversed order
1 , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , & ! striding , product of physical lenght in the 3 dimensions
tensorField_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
1 , res ( 3 ) * res ( 2 ) * res1_red , fftw_planner_flag )
plan_correction = fftw_plan_many_dft_c2r ( 3 , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , 9 , &
tensorField_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
1 , res ( 3 ) * res ( 2 ) * res1_red , &
tensorField_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
1 , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , fftw_planner_flag )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! depending on (debug) options, allocate more memory and create additional plans
2012-02-10 17:29:59 +05:30
if ( . not . simplified_algorithm ) then
print * , 'using polarization field based algorithm'
tau = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * 9_pInt , C_SIZE_T ) )
call c_f_pointer ( tau , tau_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , 3 , 3 ] )
call c_f_pointer ( tau , tau_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , 3 , 3 ] )
plan_tau = fftw_plan_many_dft_r2c ( 3 , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , 9 , &
tau_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
1 , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , &
tau_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
1 , res ( 3 ) * res ( 2 ) * res1_red , fftw_planner_flag )
endif
2012-01-13 21:48:16 +05:30
2012-02-10 17:29:59 +05:30
if ( debugDivergence ) then
divergence = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * 3_pInt , C_SIZE_T ) )
call c_f_pointer ( divergence , divergence_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , 3 ] )
call c_f_pointer ( divergence , divergence_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , 3 ] )
allocate ( divergence_postProc ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 ) ) ; divergence_postProc = 0.0_pReal
plan_divergence = fftw_plan_many_dft_c2r ( 3 , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , 3 , &
divergence_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
1 , res ( 3 ) * res ( 2 ) * res1_red , &
divergence_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
1 , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , fftw_planner_flag )
endif
2012-01-13 21:48:16 +05:30
2012-02-10 17:29:59 +05:30
if ( debugFFTW ) then
scalarField_realC = fftw_alloc_complex ( int ( res ( 1 ) * res ( 2 ) * res ( 3 ) , C_SIZE_T ) ) ! do not do an inplace transform
scalarField_fourierC = fftw_alloc_complex ( int ( res ( 1 ) * res ( 2 ) * res ( 3 ) , C_SIZE_T ) )
call c_f_pointer ( scalarField_realC , scalarField_real , [ res ( 1 ) , res ( 2 ) , res ( 3 ) ] )
call c_f_pointer ( scalarField_fourierC , scalarField_fourier , [ res ( 1 ) , res ( 2 ) , res ( 3 ) ] )
plan_scalarField_forth = fftw_plan_dft_3d ( res ( 3 ) , res ( 2 ) , res ( 1 ) , & !reversed order
scalarField_real , scalarField_fourier , - 1 , fftw_planner_flag )
plan_scalarField_back = fftw_plan_dft_3d ( res ( 3 ) , res ( 2 ) , res ( 1 ) , & !reversed order
scalarField_fourier , scalarField_real , + 1 , fftw_planner_flag )
endif
if ( debugGeneral ) print '(a)' , 'FFTW initialized'
2012-01-13 21:48:16 +05:30
2012-02-10 17:29:59 +05:30
!--------------------------------------------------------------------------------------------------
! do not correct divergence criterion (usefull to kill dimension and resolution dependenc)
correctionFactor = 1.0_pReal
2011-11-21 23:42:40 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-10 17:29:59 +05:30
! write header of output file
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
write ( 538 ) 'resolution' , res
write ( 538 ) 'dimension' , geomdim
write ( 538 ) 'materialpoint_sizeResults' , materialpoint_sizeResults
write ( 538 ) 'loadcases' , N_Loadcases
write ( 538 ) 'frequencies' , bc ( 1 : N_Loadcases ) % outputfrequency ! one entry per loadcase
write ( 538 ) 'times' , bc ( 1 : N_Loadcases ) % time ! one entry per loadcase
write ( 538 ) 'logscales' , bc ( 1 : N_Loadcases ) % logscale
write ( 538 ) 'increments' , bc ( 1 : N_Loadcases ) % incs ! one entry per loadcase
write ( 538 ) 'startingIncrement' , restartInc - 1_pInt ! start with writing out the previous inc
write ( 538 ) 'eoh' ! end of header
write ( 538 ) materialpoint_results ( 1_pInt : materialpoint_sizeResults , 1 , 1_pInt : Npoints ) ! initial (non-deformed or read-in) results
if ( debugGeneral ) print '(a)' , 'Header of result file written out'
2012-01-25 19:57:26 +05:30
2012-02-10 17:29:59 +05:30
!##################################################################################################
! Loop over loadcases defined in the loadcase file
!##################################################################################################
do loadcase = 1_pInt , N_Loadcases
time0 = time ! loadcase start time
if ( bc ( loadcase ) % followFormerTrajectory . and . &
( restartInc < totalIncsCounter . or . &
restartInc > totalIncsCounter + bc ( loadcase ) % incs ) ) then ! continue to guess along former trajectory where applicable
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc
endif
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-10 17:29:59 +05:30
! arrays for mixed boundary conditions
mask_defgrad = merge ( ones , zeroes , bc ( loadcase ) % maskDeformation )
mask_stress = merge ( ones , zeroes , bc ( loadcase ) % maskStress )
size_reduced = count ( bc ( loadcase ) % maskStressVector )
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-12-23 18:00:35 +05:30
2012-01-25 19:57:26 +05:30
2012-02-10 17:29:59 +05:30
!##################################################################################################
! loop oper incs defined in input file for current loadcase
!##################################################################################################
do inc = 1_pInt , bc ( loadcase ) % incs
totalIncsCounter = totalIncsCounter + 1_pInt
if ( totalIncsCounter > = restartInc ) then ! do calculations (otherwise just forwarding)
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-10 17:29:59 +05:30
! forwarding time
timeinc_old = timeinc
if ( bc ( loadcase ) % logscale == 0_pInt ) then ! linear scale
timeinc = bc ( loadcase ) % time / bc ( loadcase ) % incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
2011-12-06 22:28:17 +05:30
else
2012-02-10 17:29:59 +05:30
if ( loadcase == 1_pInt ) then ! 1st loadcase of logarithmic scale
if ( inc == 1_pInt ) then ! 1st inc of 1st loadcase of logarithmic scale
timeinc = bc ( 1 ) % time * ( 2.0_pReal ** real ( 1_pInt - bc ( 1 ) % incs , pReal ) ) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st loadcase of logarithmic scale
timeinc = bc ( 1 ) % time * ( 2.0_pReal ** real ( inc - 1_pInt - bc ( 1 ) % incs , pReal ) )
endif
else ! not-1st loadcase of logarithmic scale
timeinc = time0 * ( ( 1.0_pReal + bc ( loadcase ) % time / time0 ) ** ( real ( inc , pReal ) / &
real ( bc ( loadcase ) % incs , pReal ) ) &
- ( 1.0_pReal + bc ( loadcase ) % time / time0 ) ** ( real ( ( inc - 1_pInt ) , pReal ) / &
real ( bc ( loadcase ) % incs , pReal ) ) )
endif
2011-11-15 23:24:18 +05:30
endif
2012-02-10 17:29:59 +05:30
time = time + timeinc
2011-12-06 22:28:17 +05:30
2012-02-10 17:29:59 +05:30
if ( bc ( loadcase ) % velGradApplied ) then ! calculate deltaF from given L and current F
deltaF = timeinc * mask_defgrad * math_mul33x33 ( bc ( loadcase ) % deformation , defgradAim )
else ! deltaF = fDot *timeinc where applicable
deltaF = timeinc * mask_defgrad * bc ( loadcase ) % deformation
endif
2011-11-18 03:41:05 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! winding forward of deformation aim in loadcase system
2011-11-18 03:41:05 +05:30
temp33_Real = defgradAim
defgradAim = defgradAim &
2012-02-10 17:29:59 +05:30
+ guessmode * mask_stress * ( defgradAim - defgradAimOld ) * timeinc / timeinc_old &
+ deltaF
2011-11-18 03:41:05 +05:30
defgradAimOld = temp33_Real
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
2012-02-10 17:29:59 +05:30
deltaF = math_rotate_backward33 ( deltaF , bc ( loadcase ) % rotation )
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
temp33_Real = defgrad ( 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
2011-11-18 03:41:05 +05:30
+ guessmode * ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) - defgradold ( i , j , k , 1 : 3 , 1 : 3 ) ) & ! guessing...
2012-02-10 17:29:59 +05:30
* timeinc / timeinc_old &
+ ( 1.0_pReal - guessmode ) * deltaF ! if not guessing, use prescribed average deformation where applicable
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Real
enddo ; enddo ; enddo
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate reduced compliance
2012-01-26 19:20:00 +05:30
c_prev = math_rotate_forward3333 ( c_current * wgt , bc ( loadcase ) % rotation ) ! calculate stiffness from former inc
2012-01-25 19:57:26 +05:30
if ( size_reduced > 0_pInt ) then ! calculate compliance in case stress BC is applied
2011-11-18 03:41:05 +05:30
c_prev99 = math_Plain3333to99 ( c_prev )
2012-01-25 19:57:26 +05:30
k = 0_pInt ! build reduced stiffness
2011-11-18 03:41:05 +05:30
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
2012-01-25 19:57:26 +05:30
call math_invert ( size_reduced , c_reduced , s_reduced , i , errmatinv ) ! invert reduced stiffness
2012-02-13 23:11:27 +05:30
if ( errmatinv ) call IO_error ( error_ID = 400_pInt )
2012-01-25 19:57:26 +05:30
s_prev99 = 0.0_pReal ! build full compliance
2011-11-18 03:41:05 +05:30
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
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! report begin of new increment
2012-02-02 02:00:27 +05:30
print '(a)' , '##################################################################'
2012-02-16 00:28:38 +05:30
print '(A,I5.5,A,es12.5)' , 'Increment ' , totalIncsCounter , ' Time ' , time
2012-02-10 17:29:59 +05:30
2012-01-25 19:57:26 +05:30
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
!##################################################################################################
! convergence loop (looping over iterations)
!##################################################################################################
2011-11-18 03:41:05 +05:30
do while ( iter < itmax . and . &
( err_div > err_div_tol . or . &
err_stress > err_stress_tol ) )
iter = iter + 1_pInt
2011-12-06 22:28:17 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2012-02-02 02:00:27 +05:30
print '(a)' , '=================================================================='
2012-01-31 20:24:49 +05:30
print '(5(a,i6.6))' , 'Loadcase ' , loadcase , ' Increment ' , inc , '/' , bc ( loadcase ) % incs , &
' @ Iteration ' , iter , '/' , itmax
2011-11-18 03:41:05 +05:30
do n = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
2011-12-06 22:28:17 +05:30
defgrad_av_lab ( 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
2012-02-02 02:00:27 +05:30
write ( * , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) 'deformation gradient:' , &
2012-01-31 20:24:49 +05:30
math_transpose33 ( math_rotate_forward33 ( defgrad_av_lab , bc ( loadcase ) % rotation ) )
2011-12-06 22:28:17 +05:30
print '(a)' , ''
2012-02-02 02:00:27 +05:30
print '(a)' , '... update stress field P(F) .....................................'
2011-12-06 22:28:17 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
2012-02-10 17:29:59 +05:30
call deformed_fft ( res , geomdim , defgrad_av_lab , 1.0_pReal , defgrad , coordinates ) ! calculate current coordinates
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
2012-02-10 17:29:59 +05:30
call CPFEM_general ( 3_pInt , & ! collect cycle
2012-01-25 19:57:26 +05:30
coordinates ( i , j , k , 1 : 3 ) , defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) , &
2011-11-18 03:41:05 +05:30
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
2012-01-25 19:57:26 +05:30
tensorField_real = 0.0_pReal ! needed because of the padding for FFTW
2011-11-18 03:41:05 +05:30
c_current = 0.0_pReal
2012-01-13 21:48:16 +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
2012-01-25 19:57:26 +05:30
call CPFEM_general ( CPFEM_mode , & ! first element in first iteration retains CPFEM_mode 1,
coordinates ( i , j , k , 1 : 3 ) , &
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) , & ! others get 2 (saves winding forward effort)
2011-11-18 03:41:05 +05:30
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , &
cstress , dsde , pstress , dPdF )
CPFEM_mode = 2_pInt
2012-01-25 19:57:26 +05:30
tensorField_real ( i , j , k , 1 : 3 , 1 : 3 ) = pstress
2011-11-18 03:41:05 +05:30
c_current = c_current + dPdF
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
2011-12-06 22:28:17 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! copy one component of the stress field to to a single FT and check for mismatch
if ( debugFFTW ) then
row = ( mod ( totalIncsCounter + iter - 2_pInt , 9_pInt ) ) / 3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1
column = ( mod ( totalIncsCounter + iter - 2_pInt , 3_pInt ) ) + 1_pInt
scalarField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) ) = & ! store the selected component
2012-02-16 00:28:38 +05:30
cmplx ( tensorField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , row , column ) , 0.0_pReal , pReal )
2012-01-25 19:57:26 +05:30
endif
2011-12-04 15:31:32 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! build polarization field
if ( . not . simplified_algorithm ) then
tau_real = 0.0_pReal ! padding
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
tau_real ( i , j , k , 1 : 3 , 1 : 3 ) &
= tensorField_real ( i , j , k , 1 : 3 , 1 : 3 ) &
2012-02-02 18:50:09 +05:30
- math_mul3333xx33 ( c0_reference , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) )
2012-01-25 19:57:26 +05:30
enddo ; enddo ; enddo
2012-02-10 17:29:59 +05:30
call fftw_execute_dft_r2c ( plan_tau , tau_real , tau_fourier )
2012-01-25 19:57:26 +05:30
endif
!--------------------------------------------------------------------------------------------------
! call function to calculate divergence from math (for post processing) to check results
if ( debugDivergence ) &
call divergence_fft ( res , geomdim , 3_pInt , &
tensorField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) , divergence_postProc ) !padding
2012-02-02 02:00:27 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-02 02:00:27 +05:30
! doing the FT because it simplifies calculation of average stress in real space also
2012-02-10 17:29:59 +05:30
call fftw_execute_dft_r2c ( plan_stress , tensorField_real , tensorField_fourier )
pstress_av_lab = real ( tensorField_fourier ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) , pReal ) * wgt
2012-02-02 02:00:27 +05:30
pstress_av = math_rotate_forward33 ( pstress_av_lab , bc ( loadcase ) % rotation )
write ( * , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) 'Piola-Kirchhoff stress / MPa:' , &
2012-02-16 00:28:38 +05:30
math_transpose33 ( pstress_av ) / 1.e6_pReal
2012-02-02 02:00:27 +05:30
2012-02-10 17:29:59 +05:30
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 FT results
if ( debugFFTW ) then
call fftw_execute_dft ( plan_scalarField_forth , scalarField_real , scalarField_fourier )
print '(a,i1,1x,i1)' , 'checking FT results of compontent ' , row , column
2012-02-13 22:45:02 +05:30
print '(a,2(es11.4,1x))' , 'max FT relative error ' , &
2012-02-10 17:29:59 +05:30
maxval ( real ( ( scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) - &
tensorField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) , row , column ) ) / &
scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) ) ) , &
maxval ( aimag ( ( scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) - &
tensorField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) , row , column ) ) / &
scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) ) )
endif
!--------------------------------------------------------------------------------------------------
! removing highest frequencies
tensorField_fourier ( res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) &
= cmplx ( 0.0_pReal , 0.0_pReal , pReal )
tensorField_fourier ( 1 : res1_red , res ( 2 ) / 2_pInt + 1_pInt , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) &
= cmplx ( 0.0_pReal , 0.0_pReal , pReal )
if ( res ( 3 ) > 1_pInt ) &
tensorField_fourier ( 1 : res1_red , 1 : res ( 2 ) , res ( 3 ) / 2_pInt + 1_pInt , 1 : 3 , 1 : 3 ) &
= cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-02-02 02:00:27 +05:30
!--------------------------------------------------------------------------------------------------
! stress BC handling
if ( size_reduced > 0_pInt ) then ! calculate stress BC if applied
err_stress = maxval ( abs ( mask_stress * ( pstress_av - bc ( loadcase ) % stress ) ) ) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval ( abs ( pstress_av ) ) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
print '(a)' , ''
print '(a)' , '... correcting deformation gradient to fulfill BCs ...............'
2012-02-13 22:45:02 +05:30
print '(a,f6.2,a,es11.4,a)' , 'error stress = ' , err_stress / err_stress_tol , &
2012-02-13 18:08:46 +05:30
' (' , err_stress , ' Pa)'
2012-02-10 17:29:59 +05:30
defgradAim = defgradAim - math_mul3333xx33 ( s_prev , ( ( pstress_av - bc ( loadcase ) % stress ) ) ) ! residual on given stress components
2012-02-02 02:00:27 +05:30
if ( debugGeneral ) write ( * , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) 'new deformation aim:' , &
math_transpose33 ( defgradAim )
2012-02-13 22:45:02 +05:30
print '(a,1x,es11.4)' , 'determinant of new deformation: ' , math_det33 ( defgradAim )
2012-02-02 02:00:27 +05:30
else
2012-02-02 18:50:09 +05:30
err_stress_tol = + huge ( 1.0_pReal )
2012-02-02 02:00:27 +05:30
endif
2012-02-10 17:29:59 +05:30
defgradAim_lab = math_rotate_backward33 ( defgradAim , bc ( loadcase ) % rotation ) ! boundary conditions from load frame into lab (Fourier) frame
2012-02-02 02:00:27 +05:30
!--------------------------------------------------------------------------------------------------
! actual spectral method
print '(a)' , ''
print '(a)' , '... calculating equilibrium with spectral method .................'
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
2012-02-10 17:29:59 +05:30
pstress_av_L2 = sqrt ( maxval ( math_eigenvalues33 ( math_mul33x33 ( pstress_av_lab , & ! L_2 norm of average stress
math_transpose33 ( pstress_av_lab ) ) ) ) )
2012-01-25 19:57:26 +05:30
err_div_RMS = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 )
do i = 2_pInt , res1_red - 1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
err_div_RMS = err_div_RMS &
2012-02-10 17:29:59 +05:30
+ 2.0_pReal * ( sum ( real ( math_mul33x3_complex ( tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) , & ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , i , j , k ) ) * two_pi_img ) ** 2.0_pReal ) & ! --> sum squared L_2 norm of vector
2012-02-10 17:29:59 +05:30
+ sum ( aimag ( math_mul33x3_complex ( tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) , &
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , i , j , k ) ) * two_pi_img ) ** 2.0_pReal ) )
2012-01-25 19:57:26 +05:30
enddo
2012-02-10 17:29:59 +05:30
err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart
+ sum ( real ( math_mul33x3_complex ( tensorField_fourier ( 1 , j , k , 1 : 3 , 1 : 3 ) , &
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , 1 , j , k ) ) * two_pi_img ) ** 2.0_pReal ) &
2012-02-10 17:29:59 +05:30
+ sum ( aimag ( math_mul33x3_complex ( tensorField_fourier ( 1 , j , k , 1 : 3 , 1 : 3 ) , &
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , 1 , j , k ) ) * two_pi_img ) ** 2.0_pReal ) &
2012-02-10 17:29:59 +05:30
+ sum ( real ( math_mul33x3_complex ( tensorField_fourier ( res1_red , j , k , 1 : 3 , 1 : 3 ) , &
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , res1_red , j , k ) ) * two_pi_img ) ** 2.0_pReal ) &
2012-02-10 17:29:59 +05:30
+ sum ( aimag ( math_mul33x3_complex ( tensorField_fourier ( res1_red , j , k , 1 : 3 , 1 : 3 ) , &
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , res1_red , j , k ) ) * two_pi_img ) ** 2.0_pReal )
2012-01-25 19:57:26 +05:30
enddo ; enddo
2012-02-01 00:48:55 +05:30
err_div_RMS = sqrt ( err_div_RMS ) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
2012-02-10 17:29:59 +05:30
if ( err_div_RMS / pstress_av_L2 * sqrt ( wgt ) * correctionFactor > err_div &
. and . iter > 2_pInt &
. and . err_stress < err_stress_tol ) then
print * , 'Increasing divergence, stopping iterations'
iter = itmax
endif
err_div = err_div_RMS / pstress_av_L2 * sqrt ( wgt ) * correctionFactor ! criterion to stop iterations
2011-12-06 22:28:17 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate additional divergence criteria and report
if ( debugDivergence ) then ! calculate divergence again
err_div_max = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res1_red
2012-02-10 17:29:59 +05:30
temp3_Complex = math_mul33x3_complex ( tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) , &
2012-02-01 00:48:55 +05:30
xi ( 1 : 3 , i , j , k ) ) * two_pi_img
2012-01-25 19:57:26 +05:30
err_div_max = max ( err_div_max , sqrt ( sum ( abs ( temp3_Complex ) ** 2.0_pReal ) ) )
2012-02-10 17:29:59 +05:30
divergence_fourier ( i , j , k , 1 : 3 ) = temp3_Complex ! need divergence NOT squared
2012-01-25 19:57:26 +05:30
enddo ; enddo ; enddo
2011-12-23 18:00:35 +05:30
2012-02-10 17:29:59 +05:30
call fftw_execute_dft_c2r ( plan_divergence , divergence_fourier , divergence_real )
2012-01-25 19:57:26 +05:30
divergence_real = divergence_real * wgt
err_real_div_RMS = 0.0_pReal
2011-12-23 18:00:35 +05:30
err_real_div_max = 0.0_pReal
2012-01-25 19:57:26 +05:30
max_div_error = 0.0_pReal
2011-12-23 18:00:35 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2012-01-25 19:57:26 +05:30
max_div_error = max ( max_div_error , maxval ( ( divergence_real ( i , j , k , 1 : 3 ) &
- divergence_postProc ( i , j , k , 1 : 3 ) ) / divergence_real ( i , j , k , 1 : 3 ) ) )
err_real_div_RMS = err_real_div_RMS + sum ( divergence_real ( i , j , k , 1 : 3 ) ** 2.0_pReal ) ! avg of L_2 norm of div(stress) in real space
err_real_div_max = max ( err_real_div_max , sqrt ( sum ( divergence_real ( i , j , k , 1 : 3 ) ** 2.0_pReal ) ) ) ! maximum of L two norm of div(stress) in real space
2011-12-23 18:00:35 +05:30
enddo ; enddo ; enddo
2012-01-31 01:55:04 +05:30
err_real_div_RMS = sqrt ( wgt * err_real_div_RMS ) ! RMS in real space
err_div_max = err_div_max * sqrt ( wgt )
2012-02-13 22:45:02 +05:30
print '(a,es11.4)' , 'error divergence FT RMS = ' , err_div_RMS
print '(a,es11.4)' , 'error divergence FT max = ' , err_div_max
print '(a,es11.4)' , 'error divergence Real RMS = ' , err_real_div_RMS
print '(a,es11.4)' , 'error divergence Real max = ' , err_real_div_max
print '(a,es11.4)' , 'max deviat. from postProc = ' , max_div_error
2011-12-23 18:00:35 +05:30
endif
2012-02-16 00:28:38 +05:30
print '(a,f6.2,a,es11.4,3a)' , 'error divergence = ' , err_div / err_div_tol , &
' (' , err_div_RMS , ' N/m' , char ( 179 ) , ')'
2012-01-31 01:55:04 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! divergence is calculated from FT(stress), depending on algorithm use field for spectral method
2012-02-10 17:29:59 +05:30
if ( . not . simplified_algorithm ) tensorField_fourier = tau_fourier
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! to the actual spectral method calculation (mechanical equilibrium)
2012-02-16 00:28:38 +05:30
if ( memory_efficient ) then ! memory saving version, on-the-fly calculation of gamma_hat
2012-02-13 18:08:46 +05:30
2012-01-13 21:48:16 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res1_red
2012-02-13 18:08:46 +05:30
! if(k==res(3)/2 .or. k==res(3)/2+2 .or.&
! j==res(2)/2 .or. j==res(2)/2+2 .or.&
! i==res(1)/2 .or. i==res(1)/2+2) then
! forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt)&
! temp33_Complex(m,n) = sum(s0_reference(m,n, 1:3,1:3)* tensorField_fourier(i,j,k,1:3,1:3))
! else
2012-02-16 00:28:38 +05:30
if ( any ( [ i , j , k ] / = 1_pInt ) ) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
xiDyad ( l , m ) = xi ( l , i , j , k ) * xi ( m , i , j , k )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
temp33_Real ( l , m ) = sum ( c0_reference ( l , m , 1 : 3 , 1 : 3 ) * xiDyad )
temp33_Real = math_inv33 ( temp33_Real )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt , p = 1_pInt : 3_pInt ) &
gamma_hat ( 1 , 1 , 1 , l , m , n , p ) = temp33_Real ( l , n ) * xiDyad ( m , p )
forall ( l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
temp33_Complex ( l , m ) = sum ( gamma_hat ( 1 , 1 , 1 , l , m , 1 : 3 , 1 : 3 ) * &
tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) )
tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Complex
endif
2011-11-18 03:41:05 +05:30
enddo ; enddo ; enddo
2012-02-13 18:08:46 +05:30
2012-02-16 00:28:38 +05:30
else ! use precalculated gamma-operator
2012-02-13 18:08:46 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res1_red
forall ( m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt ) &
2012-02-16 00:28:38 +05:30
temp33_Complex ( m , n ) = sum ( gamma_hat ( i , j , k , m , n , 1 : 3 , 1 : 3 ) * &
tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) )
2012-02-13 18:08:46 +05:30
tensorField_fourier ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Complex
2011-11-18 03:41:05 +05:30
enddo ; enddo ; enddo
2012-02-13 18:08:46 +05:30
endif
2012-02-16 00:28:38 +05:30
if ( simplified_algorithm ) then ! do not use the polarization field based algorithm
tensorField_fourier ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) = cmplx ( ( defgrad_av_lab - defgradAim_lab ) & ! assign (negative) average deformation gradient change to zero frequency (real part)
* real ( Npoints , pReal ) , 0.0_pReal , pReal ) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
2012-02-13 18:08:46 +05:30
else
2012-02-16 00:28:38 +05:30
tensorField_fourier ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) = cmplx ( defgradAim_lab * real ( Npoints , pReal ) , & ! assign deformation aim to zero frequency (real part)
0.0_pReal , pReal )
2011-11-18 03:41:05 +05:30
endif
2011-12-06 22:28:17 +05:30
2012-02-10 17:29:59 +05:30
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results
2012-01-25 19:57:26 +05:30
if ( debugFFTW ) then
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res1_red
2012-02-10 17:29:59 +05:30
scalarField_fourier ( i , j , k ) = tensorField_fourier ( i , j , k , row , column )
2012-01-25 19:57:26 +05:30
enddo ; enddo ; enddo
2012-02-16 00:28:38 +05:30
do i = 0_pInt , res ( 1 ) / 2_pInt - 2_pInt !unpack fft data for conj complex symmetric part
2012-02-10 17:29:59 +05:30
m = 1_pInt
do k = 1_pInt , res ( 3 )
n = 1_pInt
do j = 1_pInt , res ( 2 )
scalarField_fourier ( res ( 1 ) - i , j , k ) = conjg ( scalarField_fourier ( 2 + i , n , m ) )
if ( n == 1_pInt ) n = res ( 2 ) + 1_pInt
n = n - 1_pInt
enddo
if ( m == 1_pInt ) m = res ( 3 ) + 1_pInt
m = m - 1_pInt
enddo ; enddo
2012-01-25 19:57:26 +05:30
endif
!--------------------------------------------------------------------------------------------------
! doing the inverse FT
2012-02-10 17:29:59 +05:30
call fftw_execute_dft_c2r ( plan_correction , tensorField_fourier , tensorField_real ) ! back transform of fluct deformation gradient
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results
2012-01-13 21:48:16 +05:30
if ( debugFFTW ) then
2012-02-01 00:48:55 +05:30
print '(a,i1,1x,i1)' , 'checking iFT results of compontent ' , row , column
2012-02-10 17:29:59 +05:30
call fftw_execute_dft ( plan_scalarField_back , scalarField_fourier , scalarField_real )
2012-02-13 22:45:02 +05:30
print '(a,es11.4)' , 'max iFT relative error ' , &
2012-02-10 17:29:59 +05:30
maxval ( ( real ( scalarField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) ) ) - &
2012-01-25 19:57:26 +05:30
tensorField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , row , column ) ) / &
2012-02-10 17:29:59 +05:30
real ( scalarField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) ) ) )
2012-01-25 19:57:26 +05:30
endif
!--------------------------------------------------------------------------------------------------
! calculate some additional output
if ( debugGeneral ) then
maxCorrectionSkew = 0.0_pReal
maxCorrectionSym = 0.0_pReal
temp33_Real = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
maxCorrectionSym = max ( maxCorrectionSym , &
2012-01-26 19:20:00 +05:30
maxval ( math_symmetric33 ( tensorField_real ( i , j , k , 1 : 3 , 1 : 3 ) ) ) )
2012-01-25 19:57:26 +05:30
maxCorrectionSkew = max ( maxCorrectionSkew , &
2012-01-26 19:20:00 +05:30
maxval ( math_skew33 ( tensorField_real ( i , j , k , 1 : 3 , 1 : 3 ) ) ) )
2012-01-25 19:57:26 +05:30
temp33_Real = temp33_Real + tensorField_real ( i , j , k , 1 : 3 , 1 : 3 )
enddo ; enddo ; enddo
2012-02-13 22:45:02 +05:30
print '(a,1x,es11.4)' , 'max symmetrix correction of deformation:' , &
2012-01-25 19:57:26 +05:30
maxCorrectionSym * wgt
2012-02-13 22:45:02 +05:30
print '(a,1x,es11.4)' , 'max skew correction of deformation:' , &
2012-01-25 19:57:26 +05:30
maxCorrectionSkew * wgt
2012-02-13 22:45:02 +05:30
print '(a,1x,es11.4)' , 'max sym/skew of avg correction: ' , &
2012-01-26 19:20:00 +05:30
maxval ( math_symmetric33 ( temp33_real ) ) / &
maxval ( math_skew33 ( temp33_real ) )
2012-01-13 21:48:16 +05:30
endif
2011-12-06 22:28:17 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! updated deformation gradient
2012-02-10 17:29:59 +05:30
defgrad = defgrad - tensorField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) * wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
2011-12-06 22:28:17 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! updated deformation gradient in case of fluctuation field algorithm
if ( . not . simplified_algorithm ) then
defgrad = tensorField_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) * wgt
endif
2012-01-31 01:55:04 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate bounds of det(F) and report
2012-02-02 02:00:27 +05:30
if ( debugGeneral ) then
2012-01-25 19:57:26 +05:30
defgradDetMax = - huge ( 1.0_pReal )
defgradDetMin = + huge ( 1.0_pReal )
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2012-01-26 19:20:00 +05:30
defgradDet = math_det33 ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) )
2012-01-25 19:57:26 +05:30
defgradDetMax = max ( defgradDetMax , defgradDet )
defgradDetMin = min ( defgradDetMin , defgradDet )
enddo ; enddo ; enddo
2011-12-06 22:28:17 +05:30
2012-02-13 22:45:02 +05:30
print '(a,1x,es11.4)' , 'max determinant of deformation:' , defgradDetMax
print '(a,1x,es11.4)' , 'min determinant of deformation:' , defgradDetMin
2012-01-25 19:57:26 +05:30
endif
2011-11-18 03:41:05 +05:30
enddo ! end looping when convergency is achieved
2012-01-04 23:13:26 +05:30
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2012-02-02 02:00:27 +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
2012-01-13 21:48:16 +05:30
print '(A,I5.5,A)' , 'increment ' , totalIncsCounter , ' NOT converged'
2011-12-04 15:31:32 +05:30
notConvergedCounter = notConvergedCounter + 1_pInt
2011-11-18 03:41:05 +05:30
else
2012-02-10 17:29:59 +05:30
convergedCounter = convergedCounter + 1_pInt
2012-01-13 21:48:16 +05:30
print '(A,I5.5,A)' , 'increment ' , totalIncsCounter , ' converged'
2011-11-18 03:41:05 +05:30
endif
2012-02-10 17:29:59 +05:30
2012-01-13 21:48:16 +05:30
if ( mod ( totalIncsCounter - 1_pInt , bc ( loadcase ) % outputfrequency ) == 0_pInt ) then ! at output frequency
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2012-02-02 02:00:27 +05:30
print '(a)' , '... writing results to file ......................................'
2012-02-01 00:48:55 +05:30
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
2012-02-10 17:29:59 +05:30
if ( bc ( loadcase ) % restartFrequency > 0_pInt . and . &
mod ( inc - 1_pInt , bc ( loadcase ) % restartFrequency ) == 0_pInt ) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
restartWrite = . true .
print '(A)' , 'writing converged results for restart'
2012-02-13 23:11:27 +05:30
call IO_write_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , size ( defgrad ) ) ! writing deformation gradient field to file
write ( 777 , rec = 1 ) defgrad
close ( 777 )
2012-02-10 17:29:59 +05:30
restartInc = totalIncsCounter
endif
2012-01-25 19:57:26 +05:30
if ( update_gamma ) then
print * , 'update c0_reference '
c0_reference = c_current * wgt
2012-02-10 17:29:59 +05:30
! s0_reference = math_Plain99to3333(s0_99)
!c0_99 = math_Plain3333to99(c0_reference)
! call math_invert(9_pInt, s0_99, c0_99, i, errmatinv) ! invert reduced stiffness
2012-02-13 23:11:27 +05:30
! if(errmatinv) call IO_error(error_ID=400_pInt)
2012-02-10 17:29:59 +05:30
! print*, (c0_reference - math_Plain99to3333(c0_99))/c0_reference
! pause
2012-01-25 19:57:26 +05:30
endif
2012-01-04 23:13:26 +05:30
2011-12-06 22:28:17 +05:30
endif ! end calculation/forwarding
2012-01-13 21:48:16 +05:30
enddo ! end looping over incs in current loadcase
2011-08-26 19:36:37 +05:30
deallocate ( c_reduced )
deallocate ( s_reduced )
enddo ! end looping over loadcases
2011-12-04 15:31:32 +05:30
print '(a)' , ''
2012-02-02 02:00:27 +05:30
print '(a)' , '##################################################################'
2012-01-16 20:40:16 +05:30
print '(i6.6,a,i6.6,a)' , notConvergedCounter , ' out of ' , &
2012-02-10 17:29:59 +05:30
notConvergedCounter + convergedCounter , ' increments did not converge!'
2011-01-07 18:26:45 +05:30
close ( 538 )
2012-01-25 19:57:26 +05:30
call fftw_destroy_plan ( plan_stress ) ; call fftw_destroy_plan ( plan_correction )
if ( debugDivergence ) call fftw_destroy_plan ( plan_divergence )
2012-01-13 21:48:16 +05:30
if ( debugFFTW ) then
2012-01-25 19:57:26 +05:30
call fftw_destroy_plan ( plan_scalarField_forth )
call fftw_destroy_plan ( plan_scalarField_back )
2012-01-13 21:48:16 +05:30
endif
2012-01-31 01:55:04 +05:30
stop 0
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
!
!********************************************************************
2012-02-16 00:28:38 +05:30
subroutine quit ( stop_id )
2010-06-08 15:38:15 +05:30
use prec
implicit none
2012-02-16 00:28:38 +05:30
integer ( pInt ) , intent ( in ) :: stop_id
print * , stop_id
stop 'abnormal termination of DAMASK_spectral'
2010-06-08 15:38:15 +05:30
2011-05-24 21:27:59 +05:30
end subroutine