2011-08-26 19:36:37 +05:30
! Copyright 2011 Max-Planck-Institut fuer Eisenforschung GmbH
2011-04-04 19:39:54 +05:30
!
! This file is part of DAMASK,
2011-11-04 01:02:11 +05:30
! the Duesseldorf Advanced Material Simulation Kit.
2011-04-04 19:39:54 +05:30
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
2010-06-08 15:40:57 +05:30
!* $Id$
2010-06-08 15:38:15 +05:30
!********************************************************************
! Material subroutine for BVP solution using spectral method
!
! written by P. Eisenlohr,
! F. Roters,
! L. Hantcherli,
2011-02-07 20:05:42 +05:30
! W.A. Counts,
! D.D. Tjahjanto,
! C. Kords,
! M. Diehl,
2010-06-08 15:38:15 +05:30
! R. Lebensohn
!
! MPI fuer Eisenforschung, Duesseldorf
!
!********************************************************************
! Usage:
2011-10-18 20:15:32 +05:30
! - start program with DAMASK_spectral
! -g (--geom, --geometry) PathToGeomFile/NameOfGeom.geom
! -l (--load, --loadcase) PathToLoadFile/NameOfLoadFile.load
2011-02-07 20:05:42 +05:30
! - PathToGeomFile will be the working directory
2010-06-08 15:38:15 +05:30
! - make sure the file "material.config" exists in the working
2011-11-15 23:24:18 +05:30
! directory. For further configuration use "numerics.config" and
! "numerics.config"
2010-06-10 20:21:10 +05:30
!********************************************************************
2011-05-11 22:08:45 +05:30
program DAMASK_spectral
2010-06-10 20:21:10 +05:30
!********************************************************************
2010-06-08 15:38:15 +05:30
2011-05-11 22:31:03 +05:30
use DAMASK_interface
2010-06-10 20:21:10 +05:30
use prec , only : pInt , pReal
use IO
2011-11-15 23:24:18 +05:30
use debug , only : debug_verbosity , spectral_debug_verbosity
2010-07-05 21:31:36 +05:30
use math
2011-05-26 14:53:13 +05:30
use mesh , only : mesh_ipCenterOfGravity
2011-11-11 19:47:43 +05:30
use CPFEM , only : CPFEM_general , CPFEM_initAll
use FEsolving , only : restartWrite , restartReadSpectral , restartReadStep
2011-10-25 19:08:24 +05:30
use numerics , only : err_div_tol , err_stress_tol , err_stress_tolrel , rotation_tol , &
2011-11-15 23:24:18 +05:30
itmax , memory_efficient , DAMASK_NumThreadsInt , divergence_correction , &
2011-10-25 19:08:24 +05:30
fftw_planner_flag , fftw_timelimit
2010-10-27 22:45:49 +05:30
use homogenization , only : materialpoint_sizeResults , materialpoint_results
dummy update because those messages got lost:
Date: Tue, 13 Sep 2011 17:46:44 +0200
From: m.diehl@mpie.de
To: wangleyu@msu.edu, lebenso@lanl.gov, denny.tjahjanto@imdea.org,
o.guevenc@mpie.de, n.jia@mpie.de, m.diehl@mpie.de, c.kords@mpie.de,
c.zambaldi@mpie.de, p.eisenlohr@mpie.de, f.roters@mpie.de
Subject: update: /home/svn/repos/DAMASK to 1000
Message-ID: <4e6f7ae4.Xjnd/szYAh8QCXTo%m.diehl@mpie.de>
User-Agent: Heirloom mailx 12.2 01/07/07
MIME-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
A processing/pre/FromEBSD/
A processing/pre/FromEBSD/Hex2Cub.cpp
A processing/pre/FromEBSD/SpectralMethodFromEBDS
A processing/pre/FromEBSD/patchFromReconstructedBoundaries
D processing/pre/patchFromReconstructedBoundaries
added two small (quick and dirty) tools to convert data from EBSD to input files for spectral method, put them together with patchFromReconstructedBoundaries into new folder.
Date: Tue, 13 Sep 2011 17:54:06 +0200
From: m.diehl@mpie.de
To: wangleyu@msu.edu, lebenso@lanl.gov, denny.tjahjanto@imdea.org,
o.guevenc@mpie.de, n.jia@mpie.de, m.diehl@mpie.de, c.kords@mpie.de,
c.zambaldi@mpie.de, p.eisenlohr@mpie.de, f.roters@mpie.de
Subject: update: /home/svn/repos/DAMASK to 1001
Message-ID: <4e6f7c9e.v9E4JVN2a6bg5tL8%m.diehl@mpie.de>
User-Agent: Heirloom mailx 12.2 01/07/07
MIME-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
U code/DAMASK_spectral.f90
U code/DAMASK_spectral_interface.f90
U code/IO.f90
U code/crystallite.f90
U code/homogenization_RGC.f90
U code/lattice.f90
U code/makefile
U code/material.f90
U code/mesh.f90
did a lot of polishing:
- removed unnecessary "return" before end of subroutine or function:
- changed undetermined array length (:) to (1:3)
To prevent problems with some code analysing tools:
- "3D oneliner loops" (with ";) only for "do" and "enddo" at the same time
- removed line continuation in OMP statements
made the makefile more flexible, removed heap-arrays switch
Date: Tue, 13 Sep 2011 17:57:07 +0200
From: m.diehl@mpie.de
To: wangleyu@msu.edu, lebenso@lanl.gov, denny.tjahjanto@imdea.org,
o.guevenc@mpie.de, n.jia@mpie.de, m.diehl@mpie.de, c.kords@mpie.de,
c.zambaldi@mpie.de, p.eisenlohr@mpie.de, f.roters@mpie.de
Subject: update: /home/svn/repos/DAMASK to 1002
Message-ID: <4e6f7d53.IEDDzo+JSzDWNSBr%m.diehl@mpie.de>
User-Agent: Heirloom mailx 12.2 01/07/07
MIME-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
A documentation/Compiling/
A documentation/Compiling/Stack+usage.pdf
A documentation/ParallelizationAndTuning/
A documentation/ParallelizationAndTuning/BSC_tools_Overview.pdf
A documentation/ParallelizationAndTuning/Intro_Perf.pdf
A documentation/ParallelizationAndTuning/Kcachegrind.pdf
A documentation/ParallelizationAndTuning/LRZ210703_1.pdf
A documentation/ParallelizationAndTuning/LRZ210703_2.pdf
A documentation/ParallelizationAndTuning/MUST_Overview.pdf
A documentation/ParallelizationAndTuning/NPB-MZ-MPI-BT_Exercise.pdf
A documentation/ParallelizationAndTuning/PAPI.pdf
A documentation/ParallelizationAndTuning/PSC_Exercise_BT-MPI.pdf
A documentation/ParallelizationAndTuning/Paraver_Exercise.pdf
A documentation/ParallelizationAndTuning/Periscope_Overview.pdf
A documentation/ParallelizationAndTuning/SIONlib.pdf
A documentation/ParallelizationAndTuning/Scalasca_Examples.pdf
A documentation/ParallelizationAndTuning/Scalasca_Exercise_BTMZ.pdf
A documentation/ParallelizationAndTuning/Scalasca_Overview.pdf
A documentation/ParallelizationAndTuning/Scalasca_Patterns.pdf
A documentation/ParallelizationAndTuning/TAU.pdf
A documentation/ParallelizationAndTuning/VIHPS-TW8.pdf
A documentation/ParallelizationAndTuning/Vampir_Exercise.pdf
A documentation/ParallelizationAndTuning/Vampir_Overview.pdf
A documentation/ParallelizationAndTuning/instructions_periscope.pdf
A documentation/ParallelizationAndTuning/manualf06.pdf
added some information from Tuning workshop in Aachen regarding tuning/parallelization
added slides with information how to prevent segmentation fauld
2011-09-14 13:46:42 +05:30
!$ use OMP_LIB ! the openMP function library
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
2010-06-08 15:38:15 +05:30
implicit none
2011-01-07 18:26:45 +05:30
! variables to read from loadcase and geom file
2011-11-15 23:24:18 +05:30
real ( pReal ) , dimension ( 9 ) :: temp_valueVector ! stores information temporarily from loadcase file
logical , dimension ( 9 ) :: temp_maskVector
2011-10-18 20:15:32 +05:30
integer ( pInt ) , parameter :: maxNchunksLoadcase = &
2011-10-24 23:56:34 +05:30
( 1_pInt + 9_pInt ) * 3_pInt + & ! deformation, rotation, and stress
2011-11-07 16:34:57 +05:30
( 1_pInt + 1_pInt ) * 5_pInt + & ! time, (log)incs, temp, restartfrequency, and outputfrequency
2011-10-24 23:56:34 +05:30
1_pInt ! dropguessing
2011-10-18 20:15:32 +05:30
integer ( pInt ) , dimension ( 1 + maxNchunksLoadcase * 2 ) :: posLoadcase
2011-10-25 19:08:24 +05:30
integer ( pInt ) , parameter :: maxNchunksGeom = 7_pInt ! 4 identifiers, 3 values
2011-10-18 20:15:32 +05:30
integer ( pInt ) , dimension ( 1 + maxNchunksGeom * 2 ) :: posGeom
2011-11-15 23:24:18 +05:30
integer ( pInt ) :: headerLength , N_l = 0_pInt , N_t = 0_pInt , N_n = 0_pInt , N_Fdot = 0_pInt
2011-11-07 23:55:10 +05:30
integer ( pInt ) , parameter :: myUnit = 234_pInt
2011-10-18 20:15:32 +05:30
character ( len = 1024 ) :: path , line , keyword
2011-11-07 23:55:10 +05:30
logical :: gotResolution = . false . , gotDimension = . false . , gotHomogenization = . false .
2011-11-15 23:24:18 +05:30
type bc_type
real ( pReal ) , dimension ( 3 , 3 ) :: deformation , & ! applied velocity gradient or time derivative of deformation gradient
stress , & ! stress BC (if applicable)
rotation ! rotation of BC (if applicable)
real ( pReal ) :: timeIncrement , & ! length of increment
temperature ! isothermal starting conditions
integer ( pInt ) :: steps , & ! number of steps
outputfrequency , & ! frequency of result writes
restartfrequency , & ! frequency of result writes
logscale ! linear/logaritmic time step flag
logical :: followFormerTrajectory , & ! follow trajectory of former loadcase
velGradApplied ! decide wether velocity gradient or fdot is given
logical , dimension ( 3 , 3 ) :: maskDeformation , & ! mask of boundary conditions
maskStress
logical , dimension ( 9 ) :: maskStressVector ! linear mask of boundary conditions
end type
type ( bc_type ) , allocatable , dimension ( : ) :: bc
character ( len = 3 ) :: loadcase_string
2010-07-01 20:50:06 +05:30
2011-01-07 18:26:45 +05:30
! variables storing information from geom file
2011-10-18 20:15:32 +05:30
real ( pReal ) :: wgt
2011-11-07 23:55:10 +05:30
real ( pReal ) , dimension ( 3 ) :: geomdimension = 0.0_pReal ! physical dimension of volume element in each direction
integer ( pInt ) :: homog ! homogenization scheme used
integer ( pInt ) , dimension ( 3 ) :: resolution = 1_pInt ! resolution (number of Fourier points) in each direction
logical :: spectralPictureMode = . false . ! indicating 1 to 1 mapping of FP to microstructure
2010-07-01 20:50:06 +05:30
2011-11-07 23:55:10 +05:30
! stress, stiffness and compliance average etc.
real ( pReal ) , dimension ( 3 , 3 ) :: pstress , pstress_av , defgrad_av , &
defgradAim = math_I3 , defgradAimOld = math_I3 , defgradAimCorr = math_I3 , &
mask_stress , mask_defgrad , fDot , &
2011-11-11 19:47:43 +05:30
pstress_av_load , defgradAim_lab ! quantities rotated to other coordinate system
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dPdF , c0_reference , c_current = 0.0_pReal , s_prev , c_prev ! stiffness and compliance
real ( pReal ) , dimension ( 6 ) :: cstress ! cauchy stress
real ( pReal ) , dimension ( 6 , 6 ) :: dsde ! small strain stiffness
real ( pReal ) , dimension ( 9 , 9 ) :: s_prev99 , c_prev99 ! compliance and stiffness in matrix notation
real ( pReal ) , dimension ( : , : ) , allocatable :: s_reduced , c_reduced ! reduced compliance and stiffness (only for stress BC)
2011-11-07 23:55:10 +05:30
integer ( pInt ) :: size_reduced = 0.0_pReal ! number of stress BCs
2011-10-18 20:15:32 +05:30
! pointwise data
2011-11-15 23:24:18 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: workfft , defgrad , defgradold
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: coordinates
real ( pReal ) , dimension ( : , : , : ) , allocatable :: temperature
2011-10-18 20:15:32 +05:30
2011-10-25 19:08:24 +05:30
! variables storing information for spectral method and FFTW
2011-08-26 19:36:37 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: xiDyad ! product of wave vectors
real ( pReal ) , dimension ( : , : , : , : , : , : , : ) , allocatable :: gamma_hat ! gamma operator (field) for spectral method
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: xi ! wave vector field
integer ( pInt ) , dimension ( 3 ) :: k_s
2011-11-15 23:24:18 +05:30
integer * 8 , dimension ( 3 ) :: fftw_plan ! plans for fftw (forward and backward)
2011-10-18 20:15:32 +05:30
integer * 8 :: fftw_flag ! planner flag for fftw
2011-02-09 23:17:28 +05:30
2011-02-07 20:05:42 +05:30
! loop variables, convergence etc.
2011-11-11 19:47:43 +05:30
real ( pReal ) :: time = 0.0_pReal , time0 = 0.0_pReal , timeinc ! elapsed time, begin of interval, time interval
2011-10-18 20:15:32 +05:30
real ( pReal ) :: guessmode , err_div , err_stress , p_hat_avg
2011-11-15 23:24:18 +05:30
complex ( pReal ) , parameter :: img = cmplx ( 0.0_pReal , 1.0_pReal )
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: ones = 1.0_pReal , zeroes = 0.0_pReal
complex ( pReal ) , dimension ( 3 , 3 ) :: temp33_Complex
real ( pReal ) , dimension ( 3 , 3 ) :: temp33_Real
2011-10-18 20:15:32 +05:30
integer ( pInt ) :: i , j , k , l , m , n , p
2011-11-15 23:24:18 +05:30
integer ( pInt ) :: N_Loadcases , loadcase , step , iter , ielem , CPFEM_mode , stepZero = 1_pInt , &
2011-11-11 19:47:43 +05:30
ierr , notConvergedCounter = 0_pInt , totalStepsCounter = 0_pInt
2011-11-15 23:24:18 +05:30
logical :: errmatinv , regrid = . false .
2011-11-11 19:47:43 +05:30
real ( pReal ) :: defgradDet , defgradDetMax , defgradDetMin
2011-11-15 23:24:18 +05:30
real ( pReal ) :: correctionFactor
! debuging variables
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: divergence
real ( pReal ) :: p_real_avg , err_div_max , err_real_div_avg , err_real_div_max
logical :: debugGeneral = . false . , debugDivergence = . false . , debugRestart = . false .
2011-07-25 22:00:21 +05:30
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
!Initializing
2011-11-15 23:24:18 +05:30
!$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
if ( . not . ( command_argument_count ( ) == 4 . or . command_argument_count ( ) == 6 ) ) call IO_error ( error_ID = 102_pInt ) ! check for correct number of given arguments
2011-08-26 19:36:37 +05:30
2011-11-04 01:02:11 +05:30
call DAMASK_interface_init ( )
2011-08-02 19:28:28 +05:30
2011-11-15 23:24:18 +05:30
if ( iand ( spectral_debug_verbosity , 1_pInt ) == 1_pInt ) debugGeneral = . true .
if ( iand ( spectral_debug_verbosity , 2_pInt ) == 2_pInt ) debugDivergence = . true .
if ( iand ( spectral_debug_verbosity , 4_pInt ) == 4_pInt ) debugRestart = . true .
2011-10-18 20:15:32 +05:30
!$OMP CRITICAL (write2out)
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-11-15 23:24:18 +05:30
print '(a,a)' , ' <<<+- DAMASK_spectral init -+>>>'
print '(a,a)' , ' $Id$'
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-11-15 23:24:18 +05:30
print '(a,a)' , ' Working Directory: ' , trim ( getSolverWorkingDirectoryName ( ) )
print '(a,a)' , ' Solver Job Name: ' , trim ( getSolverJobName ( ) )
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-10-18 20:15:32 +05:30
!$OMP END CRITICAL (write2out)
2011-11-04 01:02:11 +05:30
! Reading the loadcase file and allocate variables
path = getLoadcaseName ( )
2011-11-15 23:24:18 +05:30
if ( . not . IO_open_file ( myUnit , path ) ) call IO_error ( error_ID = 30_pInt , ext_msg = trim ( path ) )
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
2010-06-10 20:21:10 +05:30
do
2011-10-18 20:15:32 +05:30
read ( myUnit , '(a1024)' , END = 100 ) line
2011-11-04 01:02:11 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2011-10-18 20:15:32 +05:30
posLoadcase = IO_stringPos ( line , maxNchunksLoadcase )
2011-11-15 23:24:18 +05:30
do i = 1_pInt , maxNchunksLoadcase , 1_pInt ! reading compulsory parameters for loadcase
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posLoadcase , i ) ) )
2011-10-25 19:08:24 +05:30
case ( 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' )
2011-11-15 23:24:18 +05:30
N_l = N_l + 1_pInt
2011-07-13 22:03:12 +05:30
case ( 'fdot' )
2011-11-15 23:24:18 +05:30
N_Fdot = N_Fdot + 1_pInt
2011-07-13 22:03:12 +05:30
case ( 't' , 'time' , 'delta' )
2011-11-15 23:24:18 +05:30
N_t = N_t + 1_pInt
2011-07-13 22:03:12 +05:30
case ( 'n' , 'incs' , 'increments' , 'steps' , 'logincs' , 'logsteps' )
2011-11-15 23:24:18 +05:30
N_n = N_n + 1_pInt
2010-06-10 20:21:10 +05:30
end select
2011-11-04 01:02:11 +05:30
enddo ! count all identifiers to allocate memory and do sanity check
2010-06-10 20:21:10 +05:30
enddo
2010-06-10 14:20:04 +05:30
2011-10-18 20:15:32 +05:30
100 N_Loadcases = N_n
2011-11-04 01:02:11 +05:30
if ( ( N_l + N_Fdot / = N_n ) . or . ( N_n / = N_t ) ) & ! sanity check
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 37_pInt , ext_msg = trim ( path ) ) ! error message for incomplete loadcase
allocate ( bc ( N_Loadcases ) )
2010-06-10 20:21:10 +05:30
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
2011-07-14 15:07:31 +05:30
loadcase = 0_pInt
2010-06-10 20:21:10 +05:30
do
2011-10-18 20:15:32 +05:30
read ( myUnit , '(a1024)' , END = 101 ) line
2011-11-15 23:24:18 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
loadcase = loadcase + 1_pInt
bc ( loadcase ) % deformation = zeroes ; bc ( loadcase ) % stress = zeroes ; bc ( loadcase ) % rotation = zeroes
bc ( loadcase ) % timeIncrement = 0.0_pReal ; bc ( loadcase ) % temperature = 30 0.0_pReal
bc ( loadcase ) % steps = 0_pInt ; bc ( loadcase ) % logscale = 0_pInt
bc ( loadcase ) % outputfrequency = 1_pInt ; bc ( loadcase ) % restartfrequency = 1_pInt
bc ( loadcase ) % maskDeformation = . false . ; bc ( loadcase ) % maskStress = . false .
bc ( loadcase ) % maskStressVector = . false . ; bc ( loadcase ) % velGradApplied = . false .
bc ( loadcase ) % followFormerTrajectory = . true .
bc ( loadcase ) % rotation = math_I3 ! assume no rotation, overwrite later in case rotation of loadcase is given
2011-10-18 20:15:32 +05:30
posLoadcase = IO_stringPos ( line , maxNchunksLoadcase )
2011-11-15 23:24:18 +05:30
do j = 1_pInt , maxNchunksLoadcase
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) )
2011-11-15 23:24:18 +05:30
case ( 'fdot' , 'l' , 'velocitygrad' , 'velgrad' , 'velocitygradient' ) ! assign values for the deformation BC matrix
bc ( loadcase ) % velGradApplied = ( IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'l' . or . & ! in case of given L, set flag to true
IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'velocitygrad' . or . &
IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'velgrad' . or . &
IO_lc ( IO_stringValue ( line , posLoadcase , j ) ) == 'velocitygradient' )
temp_valueVector = 0.0_pReal
temp_maskVector = . false .
forall ( k = 1_pInt : 9_pInt ) temp_maskVector ( k ) = IO_stringValue ( line , posLoadcase , j + k ) / = '*'
do k = 1_pInt , 9_pInt
if ( temp_maskVector ( k ) ) temp_valueVector ( k ) = IO_floatValue ( line , posLoadcase , j + k )
2010-06-10 20:21:10 +05:30
enddo
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % maskDeformation = transpose ( reshape ( temp_maskVector , ( / 3 , 3 / ) ) )
bc ( loadcase ) % deformation = math_plain9to33 ( temp_valueVector )
2011-11-04 01:02:11 +05:30
case ( 'p' , 'pk1' , 'piolakirchhoff' , 'stress' )
2011-11-15 23:24:18 +05:30
temp_valueVector = 0.0_pReal
forall ( k = 1_pInt : 9_pInt ) bc ( loadcase ) % maskStressVector ( k ) = IO_stringValue ( line , posLoadcase , j + k ) / = '*'
do k = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( k ) ) temp_valueVector ( k ) = IO_floatValue ( line , posLoadcase , j + k ) ! assign values for the bc(loadcase)%stress matrix
2010-06-10 20:21:10 +05:30
enddo
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % maskStress = transpose ( reshape ( bc ( loadcase ) % maskStressVector , ( / 3 , 3 / ) ) )
bc ( loadcase ) % stress = math_plain9to33 ( temp_valueVector )
case ( 't' , 'time' , 'delta' ) ! increment time
bc ( loadcase ) % timeIncrement = IO_floatValue ( line , posLoadcase , j + 1_pInt )
case ( 'temp' , 'temperature' ) ! starting temperature
bc ( loadcase ) % temperature = IO_floatValue ( line , posLoadcase , j + 1_pInt )
case ( 'n' , 'incs' , 'increments' , 'steps' ) ! steps
bc ( loadcase ) % steps = IO_intValue ( line , posLoadcase , j + 1_pInt )
case ( 'logincs' , 'logsteps' ) ! true, if log scale
bc ( loadcase ) % steps = IO_intValue ( line , posLoadcase , j + 1_pInt )
bc ( loadcase ) % logscale = 1_pInt
case ( 'f' , 'freq' , 'frequency' , 'outputfreq' ) ! frequency of result writings
bc ( loadcase ) % outputfrequency = IO_intValue ( line , posLoadcase , j + 1_pInt )
case ( 'r' , 'restart' , 'restartwrite' ) ! frequency of writing restart information
bc ( loadcase ) % restartfrequency = IO_intValue ( line , posLoadcase , j + 1_pInt )
2011-07-14 15:07:31 +05:30
case ( 'guessreset' , 'dropguessing' )
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % followFormerTrajectory = . false . ! do not continue to predict deformation along former trajectory
case ( 'euler' ) ! rotation of loadcase given in euler angles
p = 0_pInt ! assuming values given in radians
l = 1_pInt ! assuming keyword indicating degree/radians
select case ( IO_lc ( IO_stringValue ( line , posLoadcase , j + 1_pInt ) ) )
2011-10-18 20:15:32 +05:30
case ( 'deg' , 'degree' )
2011-11-15 23:24:18 +05:30
p = 1_pInt ! for conversion from degree to radian
2011-10-18 20:15:32 +05:30
case ( 'rad' , 'radian' )
case default
2011-11-15 23:24:18 +05:30
l = 0_pInt ! imediately reading in angles, assuming radians
2011-10-18 20:15:32 +05:30
end select
2011-11-15 23:24:18 +05:30
forall ( k = 1_pInt : 3_pInt ) temp33_Real ( k , 1 ) = IO_floatValue ( line , posLoadcase , j + l + k ) * real ( p , pReal ) * inRad
bc ( loadcase ) % rotation = math_EulerToR ( temp33_Real ( : , 1 ) )
case ( 'rotation' , 'rot' ) ! assign values for the rotation of loadcase matrix
temp_valueVector = 0.0_pReal
forall ( k = 1_pInt : 9_pInt ) temp_valueVector ( k ) = IO_floatValue ( line , posLoadcase , j + k )
bc ( loadcase ) % rotation = math_plain9to33 ( temp_valueVector )
2010-06-10 20:21:10 +05:30
end select
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
enddo ; enddo
2011-10-18 20:15:32 +05:30
101 close ( myUnit )
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
2011-01-07 18:26:45 +05:30
!read header of geom file to get the information needed before the complete geom file is intepretated by mesh.f90
2011-10-18 20:15:32 +05:30
2011-02-21 20:07:38 +05:30
path = getModelName ( )
2011-10-18 20:15:32 +05:30
if ( . not . IO_open_file ( myUnit , trim ( path ) / / InputFileExtension ) ) &
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 101_pInt , ext_msg = trim ( path ) / / InputFileExtension )
2011-10-18 20:15:32 +05:30
rewind ( myUnit )
read ( myUnit , '(a1024)' ) line
2011-11-15 23:24:18 +05:30
posGeom = IO_stringPos ( line , 2_pInt )
keyword = IO_lc ( IO_StringValue ( line , posGeom , 2_pInt ) )
2011-10-18 20:15:32 +05:30
if ( keyword ( 1 : 4 ) == 'head' ) then
2011-11-15 23:24:18 +05:30
headerLength = IO_intValue ( line , posGeom , 1_pInt ) + 1_pInt
2011-10-18 20:15:32 +05:30
else
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 42_pInt )
2011-10-18 20:15:32 +05:30
endif
rewind ( myUnit )
2011-11-15 23:24:18 +05:30
do i = 1_pInt , headerLength
2011-10-18 20:15:32 +05:30
read ( myUnit , '(a1024)' ) line
2011-01-07 18:26:45 +05:30
posGeom = IO_stringPos ( line , maxNchunksGeom )
select case ( IO_lc ( IO_StringValue ( line , posGeom , 1 ) ) )
2010-06-25 17:01:05 +05:30
case ( 'dimension' )
2010-09-22 17:34:43 +05:30
gotDimension = . true .
2011-11-15 23:24:18 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posGeom , j ) ) )
2010-09-22 17:34:43 +05:30
case ( 'x' )
2011-11-15 23:24:18 +05:30
geomdimension ( 1 ) = IO_floatValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'y' )
2011-11-15 23:24:18 +05:30
geomdimension ( 2 ) = IO_floatValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'z' )
2011-11-15 23:24:18 +05:30
geomdimension ( 3 ) = IO_floatValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
end select
enddo
2010-06-25 17:01:05 +05:30
case ( 'homogenization' )
2010-09-22 17:34:43 +05:30
gotHomogenization = . true .
2011-11-15 23:24:18 +05:30
homog = IO_intValue ( line , posGeom , 2_pInt )
2010-06-25 17:01:05 +05:30
case ( 'resolution' )
2010-09-22 17:34:43 +05:30
gotResolution = . true .
2011-11-15 23:24:18 +05:30
do j = 2_pInt , 6_pInt , 2_pInt
2011-10-18 20:15:32 +05:30
select case ( IO_lc ( IO_stringValue ( line , posGeom , j ) ) )
2010-09-22 17:34:43 +05:30
case ( 'a' )
2011-11-15 23:24:18 +05:30
resolution ( 1 ) = IO_intValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'b' )
2011-11-15 23:24:18 +05:30
resolution ( 2 ) = IO_intValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
case ( 'c' )
2011-11-15 23:24:18 +05:30
resolution ( 3 ) = IO_intValue ( line , posGeom , j + 1_pInt )
2010-09-22 17:34:43 +05:30
end select
enddo
2011-10-18 20:15:32 +05:30
case ( 'picture' )
spectralPictureMode = . true .
2010-06-25 17:01:05 +05:30
end select
enddo
2011-10-18 20:15:32 +05:30
close ( myUnit )
2011-11-15 23:24:18 +05:30
if ( . not . ( gotDimension . and . gotHomogenization . and . gotResolution ) ) call IO_error ( error_ID = 45_pInt )
2011-02-14 22:51:31 +05:30
2011-08-01 15:41:32 +05:30
if ( mod ( resolution ( 1 ) , 2_pInt ) / = 0_pInt . or . &
mod ( resolution ( 2 ) , 2_pInt ) / = 0_pInt . or . &
2011-11-15 23:24:18 +05:30
( mod ( resolution ( 3 ) , 2_pInt ) / = 0_pInt . and . resolution ( 3 ) / = 1_pInt ) ) call IO_error ( error_ID = 103_pInt )
2011-08-26 19:36:37 +05:30
2011-10-24 23:56:34 +05:30
! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field
2011-11-15 23:24:18 +05:30
call CPFEM_initAll ( bc ( 1 ) % temperature , 1_pInt , 1_pInt )
2011-11-11 19:47:43 +05:30
2011-10-24 23:56:34 +05:30
!Output of geom file
2011-10-18 20:15:32 +05:30
!$OMP CRITICAL (write2out)
2011-11-04 01:02:11 +05:30
print '(a)' , ''
2011-11-15 23:24:18 +05:30
print '(a)' , '#############################################################'
2011-11-04 01:02:11 +05:30
print '(a)' , 'DAMASK spectral:'
print '(a)' , 'The spectral method boundary value problem solver for'
print '(a)' , 'the Duesseldorf Advanced Material Simulation Kit'
2011-11-15 23:24:18 +05:30
print '(a)' , '#############################################################'
2011-10-24 23:56:34 +05:30
print '(a,a)' , 'Geom File Name: ' , trim ( path ) / / '.geom'
2011-11-15 23:24:18 +05:30
print '(a)' , '============================================================='
print '(a,i12,i12,i12)' , 'resolution a b c:' , resolution
print '(a,f12.5,f12.5,f12.5)' , 'dimension x y z:' , geomdimension
2011-10-25 19:08:24 +05:30
print '(a,i5)' , 'homogenization: ' , homog
2011-10-18 20:15:32 +05:30
print '(a,L)' , 'spectralPictureMode: ' , spectralPictureMode
2011-11-15 23:24:18 +05:30
print '(a)' , '#############################################################'
2011-10-18 20:15:32 +05:30
print '(a,a)' , 'Loadcase File Name: ' , trim ( getLoadcaseName ( ) )
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( bc ( 1 ) % followFormerTrajectory ) then
2011-11-02 20:08:42 +05:30
call IO_warning ( warning_ID = 33_pInt ) ! cannot guess along trajectory for first step of first loadcase
2011-11-15 23:24:18 +05:30
bc ( 1 ) % followFormerTrajectory = . false .
2011-10-18 20:15:32 +05:30
endif
! consistency checks and output of loadcase
2011-11-11 19:47:43 +05:30
do loadcase = 1_pInt , N_Loadcases
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a)' , '============================================================='
2011-10-18 20:15:32 +05:30
print '(a,i5)' , 'Loadcase: ' , loadcase
2011-11-02 20:08:42 +05:30
write ( loadcase_string , '(i3)' ) loadcase
2011-11-15 23:24:18 +05:30
if ( . not . bc ( loadcase ) % followFormerTrajectory ) &
2011-10-18 20:15:32 +05:30
print '(a)' , 'Drop Guessing Along Trajectory'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( any ( bc ( loadcase ) % maskStress . eqv . bc ( loadcase ) % maskDeformation ) ) & ! exclusive or masking only
call IO_error ( error_ID = 31_pInt , ext_msg = loadcase_string )
if ( any ( bc ( loadcase ) % maskStress . and . transpose ( bc ( loadcase ) % maskStress ) . and . & !checking if no rotation is allowed by stress BC
2011-11-02 20:08:42 +05:30
reshape ( ( / . false . , . true . , . true . , . true . , . false . , . true . , . true . , . true . , . false . / ) , ( / 3 , 3 / ) ) ) ) &
2011-11-15 23:24:18 +05:30
call IO_error ( error_ID = 38_pInt , ext_msg = loadcase_string )
if ( bc ( loadcase ) % velGradApplied ) then
do j = 1_pInt , 3_pInt
if ( any ( bc ( loadcase ) % maskDeformation ( j , 1 : 3 ) . eqv . . true . ) . and . &
any ( bc ( loadcase ) % maskDeformation ( j , 1 : 3 ) . eqv . . false . ) ) call IO_error ( error_ID = 32_pInt , ext_msg = loadcase_string ) ! each line should be either fully or not at all defined
2011-10-18 20:15:32 +05:30
enddo
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a)' , 'Velocity Gradient:'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-10-18 20:15:32 +05:30
else
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a)' , 'Change of Deformation Gradient:'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-10-18 20:15:32 +05:30
endif
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(3(3(f12.6,x)/)\)' , merge ( math_transpose3x3 ( bc ( loadcase ) % deformation ) , &
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , ( / 3 , 3 / ) ) , transpose ( bc ( loadcase ) % maskDeformation ) )
print '(a,/,3(3(f12.6,x)/)\)' , 'Stress Boundary Condition/MPa:' , merge ( math_transpose3x3 ( bc ( loadcase ) % stress ) , &
reshape ( spread ( DAMASK_NaN , 1 , 9 ) , ( / 3 , 3 / ) ) , &
transpose ( bc ( loadcase ) % maskStress ) ) * 1e-6
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( any ( abs ( math_mul33x33 ( bc ( loadcase ) % rotation , math_transpose3x3 ( bc ( loadcase ) % rotation ) ) - math_I3 ) & ! given rotation matrix contains strain
2011-10-25 19:08:24 +05:30
> reshape ( spread ( rotation_tol , 1 , 9 ) , ( / 3 , 3 / ) ) ) &
2011-11-15 23:24:18 +05:30
. or . abs ( math_det3x3 ( bc ( loadcase ) % rotation ) ) > 1.0_pReal + rotation_tol ) call IO_error ( error_ID = 46_pInt , ext_msg = loadcase_string )
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( any ( bc ( loadcase ) % rotation / = math_I3 ) ) &
print '(a,3(3(f12.6,x)/)\)' , 'Rotation of BCs:' , math_transpose3x3 ( bc ( loadcase ) % rotation )
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % timeIncrement < 0.0_pReal ) call IO_error ( error_ID = 34_pInt , ext_msg = loadcase_string ) ! negative time increment
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a,f12.6)' , 'Temperature:' , bc ( loadcase ) % temperature
print '(a,f12.6)' , 'Time: ' , bc ( loadcase ) % timeIncrement
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % steps < 1_pInt ) call IO_error ( error_ID = 35_pInt , ext_msg = loadcase_string ) ! non-positive increment count
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a,i5)' , 'Steps: ' , bc ( loadcase ) % steps
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % outputfrequency < 1_pInt ) call IO_error ( error_ID = 36_pInt , ext_msg = loadcase_string ) ! non-positive result frequency
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a,i5)' , 'Freq. of Results Output: ' , bc ( loadcase ) % outputfrequency
2011-11-07 16:34:57 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % restartfrequency < 1_pInt ) call IO_error ( error_ID = 39_pInt , ext_msg = loadcase_string ) ! non-positive restart frequency
2011-11-07 16:34:57 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a,i5)' , 'Freq. of Restart Information Output: ' , bc ( loadcase ) % restartfrequency
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-10-18 20:15:32 +05:30
enddo
2011-07-25 22:00:21 +05:30
! Initialization of fftw (see manual on fftw.org for more details)
2011-10-24 23:56:34 +05:30
#ifdef _OPENMP
if ( DAMASK_NumThreadsInt > 0_pInt ) then
call dfftw_init_threads ( ierr )
2011-11-15 23:24:18 +05:30
if ( ierr == 0_pInt ) call IO_error ( error_ID = 104_pInt )
2011-10-24 23:56:34 +05:30
call dfftw_plan_with_nthreads ( DAMASK_NumThreadsInt )
endif
#endif
2011-07-25 22:00:21 +05:30
2011-11-15 23:24:18 +05:30
!call dfftw_timelimit(fftw_timelimit) !is not working, have to fix it in FFTW source file
2011-10-24 23:56:34 +05:30
2011-11-15 23:24:18 +05:30
select case ( IO_lc ( fftw_planner_flag ) ) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case ( 'estimate' , 'fftw_estimate' ) ! ordered from slow execution (but fast plan creation) to fast execution
2011-10-18 20:15:32 +05:30
fftw_flag = 64
case ( 'measure' , 'fftw_measure' )
fftw_flag = 0
case ( 'patient' , 'fftw_patient' )
fftw_flag = 32
case ( 'exhaustive' , 'fftw_exhaustive' )
fftw_flag = 8
2011-10-24 23:56:34 +05:30
case default
2011-11-15 23:24:18 +05:30
call IO_warning ( warning_ID = 47_pInt , ext_msg = trim ( IO_lc ( fftw_planner_flag ) ) )
2011-10-24 23:56:34 +05:30
fftw_flag = 32
2011-10-18 20:15:32 +05:30
end select
2011-11-15 23:24:18 +05:30
if ( . not . restartReadSpectral ) then ! start at first step of first loadcase
loadcase = 1_pInt
step = 1_pInt
else ! going forwarnd and use old values
do i = 1_pInt , N_Loadcases ! looping over ALL loadcases
time0 = time ! loadcase start time
timeinc = bc ( i ) % timeIncrement / bc ( i ) % steps ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
do j = 1_pInt , bc ( i ) % steps ! looping over ALL steps in current loadcase
if ( totalStepsCounter < restartReadStep ) then ! forwarding to restart step
totalStepsCounter = totalStepsCounter + 1_pInt
if ( bc ( i ) % logscale == 1_pInt ) then ! loglinear scale
if ( i == 1_pInt ) then ! 1st loadcase of loglinear scale
if ( j == 1_pInt ) then ! 1st step of 1st loadcase of loglinear scale
timeinc = bc ( 1 ) % timeIncrement * ( 2.0_pReal ** real ( 1_pInt - bc ( 1 ) % steps , pReal ) ) ! assume 1st step is equal to 2nd
else ! not-1st step of 1st loadcase of loglinear scale
timeinc = bc ( 1 ) % timeIncrement * ( 2.0_pReal ** real ( j - 1_pInt - bc ( 1 ) % steps , pReal ) )
endif
else ! not-1st loadcase of loglinear scale
timeinc = time0 * ( ( 1.0_pReal + bc ( i ) % timeIncrement / time0 ) ** real ( j / bc ( i ) % steps , pReal ) &
- ( 1.0_pReal + bc ( i ) % timeIncrement / time0 ) ** real ( ( j - 1_pInt ) / bc ( i ) % steps , pReal ) ) !ToDo: correct? how should the float casting be done
endif
endif
time = time + timeinc
endif
enddo
enddo
do i = 1_pInt , N_Loadcases ! looping over ALL loadcases
do j = 1_pInt , bc ( i ) % steps ! looping over ALL steps in current loadcase
if ( totalStepsCounter - 1_pInt < restartReadStep ) then ! forwarding to restart step
step = j
loadcase = i
endif
enddo
enddo
2011-02-07 20:05:42 +05:30
2011-11-15 23:24:18 +05:30
endif
print * , totalStepsCounter
print * , loadcase
print * , step
pause
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
2011-04-06 15:28:17 +05:30
! Loop over loadcases defined in the loadcase file
2011-11-11 19:47:43 +05:30
do loadcase = loadcase , N_Loadcases
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
2011-07-13 22:03:12 +05:30
time0 = time ! loadcase start time
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % followFormerTrajectory ) then ! continue to guess along former trajectory where applicable
2011-07-13 22:03:12 +05:30
guessmode = 1.0_pReal
else
guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step
endif
2011-07-06 18:40:18 +05:30
2011-11-15 23:24:18 +05:30
mask_defgrad = merge ( ones , zeroes , bc ( loadcase ) % maskDeformation )
mask_stress = merge ( ones , zeroes , bc ( loadcase ) % maskStress )
size_reduced = count ( bc ( loadcase ) % maskStressVector )
2011-08-26 19:36:37 +05:30
allocate ( c_reduced ( size_reduced , size_reduced ) ) ; c_reduced = 0.0_pReal
allocate ( s_reduced ( size_reduced , size_reduced ) ) ; s_reduced = 0.0_pReal
2011-08-10 21:32:13 +05:30
2011-11-15 23:24:18 +05:30
timeinc = bc ( loadcase ) % timeIncrement / bc ( loadcase ) % steps ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
fDot = bc ( loadcase ) % deformation ! only valid for given fDot. will be overwritten later in case L is given
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
! loop oper steps defined in input file for current loadcase
2011-11-15 23:24:18 +05:30
do step = step , bc ( loadcase ) % steps
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
2011-11-15 23:24:18 +05:30
!*************************************************************
! Initialization Start
!*************************************************************
if ( stepZero == 1_pInt ) then ! we start
stepZero = 0_pInt
if ( regrid == . true . ) then ! 'real' start vs. regrid start
call dfftw_destroy_plan ( fftw_plan ( 1 ) ) ; call dfftw_destroy_plan ( fftw_plan ( 2 ) )
if ( debugDivergence ) call dfftw_destroy_plan ( fftw_plan ( 3 ) )
deallocate ( defgrad )
deallocate ( defgradold )
deallocate ( coordinates )
deallocate ( temperature )
deallocate ( xi )
deallocate ( workfft )
! here we have to create the new geometry and assign the values from the previous step
endif
allocate ( defgrad ( resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) , 3 , 3 ) ) ; defgrad = 0.0_pReal
allocate ( defgradold ( resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) , 3 , 3 ) ) ; defgradold = 0.0_pReal
allocate ( coordinates ( 3 , resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) ) ) ; coordinates = 0.0_pReal
allocate ( temperature ( resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) ) ) ; temperature = bc ( 1 ) % temperature ! start out isothermally
allocate ( xi ( 3 , resolution ( 1 ) / 2 + 1 , resolution ( 2 ) , resolution ( 3 ) ) ) ; xi = 0.0_pReal
allocate ( workfft ( resolution ( 1 ) + 2 , resolution ( 2 ) , resolution ( 3 ) , 3 , 3 ) ) ; workfft = 0.0_pReal
if ( debugDivergence ) allocate ( divergence ( resolution ( 1 ) + 2 , resolution ( 2 ) , resolution ( 3 ) , 3 ) ) ; divergence = 0.0_pReal
wgt = 1.0_pReal / real ( resolution ( 1 ) * resolution ( 2 ) * resolution ( 3 ) , pReal )
call dfftw_plan_many_dft_r2c ( fftw_plan ( 1 ) , 3 , ( / resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) / ) , 9 , &
workfft , ( / resolution ( 1 ) + 2_pInt , resolution ( 2 ) , resolution ( 3 ) / ) , 1 , ( resolution ( 1 ) + 2_pInt ) * resolution ( 2 ) * resolution ( 3 ) , &
workfft , ( / resolution ( 1 ) / 2_pInt + 1_pInt , resolution ( 2 ) , resolution ( 3 ) / ) , 1 , ( resolution ( 1 ) / 2_pInt + 1_pInt ) * resolution ( 2 ) * resolution ( 3 ) , fftw_flag )
call dfftw_plan_many_dft_c2r ( fftw_plan ( 2 ) , 3 , ( / resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) / ) , 9 , &
workfft , ( / resolution ( 1 ) / 2_pInt + 1_pInt , resolution ( 2 ) , resolution ( 3 ) / ) , 1 , ( resolution ( 1 ) / 2_pInt + 1_pInt ) * resolution ( 2 ) * resolution ( 3 ) , &
workfft , ( / resolution ( 1 ) + 2_pInt , resolution ( 2 ) , resolution ( 3 ) / ) , 1 , ( resolution ( 1 ) + 2_pInt ) * resolution ( 2 ) * resolution ( 3 ) , fftw_flag )
if ( debugDivergence ) &
call dfftw_plan_many_dft_c2r ( fftw_plan ( 3 ) , 3 , ( / resolution ( 1 ) , resolution ( 2 ) , resolution ( 3 ) / ) , 3 , &
divergence , ( / resolution ( 1 ) / 2_pInt + 1_pInt , resolution ( 2 ) , resolution ( 3 ) / ) , 1 , ( resolution ( 1 ) / 2_pInt + 1_pInt ) * resolution ( 2 ) * resolution ( 3 ) , &
divergence , ( / resolution ( 1 ) + 2_pInt , resolution ( 2 ) , resolution ( 3 ) / ) , 1 , ( resolution ( 1 ) + 2_pInt ) * resolution ( 2 ) * resolution ( 3 ) , fftw_flag )
if ( debugGeneral ) then
!$OMP CRITICAL (write2out)
write ( 6 , * ) 'FFTW initialized'
!$OMP END CRITICAL (write2out)
endif
if ( . not . restartReadSpectral ) then ! no deformation at the beginning
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = math_I3
enddo ; enddo ; enddo
else ! using old values
if ( IO_read_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , trim ( getModelName ( ) ) , size ( defgrad ) ) ) then
read ( 777 , rec = 1 ) defgrad
2011-11-11 19:47:43 +05:30
close ( 777 )
endif
2011-11-15 23:24:18 +05:30
defgradold = defgrad
defgradAim = 0.0_pReal
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
defgradAim = defgradAim + defgrad ( i , j , k , 1 : 3 , 1 : 3 ) ! calculating old average deformation
enddo ; enddo ; enddo
defgradAim = defgradAim * wgt
defgradAimOld = defgradAim
2011-11-07 16:34:57 +05:30
endif
2011-11-15 23:24:18 +05:30
ielem = 0_pInt
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
ielem = ielem + 1_pInt
coordinates ( 1 : 3 , i , j , k ) = mesh_ipCenterOfGravity ( 1 : 3 , 1 , ielem ) ! set to initial coordinates ToDo: SHOULD BE UPDATED TO CURRENT POSITION IN FUTURE REVISIONS!!! But do we know them? I don't think so. Otherwise we don't need geometry reconstruction
call CPFEM_general ( 2_pInt , coordinates ( 1 : 3 , i , j , k ) , math_I3 , math_I3 , temperature ( i , j , k ) , 0.0_pReal , ielem , 1_pInt , cstress , dsde , pstress , dPdF )
c_current = c_current + dPdF
enddo ; enddo ; enddo
c0_reference = c_current * wgt ! linear reference material stiffness
c_prev = math_rotate_forward3x3x3x3 ( c0_reference , bc ( loadcase ) % rotation ) ! rotate_forward: lab -> load system
if ( debugGeneral ) then
!$OMP CRITICAL (write2out)
write ( 6 , * ) 'First Call to CPFEM_general finished'
!$OMP END CRITICAL (write2out)
endif
do k = 1_pInt , resolution ( 3 ) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
k_s ( 3 ) = k - 1_pInt
if ( k > resolution ( 3 ) / 2_pInt + 1_pInt ) k_s ( 3 ) = k_s ( 3 ) - resolution ( 3 )
do j = 1_pInt , resolution ( 2 )
k_s ( 2 ) = j - 1_pInt
if ( j > resolution ( 2 ) / 2_pInt + 1_pInt ) k_s ( 2 ) = k_s ( 2 ) - resolution ( 2 )
do i = 1 , resolution ( 1 ) / 2_pInt + 1_pInt
k_s ( 1 ) = i - 1_pInt
xi ( 3 , i , j , k ) = 0.0_pReal ! 2D case
if ( resolution ( 3 ) > 1_pInt ) xi ( 3 , i , j , k ) = real ( k_s ( 3 ) , pReal ) / geomdimension ( 3 ) ! 3D case
xi ( 2 , i , j , k ) = real ( k_s ( 2 ) , pReal ) / geomdimension ( 2 )
xi ( 1 , i , j , k ) = real ( k_s ( 1 ) , pReal ) / geomdimension ( 1 )
enddo ; enddo ; enddo
! remove highest frequencies for calculation of divergence (CAREFULL, they will be used for pre calculatet gamma operator!)
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 ) / 2_pInt + 1_pInt
if ( k == resolution ( 3 ) / 2_pInt + 1_pInt ) xi ( 3 , i , j , k ) = 0.0_pReal
if ( j == resolution ( 2 ) / 2_pInt + 1_pInt ) xi ( 2 , i , j , k ) = 0.0_pReal
if ( i == resolution ( 1 ) / 2_pInt + 1_pInt ) xi ( 1 , i , j , k ) = 0.0_pReal
enddo ; enddo ; enddo
if ( memory_efficient ) then ! allocate just single fourth order tensor
allocate ( gamma_hat ( 1 , 1 , 1 , 3 , 3 , 3 , 3 ) ) ; gamma_hat = 0.0_pReal
else ! precalculation of gamma_hat field
allocate ( gamma_hat ( resolution ( 1 ) / 2_pInt + 1_pInt , resolution ( 2 ) , resolution ( 3 ) , 3 , 3 , 3 , 3 ) ) ; gamma_hat = 0.0_pReal
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 ) / 2_pInt + 1_pInt
if ( any ( xi ( : , i , j , k ) / = 0.0_pReal ) ) then
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
xiDyad ( l , m ) = xi ( l , i , j , k ) * xi ( m , i , j , k )
enddo ; enddo
temp33_Real = math_inv3x3 ( math_mul3333xx33 ( c0_reference , xiDyad ) )
else
xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt ; do p = 1_pInt , 3_pInt
gamma_hat ( i , j , k , l , m , n , p ) = - 0.25 * ( temp33_Real ( l , n ) + temp33_Real ( n , l ) ) * &
( xiDyad ( m , p ) + xiDyad ( p , m ) )
enddo ; enddo ; enddo ; enddo
enddo ; enddo ; enddo
endif
! write header of output file
!$OMP CRITICAL (write2out)
open ( 538 , file = trim ( getSolverWorkingDirectoryName ( ) ) / / trim ( getSolverJobName ( ) ) &
/ / '.spectralOut' , form = 'UNFORMATTED' , status = 'REPLACE' )
write ( 538 ) , 'load' , trim ( getLoadcaseName ( ) )
write ( 538 ) , 'workingdir' , trim ( getSolverWorkingDirectoryName ( ) )
write ( 538 ) , 'geometry' , trim ( getSolverJobName ( ) ) / / InputFileExtension
write ( 538 ) , 'resolution' , resolution
write ( 538 ) , 'dimension' , geomdimension
write ( 538 ) , 'materialpoint_sizeResults' , materialpoint_sizeResults
write ( 538 ) , 'loadcases' , N_Loadcases
write ( 538 ) , 'logscale' , bc ( loadcase ) % logscale ! one entry per loadcase (0: linear, 1: log)
write ( 538 ) , 'frequencies' , bc ( loadcase ) % outputfrequency ! one entry per loadcase
write ( 538 ) , 'times' , bc ( loadcase ) % timeIncrement ! one entry per loadcase
bc ( 1 ) % steps = bc ( 1 ) % steps + 1_pInt
write ( 538 ) , 'increments' , bc ( loadcase ) % steps ! one entry per loadcase ToDo: rename keyword to steps
bc ( 1 ) % steps = bc ( 1 ) % steps - 1_pInt
write ( 538 ) , 'startingIncrement' , totalStepsCounter
write ( 538 ) , 'eoh' ! end of header
write ( 538 ) , materialpoint_results ( : , 1 , : ) ! initial (non-deformed) results !ToDo: define array size
!$OMP END CRITICAL (write2out)
endif
!*************************************************************
! Initialization End
!*************************************************************
totalStepsCounter = totalStepsCounter + 1_pInt
if ( mod ( step - 1_pInt , bc ( loadcase ) % restartFrequency ) == 0_pInt ) then ! at frequency of writing restart information
restartWrite = . true . ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?)
if ( IO_write_jobBinaryFile ( 777 , 'convergedSpectralDefgrad' , size ( defgrad ) ) ) then ! and writing deformation gradient field to file
write ( 777 , rec = 1 ) defgrad
close ( 777 )
endif
else
restartWrite = . false .
endif
if ( bc ( loadcase ) % logscale == 1_pInt ) then ! loglinear scale
if ( loadcase == 1_pInt ) then ! 1st loadcase of loglinear scale
if ( step == 1_pInt ) then ! 1st step of 1st loadcase of loglinear scale
timeinc = bc ( 1 ) % timeIncrement * ( 2.0_pReal ** real ( 1_pInt - bc ( 1 ) % steps , pReal ) ) ! assume 1st step is equal to 2nd
else ! not-1st step of 1st loadcase of loglinear scale
timeinc = bc ( 1 ) % timeIncrement * ( 2.0_pReal ** real ( step - 1_pInt - bc ( 1 ) % steps , pReal ) )
endif
else ! not-1st loadcase of loglinear scale
timeinc = time0 * ( ( 1.0_pReal + bc ( loadcase ) % timeIncrement / time0 ) ** real ( step / bc ( loadcase ) % steps , pReal ) &
- ( 1.0_pReal + bc ( loadcase ) % timeIncrement / time0 ) ** real ( ( step - 1_pInt ) / bc ( loadcase ) % steps , pReal ) )
2011-07-06 18:40:18 +05:30
endif
2011-11-15 23:24:18 +05:30
endif
time = time + timeinc
2011-07-14 15:07:31 +05:30
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % velGradApplied ) & ! calculate fDot from given L and current F
fDot = math_mul33x33 ( bc ( loadcase ) % deformation , defgradAim )
2011-07-14 15:07:31 +05:30
2011-10-24 23:56:34 +05:30
!winding forward of deformation aim in loadcase system
2011-07-14 15:07:31 +05:30
temp33_Real = defgradAim
defgradAim = defgradAim &
+ guessmode * mask_stress * ( defgradAim - defgradAimOld ) &
2011-08-31 20:07:01 +05:30
+ mask_defgrad * fDot * timeinc
2011-07-13 22:03:12 +05:30
defgradAimOld = temp33_Real
2011-07-14 15:07:31 +05:30
2011-07-13 22:03:12 +05:30
! update local deformation gradient
2011-11-15 23:24:18 +05:30
if ( any ( bc ( loadcase ) % rotation / = math_I3 ) ) then ! lab and loadcase coordinate system are NOT the same
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
2011-10-25 19:08:24 +05:30
temp33_Real = defgrad ( i , j , k , 1 : 3 , 1 : 3 )
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % velGradApplied ) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33 ( bc ( loadcase ) % deformation , &
math_rotate_forward3x3 ( defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , bc ( loadcase ) % rotation ) )
2011-10-25 19:08:24 +05:30
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = defgrad ( i , j , k , 1 : 3 , 1 : 3 ) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) - defgradold ( i , j , k , 1 : 3 , 1 : 3 ) ) & ! guessing...
+ math_rotate_backward3x3 ( ( 1.0_pReal - guessmode ) * mask_defgrad * fDot , &
2011-11-15 23:24:18 +05:30
bc ( loadcase ) % rotation ) * timeinc ! apply the prescribed value where deformation is given if not guessing
2011-10-25 19:08:24 +05:30
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Real
enddo ; enddo ; enddo
2011-11-15 23:24:18 +05:30
else ! one coordinate system for lab and loadcase, save some multiplications
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
2011-10-25 19:08:24 +05:30
temp33_Real = defgrad ( i , j , k , 1 : 3 , 1 : 3 )
2011-11-15 23:24:18 +05:30
if ( bc ( loadcase ) % velGradApplied ) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33 ( bc ( loadcase ) % deformation , defgradold ( i , j , k , 1 : 3 , 1 : 3 ) )
2011-10-25 19:08:24 +05:30
defgrad ( i , j , k , 1 : 3 , 1 : 3 ) = defgrad ( i , j , k , 1 : 3 , 1 : 3 ) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) - defgradold ( i , j , k , 1 : 3 , 1 : 3 ) ) & ! guessing...
+ ( 1.0_pReal - guessmode ) * mask_defgrad * fDot * timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) = temp33_Real
enddo ; enddo ; enddo
endif
2011-08-26 19:36:37 +05:30
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
2010-07-02 22:45:53 +05:30
2011-08-26 19:36:37 +05:30
CPFEM_mode = 1_pInt ! winding forward
2010-07-02 22:45:53 +05:30
iter = 0_pInt
2011-08-26 19:36:37 +05:30
err_div = 2.0_pReal * err_div_tol ! go into loop
if ( size_reduced > 0_pInt ) then ! calculate compliance in case stress BC is applied
c_prev99 = math_Plain3333to99 ( c_prev )
k = 0_pInt ! build reduced stiffness
2011-11-15 23:24:18 +05:30
do n = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( n ) ) then
2011-08-26 19:36:37 +05:30
k = k + 1_pInt
j = 0_pInt
2011-11-15 23:24:18 +05:30
do m = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( m ) ) then
2011-08-26 19:36:37 +05:30
j = j + 1_pInt
c_reduced ( k , j ) = c_prev99 ( n , m )
endif ; enddo ; endif ; enddo
call math_invert ( size_reduced , c_reduced , s_reduced , i , errmatinv ) ! invert reduced stiffness
2011-11-02 20:08:42 +05:30
if ( errmatinv ) call IO_error ( error_ID = 800 )
2011-08-26 19:36:37 +05:30
s_prev99 = 0.0_pReal ! build full compliance
k = 0_pInt
2011-11-15 23:24:18 +05:30
do n = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( n ) ) then
2011-08-26 19:36:37 +05:30
k = k + 1_pInt
j = 0_pInt
2011-11-15 23:24:18 +05:30
do m = 1_pInt , 9_pInt
if ( bc ( loadcase ) % maskStressVector ( m ) ) then
2011-08-26 19:36:37 +05:30
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
!$OMP CRITICAL (write2out)
print '(a)' , '#############################################################'
print '(A,I5.5,A,es12.6)' , 'Increment ' , totalStepsCounter , ' Time ' , time
if ( restartWrite . eq . . true . ) print '(A)' , 'Writing converged Results of previous Step for Restart'
!$OMP END CRITICAL (write2out)
added fftw3 as fft(library will not versioned, should be in a linkable folder) , did some corrections on the code, splitted main file up (allows use of makefile), added makefile
changes on mpie_spectral.f90:
new structure, changed variable names, now using defgrad instead of disgrad, cleaned up, removed augmented Lagrange.
ToDo: Implement Augmented Lagrange again (but then a working version), implement Large strain, think about complex-to real-transform backwards, try to implement MP-support
2010-08-27 22:09:38 +05:30
!*************************************************************
2010-10-13 21:34:44 +05:30
! convergence loop
2011-04-06 15:28:17 +05:30
do while ( iter < itmax . and . &
2011-06-06 20:50:28 +05:30
( err_div > err_div_tol . or . &
2011-08-26 19:36:37 +05:30
err_stress > err_stress_tol ) )
2011-02-07 20:05:42 +05:30
iter = iter + 1_pInt
2010-10-13 21:34:44 +05:30
!*************************************************************
2011-11-15 23:24:18 +05:30
print '(a)' , '============================================================='
print '(5(A,I5.5))' , 'Loadcase ' , loadcase , ' Step ' , step , '/' , bc ( loadcase ) % steps , '@Iteration ' , iter , '/' , itmax
do n = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
defgrad_av ( m , n ) = sum ( defgrad ( 1 : resolution ( 1 ) , 1 : resolution ( 2 ) , 1 : resolution ( 3 ) , m , n ) ) * wgt
2011-08-31 20:07:01 +05:30
enddo ; enddo
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a,/,3(3(f12.7,x)/)\)' , 'Deformation Gradient:' , math_transpose3x3 ( defgrad_av )
print '(A)' , '... Update Stress Field (Constitutive Evaluation P(F)) ......'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-08-26 19:36:37 +05:30
ielem = 0_pInt
2011-11-15 23:24:18 +05:30
defgradDetMax = - 99 9.0_pReal
defgradDetMin = 99 9.0_pReal
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
2011-11-11 19:47:43 +05:30
defgradDet = math_det3x3 ( defgrad ( i , j , k , 1 : 3 , 1 : 3 ) )
defgradDetMax = max ( defgradDetMax , defgradDet )
defgradDetMin = min ( defgradDetMin , defgradDet )
2011-11-15 23:24:18 +05:30
ielem = ielem + 1_pInt
call CPFEM_general ( 3_pInt , & ! collect cycle
coordinates ( 1 : 3 , i , j , k ) , defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) , &
2011-08-26 19:36:37 +05:30
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , &
cstress , dsde , pstress , dPdF )
enddo ; enddo ; enddo
2011-11-15 23:24:18 +05:30
print '(a,x,es10.4)' , 'Maximum Determinant of Deformation:' , defgradDetMax
print '(a,x,es10.4)' , 'Minimum Determinant of Deformation:' , defgradDetMin
workfft = 0.0_pReal ! needed because of the padding for FFTW
2011-08-26 19:36:37 +05:30
c_current = 0.0_pReal
ielem = 0_pInt
2011-11-15 23:24:18 +05:30
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 )
2011-08-26 19:36:37 +05:30
ielem = ielem + 1_pInt
call CPFEM_general ( CPFEM_mode , & ! first element in first iteration retains CPFEM_mode 1,
coordinates ( 1 : 3 , i , j , k ) , &
2011-11-15 23:24:18 +05:30
defgradold ( i , j , k , 1 : 3 , 1 : 3 ) , defgrad ( i , j , k , 1 : 3 , 1 : 3 ) , & ! others get 2 (saves winding forward effort)
2011-08-26 19:36:37 +05:30
temperature ( i , j , k ) , timeinc , ielem , 1_pInt , &
cstress , dsde , pstress , dPdF )
CPFEM_mode = 2_pInt
2011-11-15 23:24:18 +05:30
workfft ( i , j , k , 1 : 3 , 1 : 3 ) = pstress ! build up average P-K stress
2011-08-26 19:36:37 +05:30
c_current = c_current + dPdF
enddo ; enddo ; enddo
2011-11-15 23:24:18 +05:30
restartWrite = . false . ! ToDo: don't know if we need it. Depends on how CPFEM_general is writing results
do n = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
pstress_av ( m , n ) = sum ( workfft ( 1 : resolution ( 1 ) , 1 : resolution ( 2 ) , 1 : resolution ( 3 ) , m , n ) ) * wgt
2011-08-26 19:36:37 +05:30
enddo ; enddo
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a,/,3(3(f12.7,x)/)\)' , 'Piola-Kirchhoff Stress / MPa: ' , math_transpose3x3 ( pstress_av ) / 1.e6
2011-11-02 20:08:42 +05:30
2011-08-26 19:36:37 +05:30
err_stress_tol = 0.0_pReal
2011-11-15 23:24:18 +05:30
pstress_av_load = math_rotate_forward3x3 ( pstress_av , bc ( loadcase ) % rotation )
if ( size_reduced > 0_pInt ) then ! calculate stress BC if applied
err_stress = maxval ( abs ( mask_stress * ( pstress_av_load - bc ( loadcase ) % stress ) ) ) ! maximum deviaton (tensor norm not applicable)
err_stress_tol = maxval ( abs ( mask_defgrad * pstress_av_load ) ) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
print '(A)' , '... Correcting Deformation Gradient to Fullfill BCs .........'
print '(2(a,es10.4))' , 'Error Stress = ' , err_stress , ', Tol. = ' , err_stress_tol
defgradAimCorr = - math_mul3333xx33 ( s_prev , ( ( pstress_av_load - bc ( loadcase ) % stress ) ) ) ! residual on given stress components
2011-08-26 19:36:37 +05:30
defgradAim = defgradAim + defgradAimCorr
2011-11-15 23:24:18 +05:30
print '(a,/,3(3(f12.7,x)/)\)' , 'New Deformation Aim: ' , math_transpose3x3 ( math_rotate_backward3x3 ( &
defgradAim , bc ( loadcase ) % rotation ) )
print '(a,x,es10.4)' , 'Determinant of New Deformation Aim:' , math_det3x3 ( defgradAim )
2011-08-26 19:36:37 +05:30
endif
2011-11-15 23:24:18 +05:30
print '(A)' , '... Calculating Equilibrium Using Spectral Method ...........'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
call dfftw_execute_dft_r2c ( fftw_plan ( 1 ) , workfft , workfft ) ! FFT of pstress
2011-08-31 20:07:01 +05:30
p_hat_avg = sqrt ( maxval ( math_eigenvalues3x3 ( math_mul33x33 ( workfft ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) , & ! L_2 norm of average stress in fourier space,
math_transpose3x3 ( workfft ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) ) ) ) ) ) ! ignore imaginary part as it is always zero for real only input))
err_div = 0.0_pReal
2011-11-15 23:24:18 +05:30
err_div_max = 0.0_pReal
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 ) / 2_pInt + 1_pInt
err_div = err_div + sqrt ( sum ( ( & ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain)
math_mul33x3_complex ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) + &
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img , &
xi ( 1 : 3 , i , j , k ) ) &
) ** 2.0_pReal ) )
if ( debugDivergence ) &
err_div_max = max ( err_div_max , abs ( sqrt ( sum ( ( & ! maximum of L two norm of div(stress) in fourier space (Suquet large strain)
math_mul33x3_complex ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) + &
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img , &
xi ( 1 : 3 , i , j , k ) ) &
) ** 2.0_pReal ) ) ) )
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
2011-11-15 23:24:18 +05:30
correctionFactor = minval ( geomdimension ) * wgt ** ( - 1.0_pReal / 4.0_pReal ) ! multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency
if ( resolution ( 3 ) == 1_pInt ) correctionFactor = minval ( geomdimension ( 1 : 2 ) ) * wgt ** ( - 1.0_pReal / 4.0_pReal ) ! 2D case, ToDo: correct?
if ( . not . divergence_correction ) correctionFactor = 1.0_pReal
err_div = err_div * wgt / p_hat_avg * correctionFactor ! weighting by points and average stress and multiplying with correction factor
err_div_max = err_div_max / p_hat_avg * correctionFactor ! weighting by average stress and multiplying with correction factor
2011-02-09 23:17:28 +05:30
2011-08-26 19:36:37 +05:30
if ( memory_efficient ) then ! memory saving version, on-the-fly calculation of gamma_hat
2011-11-15 23:24:18 +05:30
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 ) / 2_pInt + 1_pInt
2011-07-07 15:33:55 +05:30
if ( any ( xi ( : , i , j , k ) / = 0.0_pReal ) ) then
2011-11-15 23:24:18 +05:30
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt
2011-07-07 20:57:35 +05:30
xiDyad ( l , m ) = xi ( l , i , j , k ) * xi ( m , i , j , k )
2011-07-07 15:33:55 +05:30
enddo ; enddo
2011-08-10 21:32:13 +05:30
temp33_Real = math_inv3x3 ( math_mul3333xx33 ( c0_reference , xiDyad ) )
2011-07-07 15:33:55 +05:30
else
2011-07-07 20:57:35 +05:30
xiDyad = 0.0_pReal
2011-07-07 15:33:55 +05:30
temp33_Real = 0.0_pReal
endif
2011-11-15 23:24:18 +05:30
do l = 1_pInt , 3_pInt ; do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt ; do p = 1_pInt , 3_pInt
2011-07-07 15:33:55 +05:30
gamma_hat ( 1 , 1 , 1 , l , m , n , p ) = - 0.25_pReal * ( temp33_Real ( l , n ) + temp33_Real ( n , l ) ) * &
2011-07-07 20:57:35 +05:30
( xiDyad ( m , p ) + xiDyad ( p , m ) )
2011-07-07 15:33:55 +05:30
enddo ; enddo ; enddo ; enddo
2011-11-15 23:24:18 +05:30
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
temp33_Complex ( m , n ) = sum ( gamma_hat ( 1 , 1 , 1 , m , n , 1 : 3 , 1 : 3 ) * ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) &
+ workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img ) )
2011-07-07 15:33:55 +05:30
enddo ; enddo
2011-11-15 23:24:18 +05:30
workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) = real ( temp33_Complex )
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) = aimag ( temp33_Complex )
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
else ! use precalculated gamma-operator
2011-11-15 23:24:18 +05:30
do k = 1_pInt , resolution ( 3 ) ; do j = 1_pInt , resolution ( 2 ) ; do i = 1_pInt , resolution ( 1 ) / 2_pInt + 1_pInt
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt
temp33_Complex ( m , n ) = sum ( gamma_hat ( i , j , k , m , n , 1 : 3 , 1 : 3 ) * ( workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) &
+ workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) * img ) )
2011-02-21 20:07:38 +05:30
enddo ; enddo
2011-11-15 23:24:18 +05:30
workfft ( i * 2_pInt - 1_pInt , j , k , 1 : 3 , 1 : 3 ) = real ( temp33_Complex )
workfft ( i * 2_pInt , j , k , 1 : 3 , 1 : 3 ) = aimag ( temp33_Complex )
2011-08-26 19:36:37 +05:30
enddo ; enddo ; enddo
endif
2011-11-15 23:24:18 +05:30
if ( debugDivergence ) then
divergence = 0.0
do k = 1 , resolution ( 3 ) ; do j = 1 , resolution ( 2 ) ; do i = 1 , resolution ( 1 ) / 2 + 1
divergence ( i , j , k , 1 ) = ( workfft ( i * 2 - 1 , j , k , 1 , 1 ) + workfft ( i * 2 , j , k , 1 , 1 ) * img ) * xi ( 1 , i , j , k ) * img * pi * 2.0 &
+ ( workfft ( i * 2 - 1 , j , k , 2 , 1 ) + workfft ( i * 2 , j , k , 2 , 1 ) * img ) * xi ( 2 , i , j , k ) * img * pi * 2.0 &
+ ( workfft ( i * 2 - 1 , j , k , 3 , 1 ) + workfft ( i * 2 , j , k , 3 , 1 ) * img ) * xi ( 3 , i , j , k ) * img * pi * 2.0
divergence ( i , j , k , 2 ) = ( workfft ( i * 2 - 1 , j , k , 1 , 2 ) + workfft ( i * 2 , j , k , 1 , 2 ) * img ) * xi ( 1 , i , j , k ) * img * pi * 2.0 &
+ ( workfft ( i * 2 - 1 , j , k , 2 , 2 ) + workfft ( i * 2 , j , k , 2 , 2 ) * img ) * xi ( 2 , i , j , k ) * img * pi * 2.0 &
+ ( workfft ( i * 2 - 1 , j , k , 3 , 2 ) + workfft ( i * 2 , j , k , 3 , 2 ) * img ) * xi ( 3 , i , j , k ) * img * pi * 2.0
divergence ( i , j , k , 3 ) = ( workfft ( i * 2 - 1 , j , k , 1 , 3 ) + workfft ( i * 2 , j , k , 1 , 3 ) * img ) * xi ( 1 , i , j , k ) * img * pi * 2.0 &
+ ( workfft ( i * 2 - 1 , j , k , 2 , 3 ) + workfft ( i * 2 , j , k , 2 , 3 ) * img ) * xi ( 2 , i , j , k ) * img * pi * 2.0 &
+ ( workfft ( i * 2 - 1 , j , k , 3 , 3 ) + workfft ( i * 2 , j , k , 3 , 3 ) * img ) * xi ( 3 , i , j , k ) * img * pi * 2.0
enddo ; enddo ; enddo
call dfftw_execute_dft_c2r ( fftw_plan ( 3 ) , divergence , divergence )
endif
2011-08-26 19:36:37 +05:30
! average strain
2011-11-15 23:24:18 +05:30
workfft ( 1 , 1 , 1 , 1 : 3 , 1 : 3 ) = defgrad_av - math_I3 ! zero frequency (real part)
workfft ( 2 , 1 , 1 , 1 : 3 , 1 : 3 ) = 0.0_pReal ! zero frequency (imaginary part)
2011-08-26 19:36:37 +05:30
2011-10-18 20:15:32 +05:30
call dfftw_execute_dft_c2r ( fftw_plan ( 2 ) , workfft , workfft )
2011-08-26 19:36:37 +05:30
defgrad = defgrad + workfft ( 1 : resolution ( 1 ) , : , : , : , : ) * wgt
2011-10-24 23:56:34 +05:30
do m = 1 , 3 ; do n = 1 , 3
2011-08-26 19:36:37 +05:30
defgrad_av ( m , n ) = sum ( defgrad ( : , : , : , m , n ) ) * wgt
2011-10-24 23:56:34 +05:30
enddo ; enddo
2011-11-15 23:24:18 +05:30
defgradAim_lab = math_rotate_backward3x3 ( defgradAim , bc ( loadcase ) % rotation )
2011-10-24 23:56:34 +05:30
do m = 1 , 3 ; do n = 1 , 3
defgrad ( : , : , : , m , n ) = defgrad ( : , : , : , m , n ) + ( defgradAim_lab ( m , n ) - defgrad_av ( m , n ) ) ! anticipated target minus current state
2011-08-26 19:36:37 +05:30
enddo ; enddo
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(2(a,es10.4))' , 'Error Divergence = ' , err_div , ', Tol. = ' , err_div_tol
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-08-26 19:36:37 +05:30
2010-10-13 21:34:44 +05:30
enddo ! end looping when convergency is achieved
2010-10-27 22:45:49 +05:30
2011-11-15 23:24:18 +05:30
c_prev = math_rotate_forward3x3x3x3 ( c_current * wgt , bc ( loadcase ) % rotation ) ! calculate stiffness for next step
2011-10-24 23:56:34 +05:30
!ToDo: Incfluence for next loadcase
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a)' , '============================================================='
2011-10-18 20:15:32 +05:30
if ( err_div < = err_div_tol . and . err_stress < = err_stress_tol ) then
2011-11-15 23:24:18 +05:30
print '(A,I5.5,A)' , 'Increment ' , totalStepsCounter , ' Converged'
2011-08-31 20:07:01 +05:30
else
2011-11-15 23:24:18 +05:30
print '(A,I5.5,A)' , 'Increment ' , totalStepsCounter , ' NOT Converged'
2011-10-18 20:15:32 +05:30
notConvergedCounter = notConvergedCounter + 1
2011-08-31 20:07:01 +05:30
endif
2011-11-15 23:24:18 +05:30
if ( mod ( totalStepsCounter - 1_pInt , bc ( loadcase ) % outputfrequency ) == 0_pInt ) then ! at output frequency
print '(A)' , '... Writing Results to File .................................'
write ( 538 ) , materialpoint_results ( : , 1 , : ) ! write result to file
endif
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2010-09-06 15:30:59 +05:30
enddo ! end looping over steps in current loadcase
2011-08-26 19:36:37 +05:30
deallocate ( c_reduced )
deallocate ( s_reduced )
2011-11-15 23:24:18 +05:30
step = 1_pInt ! Reset Step Counter
2011-08-26 19:36:37 +05:30
enddo ! end looping over loadcases
2011-11-02 20:08:42 +05:30
!$OMP CRITICAL (write2out)
2011-11-15 23:24:18 +05:30
print '(a)' , '#############################################################'
print '(a,i5.5,a,i5.5,a)' , 'Of ' , totalStepsCounter - restartReadStep , ' Calculated Steps, ' , notConvergedCounter , ' Steps did not Converge!'
2011-11-02 20:08:42 +05:30
!$OMP END CRITICAL (write2out)
2011-01-07 18:26:45 +05:30
close ( 538 )
2011-10-18 20:15:32 +05:30
call dfftw_destroy_plan ( fftw_plan ( 1 ) ) ; call dfftw_destroy_plan ( fftw_plan ( 2 ) )
2011-11-15 23:24:18 +05:30
if ( debugDivergence ) call dfftw_destroy_plan ( fftw_plan ( 3 ) )
2011-05-11 22:08:45 +05:30
end program DAMASK_spectral
2010-06-10 20:21:10 +05:30
!********************************************************************
! quit subroutine to satisfy IO_error
!
!********************************************************************
2010-06-08 15:38:15 +05:30
subroutine quit ( id )
use prec
implicit none
integer ( pInt ) id
stop
2011-05-24 21:27:59 +05:30
end subroutine