2012-03-09 20:52:52 +05:30
! Copyright 2012 Max-Planck-Institut für Eisenforschung GmbH
2011-04-04 19:39:54 +05:30
!
! This file is part of DAMASK,
2012-03-09 20:52:52 +05:30
! the Düsseldorf 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-03-09 20:52:52 +05:30
2011-05-11 22:08:45 +05:30
program DAMASK_spectral
2010-06-08 15:38:15 +05:30
2012-02-23 22:50:57 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
2012-03-09 20:52:52 +05:30
use DAMASK_interface , only : &
DAMASK_interface_init , &
getLoadcaseName , &
getSolverWorkingDirectoryName , &
getSolverJobName , &
getModelName , &
inputFileExtension
use prec , only : &
pInt , &
pReal , &
DAMASK_NaN
use IO , only : &
IO_isBlank , &
IO_open_file , &
IO_stringPos , &
IO_stringValue , &
IO_floatValue , &
IO_intValue , &
IO_error , &
IO_lc , &
IO_read_jobBinaryFile , &
IO_write_jobBinaryFile
use debug , only : &
debug_what , &
debug_spectral , &
debug_levelBasic , &
debug_spectralDivergence , &
debug_spectralRestart , &
2012-04-05 14:47:09 +05:30
debug_spectralFFTW , &
debug_reset , &
debug_info
2012-03-09 20:52:52 +05:30
2010-07-05 21:31:36 +05:30
use math
2012-03-09 20:52:52 +05:30
use CPFEM , only : &
CPFEM_general , &
CPFEM_initAll
use FEsolving , only : &
restartWrite , &
restartInc
use numerics , only : &
err_div_tol , &
err_stress_tolrel , &
rotation_tol , &
itmax , &
itmin , &
memory_efficient , &
update_gamma , &
divergence_correction , &
DAMASK_NumThreadsInt , &
fftw_planner_flag , &
fftw_timelimit
use homogenization , only : &
materialpoint_sizeResults , &
materialpoint_results
2012-01-25 19:57:26 +05:30
!$ use OMP_LIB ! the openMP function library
2012-03-09 20:52:52 +05:30
!##################################################################################################
2012-01-25 19:57:26 +05:30
! variable declaration
!##################################################################################################
2010-06-08 15:38:15 +05:30
implicit none
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 20:52:52 +05:30
! variables related to information from load case and geom file
real ( pReal ) , dimension ( 9 ) :: &
temp_valueVector !> temporarily from loadcase file when reading in tensors
logical , dimension ( 9 ) :: &
temp_maskVector !> temporarily from loadcase file when reading in tensors
2012-01-25 19:57:26 +05:30
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
2012-03-09 20:52:52 +05:30
integer ( pInt ) :: &
headerLength , &
N_l = 0_pInt , &
N_t = 0_pInt , &
N_n = 0_pInt , &
N_Fdot = 0_pInt , &
Npoints , & ! number of Fourier points
homog , & ! homogenization scheme used
res1_red ! to store res(1)/2 +1
character ( len = 1024 ) :: &
line , &
keyword
logical :: &
gotResolution = . false . , &
gotDimension = . false . , &
gotHomogenization = . false .
2012-02-23 22:50:57 +05:30
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
2010-07-01 20:50:06 +05:30
2012-03-09 20:52:52 +05:30
2011-12-04 15:31:32 +05:30
real ( pReal ) :: wgt
2012-02-23 22:50:57 +05:30
real ( pReal ) , dimension ( 3 ) :: geomdim = 0.0_pReal , virt_dim = 0.0_pReal ! physical dimension of volume element per direction
2012-01-25 19:57:26 +05:30
integer ( pInt ) , dimension ( 3 ) :: res = 1_pInt ! resolution (number of Fourier points) in each direction
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-03-09 20:52:52 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
P_av , &
F_aim = math_I3 , &
F_aim_lastInc = math_I3 , &
mask_stress , &
mask_defgrad , &
deltaF_aim , &
F_aim_lab , &
F_aim_lab_lastIter , &
P_av_lab
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dPdF , &
C_ref = 0.0_pReal , &
C = 0.0_pReal , &
S_lastInc , &
C_lastInc ! stiffness and compliance
real ( pReal ) , dimension ( 6 ) :: sigma ! cauchy stress
real ( pReal ) , dimension ( 6 , 6 ) :: dsde
real ( pReal ) , dimension ( 9 , 9 ) :: temp99_Real ! 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-02-23 22:50:57 +05:30
type ( C_PTR ) :: tensorField ! field in real an fourier space
2012-03-09 20:52:52 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: P_real , deltaF_real ! field in real space (pointer)
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: P_fourier , deltaF_fourier ! field in fourier space (pointer)
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: F , F_lastInc
2012-01-25 19:57:26 +05:30
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-02-23 22:50:57 +05:30
type ( C_PTR ) :: plan_stress , plan_correction ! plans for fftw
2012-01-25 19:57:26 +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 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-03-09 20:52:52 +05:30
integer ( pInt ) :: N_Loadcases , loadcase = 0_pInt , 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-02-23 22:50:57 +05:30
real ( pReal ) :: defgradDet
character ( len = 6 ) :: loadcase_string
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
!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-03-19 22:11:55 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: divergence_post
real ( pReal ) :: pstress_av_L2 , err_div_RMS , err_real_div_RMS , err_post_div_RMS , &
err_div_max , err_real_div_max
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! 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
2012-03-09 20:52:52 +05:30
open ( 6 , encoding = 'UTF-8' )
call DAMASK_interface_init
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) ' <<<+- DAMASK_spectral init -+>>>'
write ( 6 , '(a)' ) ' $Id$'
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ' Working Directory: ' , trim ( getSolverWorkingDirectoryName ( ) )
write ( 6 , '(a)' ) ' Solver Job Name: ' , trim ( getSolverJobName ( ) )
write ( 6 , '(a)' ) ''
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! reading the load case file and allocate data structure containing load cases
2012-03-09 20:52:52 +05:30
call IO_open_file ( myUnit , trim ( getLoadcaseName ( ) ) )
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-03-09 20:52:52 +05:30
call IO_error ( error_ID = 837_pInt , ext_msg = trim ( getLoadcaseName ( ) ) ) ! 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 )
2012-03-09 20:52:52 +05:30
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 )
2012-03-30 01:24:31 +05:30
if ( sum ( bc ( 1 : N_Loadcases ) % incs ) > 9000_pInt ) stop 'to many incs' !discuss with Philip, stop code trouble. suggesting warning
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
2012-03-09 20:52:52 +05:30
call IO_open_file ( myUnit , trim ( getModelName ( ) ) / / 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-02-23 22:50:57 +05:30
( mod ( res ( 3 ) , 2_pInt ) / = 0_pInt . and . res ( 3 ) / = 1_pInt ) ) 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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '#############################################################'
write ( 6 , '(a)' ) 'DAMASK spectral:'
write ( 6 , '(a)' ) 'The spectral method boundary value problem solver for'
write ( 6 , '(a)' ) 'the Duesseldorf Advanced Material Simulation Kit'
write ( 6 , '(a)' ) '#############################################################'
write ( 6 , '(a)' ) 'geometry file: ' , trim ( getModelName ( ) ) / / InputFileExtension
write ( 6 , '(a)' ) '============================================================='
write ( 6 , '(a,3(i12 ))' ) 'resolution a b c:' , res
write ( 6 , '(a,3(f12.5))' ) 'dimension x y z:' , geomdim
write ( 6 , '(a,i5)' ) 'homogenization: ' , homog
write ( 6 , '(a)' ) '#############################################################'
write ( 6 , '(a)' ) 'loadcase file: ' , trim ( getLoadcaseName ( ) )
2011-12-04 15:31:32 +05:30
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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) '============================================================='
write ( 6 , '(a,i6)' ) 'loadcase: ' , loadcase
2011-12-04 15:31:32 +05:30
2012-03-09 20:52:52 +05:30
if ( . not . bc ( loadcase ) % followFormerTrajectory ) write ( 6 , '(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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) 'velocity gradient:'
2011-10-18 20:15:32 +05:30
else
2012-03-09 20:52:52 +05:30
write ( 6 , '(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 )
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,f12.6)' ) 'temperature:' , bc ( loadcase ) % temperature
write ( 6 , '(a,f12.6)' ) 'time: ' , bc ( loadcase ) % time
write ( 6 , '(a,i5)' ) 'increments: ' , bc ( loadcase ) % incs
write ( 6 , '(a,i5)' ) 'output frequency: ' , bc ( loadcase ) % outputfrequency
write ( 6 , '(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-23 22:50:57 +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
2012-03-09 20:52:52 +05:30
debugGeneral = iand ( debug_what ( debug_spectral ) , debug_levelBasic ) / = 0
debugDivergence = iand ( debug_what ( debug_spectral ) , debug_spectralDivergence ) / = 0
debugRestart = iand ( debug_what ( debug_spectral ) , debug_spectralRestart ) / = 0
debugFFTW = iand ( debug_what ( debug_spectral ) , debug_spectralFFTW ) / = 0
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-03-09 20:52:52 +05:30
allocate ( F ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 , 3 ) , source = 0.0_pReal )
allocate ( F_lastInc ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 , 3 ) , source = 0.0_pReal )
2012-02-23 22:50:57 +05:30
allocate ( xi ( 3 , res1_red , res ( 2 ) , res ( 3 ) ) , source = 0.0_pReal )
allocate ( coordinates ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 ) , source = 0.0_pReal )
allocate ( temperature ( res ( 1 ) , res ( 2 ) , res ( 3 ) ) , source = bc ( 1 ) % temperature ) ! start out isothermally
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)
2012-03-09 20:52:52 +05:30
call c_f_pointer ( tensorField , P_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , 3 , 3 ] ) ! place a pointer for a real representation on tensorField
call c_f_pointer ( tensorField , deltaF_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , 3 , 3 ] ) ! place a pointer for a real representation on tensorField
call c_f_pointer ( tensorField , P_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , 3 , 3 ] ) ! place a pointer for a complex representation on tensorField
call c_f_pointer ( tensorField , deltaF_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , 3 , 3 ] ) ! place a pointer for a complex representation on tensorField
2012-02-10 17:29:59 +05:30
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-23 22:50:57 +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
2012-03-09 20:52:52 +05:30
P_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , & ! input data , physical length in each dimension in reversed order
2012-02-23 22:50:57 +05:30
1 , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , & ! striding , product of physical lenght in the 3 dimensions
2012-03-09 20:52:52 +05:30
P_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
2012-02-10 17:29:59 +05:30
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 , &
2012-03-09 20:52:52 +05:30
deltaF_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
2012-02-10 17:29:59 +05:30
1 , res ( 3 ) * res ( 2 ) * res1_red , &
2012-03-09 20:52:52 +05:30
deltaF_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
2012-02-10 17:29:59 +05:30
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 ( 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 ] )
2012-03-19 22:11:55 +05:30
allocate ( divergence_post ( res ( 1 ) , res ( 2 ) , res ( 3 ) , 3 ) ) ; divergence_post = 0.0_pReal
2012-02-10 17:29:59 +05:30
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
2012-03-09 20:52:52 +05:30
if ( debugGeneral ) write ( 6 , '(a)' ) 'FFTW initialized'
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
if ( divergence_correction ) then
2012-03-19 18:49:15 +05:30
do i = 1_pInt , 3_pInt
if ( i / = minloc ( geomdim , 1 ) . and . i / = maxloc ( geomdim , 1 ) ) virt_dim = geomdim / geomdim ( i )
enddo
2012-03-09 20:52:52 +05:30
else
virt_dim = geomdim
endif
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 )
do i = 1_pInt , res1_red
k_s ( 1 ) = i - 1_pInt
xi ( 1 : 3 , i , j , k ) = real ( k_s , pReal ) / virt_dim
enddo ; enddo ; enddo
!--------------------------------------------------------------------------------------------------
! get reference material stifness and 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
F ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
F_lastInc ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
coordinates ( i , j , k , 1 : 3 ) = geomdim / real ( res * [ i , j , k ] , pReal ) - geomdim / real ( 2_pInt * res , pReal )
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 , sigma , dsde , P_real ( i , j , k , 1 : 3 , 1 : 3 ) , dPdF )
C = C + dPdF
enddo ; enddo ; enddo
2012-03-19 18:49:15 +05:30
C_ref = C * wgt ! linear reference material stiffness
2012-03-09 20:52:52 +05:30
!--------------------------------------------------------------------------------------------------
! calculate the gamma operator
if ( memory_efficient ) then ! allocate just single fourth order tensor
allocate ( gamma_hat ( 1 , 1 , 1 , 3 , 3 , 3 , 3 ) , source = 0.0_pReal )
else ! precalculation of gamma_hat field
allocate ( gamma_hat ( res1_red , res ( 2 ) , res ( 3 ) , 3 , 3 , 3 , 3 ) , source = 0.0_pReal )
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res1_red
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 ( C_ref ( 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
enddo ; enddo ; enddo
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
endif
!--------------------------------------------------------------------------------------------------
! possible restore deformation gradient from saved state
if ( restartInc > 1_pInt ) then ! using old values from file
if ( debugRestart ) write ( 6 , '(a,i6,a)' ) 'Reading values of increment ' , &
restartInc - 1_pInt , ' from file'
call IO_read_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , &
trim ( getSolverJobName ( ) ) , size ( F ) )
read ( 777 , rec = 1 ) F
close ( 777 )
F_lastInc = F
F_aim = 0.0_pReal
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
F_aim = F_aim + F ( i , j , k , 1 : 3 , 1 : 3 ) ! calculating old average deformation
enddo ; enddo ; enddo
F_aim = F_aim * wgt
F_aim_lastInc = F_aim
endif
2012-01-13 21:48:16 +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
2012-02-23 22:50:57 +05:30
write ( 538 ) 'frequencies' , bc ( 1 : N_Loadcases ) % outputfrequency ! one entry per loadcase
write ( 538 ) 'times' , bc ( 1 : N_Loadcases ) % time ! one entry per loadcase
2012-02-10 17:29:59 +05:30
write ( 538 ) 'logscales' , bc ( 1 : N_Loadcases ) % logscale
2012-02-23 22:50:57 +05:30
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
2012-03-09 20:52:52 +05:30
if ( debugGeneral ) write ( 6 , '(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 )
2012-02-23 22:50:57 +05:30
size_reduced = int ( count ( bc ( loadcase ) % maskStressVector ) , pInt )
allocate ( c_reduced ( size_reduced , size_reduced ) , source = 0.0_pReal )
allocate ( s_reduced ( size_reduced , size_reduced ) , source = 0.0_pReal )
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-03-09 20:52:52 +05:30
if ( bc ( loadcase ) % velGradApplied ) then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * mask_defgrad * math_mul33x33 ( bc ( loadcase ) % deformation , F_aim )
else ! deltaF_aim = fDot *timeinc where applicable
deltaF_aim = timeinc * mask_defgrad * bc ( loadcase ) % deformation
2012-02-10 17:29:59 +05:30
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
2012-03-09 20:52:52 +05:30
temp33_Real = F_aim
F_aim = F_aim &
+ guessmode * mask_stress * ( F_aim - F_aim_lastInc ) * timeinc / timeinc_old &
+ deltaF_aim
F_aim_lastInc = temp33_Real
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 20:52:52 +05:30
! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33 ( deltaF_aim , bc ( loadcase ) % rotation )
2012-02-10 17:29:59 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2012-03-09 20:52:52 +05:30
temp33_Real = F ( i , j , k , 1 : 3 , 1 : 3 )
F ( i , j , k , 1 : 3 , 1 : 3 ) = F ( i , j , k , 1 : 3 , 1 : 3 ) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * ( F ( i , j , k , 1 : 3 , 1 : 3 ) - F_lastInc ( i , j , k , 1 : 3 , 1 : 3 ) ) & ! guessing...
2012-02-10 17:29:59 +05:30
* timeinc / timeinc_old &
2012-03-09 20:52:52 +05:30
+ ( 1.0_pReal - guessmode ) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable
F_lastInc ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Real
2012-02-10 17:29:59 +05:30
enddo ; enddo ; enddo
2012-03-09 20:52:52 +05:30
call deformed_fft ( res , geomdim , math_rotate_backward33 ( F_aim , bc ( loadcase ) % rotation ) , & ! calculate current coordinates
1.0_pReal , F_lastInc , coordinates )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculate reduced compliance
if ( size_reduced > 0_pInt ) then ! calculate compliance in case stress BC is applied
2012-03-09 20:52:52 +05:30
C_lastInc = math_rotate_forward3333 ( C * wgt , bc ( loadcase ) % rotation ) ! calculate stiffness from former inc
temp99_Real = math_Plain3333to99 ( C_lastInc )
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
2012-03-09 20:52:52 +05:30
c_reduced ( k , j ) = temp99_Real ( n , m )
2011-11-18 03:41:05 +05:30
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-03-09 20:52:52 +05:30
temp99_Real = 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
2012-03-09 20:52:52 +05:30
temp99_Real ( n , m ) = s_reduced ( k , j )
2011-11-18 03:41:05 +05:30
endif ; enddo ; endif ; enddo
2012-03-09 20:52:52 +05:30
S_lastInc = ( math_Plain99to3333 ( temp99_Real ) )
2011-11-18 03:41:05 +05:30
endif
2011-11-15 23:24:18 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! report begin of new increment
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) '##################################################################'
write ( 6 , '(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
2012-03-19 18:49:15 +05:30
err_div = huge ( err_div_tol ) ! go into loop
2012-01-25 19:57:26 +05:30
!##################################################################################################
! convergence loop (looping over iterations)
!##################################################################################################
2012-02-23 22:50:57 +05:30
do while ( ( iter < itmax . and . ( err_div > err_div_tol . or . err_stress > err_stress_tol ) ) &
. or . iter < itmin )
2011-11-18 03:41:05 +05:30
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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '=================================================================='
write ( 6 , '(5(a,i6.6))' ) 'Loadcase ' , loadcase , ' Increment ' , inc , '/' , bc ( loadcase ) % incs , &
2012-01-31 20:24:49 +05:30
' @ Iteration ' , iter , '/' , itmax
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) 'deformation gradient aim =' , &
math_transpose33 ( F_aim )
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '... update stress field P(F) .....................................'
F_aim_lab_lastIter = math_rotate_backward33 ( F_aim , bc ( loadcase ) % rotation )
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
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-03-09 20:52:52 +05:30
coordinates ( i , j , k , 1 : 3 ) , F_lastInc ( i , j , k , 1 : 3 , 1 : 3 ) , F ( i , j , k , 1 : 3 , 1 : 3 ) , &
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , sigma , dsde , &
P_real ( i , j , k , 1 : 3 , 1 : 3 ) , dPdF )
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
2011-11-18 03:41:05 +05:30
2012-03-09 20:52:52 +05:30
P_real = 0.0_pReal ! needed because of the padding for FFTW
C = 0.0_pReal
2012-01-13 21:48:16 +05:30
ielem = 0_pInt
2012-04-05 14:47:09 +05:30
call debug_reset ( )
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,
2012-03-09 20:52:52 +05:30
coordinates ( i , j , k , 1 : 3 ) , F_lastInc ( i , j , k , 1 : 3 , 1 : 3 ) , F ( i , j , k , 1 : 3 , 1 : 3 ) , & ! others get 2 (saves winding forward effort)
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , sigma , dsde , &
P_real ( i , j , k , 1 : 3 , 1 : 3 ) , dPdF )
2011-11-18 03:41:05 +05:30
CPFEM_mode = 2_pInt
2012-03-09 20:52:52 +05:30
C = C + dPdF
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
2012-04-05 14:47:09 +05:30
call debug_info ( )
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-03-09 20:52:52 +05:30
cmplx ( P_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
!--------------------------------------------------------------------------------------------------
! call function to calculate divergence from math (for post processing) to check results
if ( debugDivergence ) &
2012-03-19 22:11:55 +05:30
call divergence_fft ( res , virt_dim , 3_pInt , &
P_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) , divergence_post ) ! 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-03-09 20:52:52 +05:30
call fftw_execute_dft_r2c ( plan_stress , P_real , P_fourier )
2012-02-10 17:29:59 +05:30
2012-03-09 20:52:52 +05:30
P_av_lab = real ( P_fourier ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) , pReal ) * wgt
P_av = math_rotate_forward33 ( P_av_lab , bc ( loadcase ) % rotation )
2012-02-23 22:50:57 +05:30
write ( * , '(a,/,3(3(f12.7,1x)/))' , advance = 'no' ) 'Piola-Kirchhoff stress / MPa =' , &
2012-03-09 20:52:52 +05:30
math_transpose33 ( P_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 )
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,i1,1x,i1)' ) 'checking FT results of compontent ' , row , column
write ( 6 , '(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 ) ) - &
2012-03-09 20:52:52 +05:30
P_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) , row , column ) ) / &
2012-02-10 17:29:59 +05:30
scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) ) ) , &
maxval ( aimag ( ( scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) - &
2012-03-09 20:52:52 +05:30
P_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) , row , column ) ) / &
2012-02-10 17:29:59 +05:30
scalarField_fourier ( 1 : res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) ) ) )
endif
!--------------------------------------------------------------------------------------------------
! removing highest frequencies
2012-03-09 20:52:52 +05:30
P_fourier ( res1_red , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) &
2012-02-10 17:29:59 +05:30
= cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-03-09 20:52:52 +05:30
P_fourier ( 1 : res1_red , res ( 2 ) / 2_pInt + 1_pInt , 1 : res ( 3 ) , 1 : 3 , 1 : 3 ) &
2012-02-10 17:29:59 +05:30
= cmplx ( 0.0_pReal , 0.0_pReal , pReal )
if ( res ( 3 ) > 1_pInt ) &
2012-03-09 20:52:52 +05:30
P_fourier ( 1 : res1_red , 1 : res ( 2 ) , res ( 3 ) / 2_pInt + 1_pInt , 1 : 3 , 1 : 3 ) &
2012-02-10 17:29:59 +05:30
= cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-03-09 20:52:52 +05:30
2012-02-02 02:00:27 +05:30
!--------------------------------------------------------------------------------------------------
! stress BC handling
2012-03-09 20:52:52 +05:30
if ( size_reduced > 0_pInt ) then ! calculate stress BC if applied
err_stress = maxval ( abs ( mask_stress * ( P_av - bc ( loadcase ) % stress ) ) ) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval ( abs ( P_av ) ) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '... correcting deformation gradient to fulfill BCs ...............'
write ( 6 , '(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-03-09 20:52:52 +05:30
F_aim = F_aim - math_mul3333xx33 ( S_lastInc , ( ( P_av - bc ( loadcase ) % stress ) ) ) ! residual on given stress components
write ( 6 , '(a,1x,es11.4)' ) 'determinant of new deformation = ' , math_det33 ( F_aim )
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-03-09 20:52:52 +05:30
F_aim_lab = math_rotate_backward33 ( F_aim , bc ( loadcase ) % rotation ) ! boundary conditions from load frame into lab (Fourier) frame
2012-02-02 02:00:27 +05:30
!--------------------------------------------------------------------------------------------------
! actual spectral method
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '... calculating equilibrium with spectral method .................'
2012-02-02 02:00:27 +05:30
2012-01-25 19:57:26 +05:30
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
2012-03-19 22:11:55 +05:30
pstress_av_L2 = sqrt ( maxval ( math_eigenvalues33 ( math_mul33x33 ( P_av_lab , & ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html)
math_transpose33 ( P_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-03-09 20:52:52 +05:30
+ 2.0_pReal * ( sum ( real ( math_mul33x3_complex ( P_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-03-19 18:49:15 +05:30
xi ( 1 : 3 , i , j , k ) ) * TWOPIIMG ) ** 2.0_pReal ) & ! --> sum squared L_2 norm of vector
2012-03-09 20:52:52 +05:30
+ sum ( aimag ( math_mul33x3_complex ( P_fourier ( i , j , k , 1 : 3 , 1 : 3 ) , &
2012-03-19 22:11:55 +05:30
xi ( 1 : 3 , i , j , k ) ) * TWOPIIMG ) ** 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
2012-03-19 22:11:55 +05:30
+ sum ( real ( math_mul33x3_complex ( P_fourier ( 1 , j , k , 1 : 3 , 1 : 3 ) , &
xi ( 1 : 3 , 1 , j , k ) ) * TWOPIIMG ) ** 2.0_pReal ) &
+ sum ( aimag ( math_mul33x3_complex ( P_fourier ( 1 , j , k , 1 : 3 , 1 : 3 ) , &
xi ( 1 : 3 , 1 , j , k ) ) * TWOPIIMG ) ** 2.0_pReal ) &
+ sum ( real ( math_mul33x3_complex ( P_fourier ( res1_red , j , k , 1 : 3 , 1 : 3 ) , &
xi ( 1 : 3 , res1_red , j , k ) ) * TWOPIIMG ) ** 2.0_pReal ) &
+ sum ( aimag ( math_mul33x3_complex ( P_fourier ( res1_red , j , k , 1 : 3 , 1 : 3 ) , &
xi ( 1 : 3 , res1_red , j , k ) ) * TWOPIIMG ) ** 2.0_pReal )
2012-01-25 19:57:26 +05:30
enddo ; enddo
2012-03-19 22:11:55 +05:30
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-03-19 22:11:55 +05:30
if ( err_div_RMS / pstress_av_L2 > err_div &
2012-03-19 18:49:15 +05:30
. and . err_stress < err_stress_tol ) then
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) 'Increasing divergence, stopping iterations'
2012-02-10 17:29:59 +05:30
iter = itmax
endif
2012-03-19 18:49:15 +05:30
err_div = err_div_RMS / pstress_av_L2 ! 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
2012-03-19 22:11:55 +05:30
if ( debugDivergence ) then ! calculate divergence again
2012-01-25 19:57:26 +05:30
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-03-19 18:49:15 +05:30
temp3_Complex = math_mul33x3_complex ( P_fourier ( i , j , k , 1 : 3 , 1 : 3 ) * wgt , & ! weighting P_fourier
2012-03-19 22:11:55 +05:30
xi ( 1 : 3 , i , j , k ) ) * TWOPIIMG
err_div_max = max ( err_div_max , 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
2012-03-19 18:49:15 +05:30
call fftw_execute_dft_c2r ( plan_divergence , divergence_fourier , divergence_real ) ! already weighted
2011-12-23 18:00:35 +05:30
2012-01-25 19:57:26 +05:30
err_real_div_RMS = 0.0_pReal
2012-03-19 22:11:55 +05:30
err_post_div_RMS = 0.0_pReal
2011-12-23 18:00:35 +05:30
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 )
2012-03-19 22:11:55 +05:30
err_real_div_RMS = err_real_div_RMS + sum ( divergence_real ( i , j , k , 1 : 3 ) ** 2.0_pReal ) ! avg of squared L_2 norm of div(stress) in real space
err_post_div_RMS = err_post_div_RMS + sum ( divergence_post ( i , j , k , 1 : 3 ) ** 2.0_pReal ) ! avg of squared L_2 norm of div(stress) in real space
err_real_div_max = max ( err_real_div_max , sum ( divergence_real ( i , j , k , 1 : 3 ) ** 2.0_pReal ) ) ! max of squared L_2 norm of div(stress) in real space
2011-12-23 18:00:35 +05:30
enddo ; enddo ; enddo
2012-03-19 22:11:55 +05:30
2012-02-23 22:50:57 +05:30
err_real_div_RMS = sqrt ( wgt * err_real_div_RMS ) ! RMS in real space
2012-03-19 22:11:55 +05:30
err_post_div_RMS = sqrt ( wgt * err_post_div_RMS ) ! RMS in real space
err_real_div_max = sqrt ( err_real_div_max ) ! max in real space
err_div_max = sqrt ( err_div_max ) ! max in Fourier space
2012-01-31 01:55:04 +05:30
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,es11.4)' ) 'error divergence FT RMS = ' , err_div_RMS
write ( 6 , '(a,es11.4)' ) 'error divergence Real RMS = ' , err_real_div_RMS
2012-03-19 22:11:55 +05:30
write ( 6 , '(a,es11.4)' ) 'error divergence post RMS = ' , err_post_div_RMS
2012-03-19 18:49:15 +05:30
write ( 6 , '(a,es11.4)' ) 'error divergence FT max = ' , err_div_max
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,es11.4)' ) 'error divergence Real max = ' , err_real_div_max
2011-12-23 18:00:35 +05:30
endif
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,f6.2,a,es11.4,a)' ) 'error divergence = ' , err_div / err_div_tol , &
2012-03-30 01:24:31 +05:30
' (' , err_div , ' N/m³)'
2012-02-23 22:50:57 +05:30
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-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 ) &
2012-03-09 20:52:52 +05:30
temp33_Real ( l , m ) = sum ( C_ref ( l , m , 1 : 3 , 1 : 3 ) * xiDyad )
2012-02-16 00:28:38 +05:30
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 ) * &
2012-03-09 20:52:52 +05:30
P_fourier ( i , j , k , 1 : 3 , 1 : 3 ) )
deltaF_fourier ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Complex
2012-02-16 00:28:38 +05:30
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 ) * &
2012-03-09 20:52:52 +05:30
P_fourier ( i , j , k , 1 : 3 , 1 : 3 ) )
deltaF_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-03-09 20:52:52 +05:30
deltaF_fourier ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) = cmplx ( ( F_aim_lab_lastIter - F_aim_lab ) & ! assign (negative) average deformation gradient change to zero frequency (real part)
2012-02-16 00:28:38 +05:30
* real ( Npoints , pReal ) , 0.0_pReal , pReal ) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
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-03-09 20:52:52 +05:30
scalarField_fourier ( i , j , k ) = deltaF_fourier ( i , j , k , row , column )
2012-01-25 19:57:26 +05:30
enddo ; enddo ; enddo
2012-03-09 20:52:52 +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-03-09 20:52:52 +05:30
call fftw_execute_dft_c2r ( plan_correction , deltaF_fourier , deltaF_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-03-09 20:52:52 +05:30
write ( 6 , '(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-03-09 20:52:52 +05:30
write ( 6 , '(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-03-09 20:52:52 +05:30
deltaF_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-03-09 20:52:52 +05:30
maxval ( math_symmetric33 ( deltaF_real ( i , j , k , 1 : 3 , 1 : 3 ) ) ) )
2012-01-25 19:57:26 +05:30
maxCorrectionSkew = max ( maxCorrectionSkew , &
2012-03-09 20:52:52 +05:30
maxval ( math_skew33 ( deltaF_real ( i , j , k , 1 : 3 , 1 : 3 ) ) ) )
temp33_Real = temp33_Real + deltaF_real ( i , j , k , 1 : 3 , 1 : 3 )
2012-01-25 19:57:26 +05:30
enddo ; enddo ; enddo
2012-03-19 22:11:55 +05:30
write ( 6 , '(a,1x,es11.4)' ) 'max symmetric correction of deformation =' , &
2012-01-25 19:57:26 +05:30
maxCorrectionSym * wgt
2012-03-09 20:52:52 +05:30
write ( 6 , '(a,1x,es11.4)' ) 'max skew correction of deformation =' , &
2012-01-25 19:57:26 +05:30
maxCorrectionSkew * wgt
2012-03-09 20:52:52 +05:30
write ( 6 , '(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-03-09 20:52:52 +05:30
F = F - deltaF_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
!--------------------------------------------------------------------------------------------------
! 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-03-09 20:52:52 +05:30
defgradDet = math_det33 ( F ( 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-03-09 20:52:52 +05:30
write ( 6 , '(a,1x,es11.4)' ) 'max determinant of deformation =' , defgradDetMax
write ( 6 , '(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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '=================================================================='
2011-12-04 15:31:32 +05:30
if ( err_div > err_div_tol . or . err_stress > err_stress_tol ) then
2012-03-09 20:52:52 +05:30
write ( 6 , '(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-03-09 20:52:52 +05:30
write ( 6 , '(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-02-23 22:50:57 +05:30
if ( mod ( inc , bc ( loadcase ) % outputFrequency ) == 0_pInt ) then ! at output frequency
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '... writing results to file ......................................'
2012-02-23 22:50:57 +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 . &
2012-02-23 22:50:57 +05:30
mod ( inc , 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?)
2012-02-10 17:29:59 +05:30
restartWrite = . true .
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) 'writing converged results for restart'
call IO_write_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , size ( F ) ) ! writing deformation gradient field to file
write ( 777 , rec = 1 ) F
2012-02-13 23:11:27 +05:30
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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) 'update C_ref '
C_ref = C * wgt
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
2012-03-09 20:52:52 +05:30
write ( 6 , '(a)' ) ''
write ( 6 , '(a)' ) '##################################################################'
write ( 6 , '(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-02-23 22:50:57 +05:30
call quit ( 0_pInt )
2011-05-11 22:08:45 +05:30
end program DAMASK_spectral
2010-06-10 20:21:10 +05:30
2012-04-10 19:00:34 +05:30
#include "DAMASK_quit.f90"