From 5ef73e164af70cc0faaf4ef25f58912af61fa5eb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Nov 2011 17:54:18 +0000 Subject: [PATCH] restructured algorithm: moved into loop to reallocate fields and replan FFTW in case resolution changes during runtime ==> regridding introduced parameters for selective debugging of spectral code and partly introduced the advanced divergence calculation again which is controlled by debug.config added switch in numerics to control divergence behavior (uncorrected and corrected by phenomenological factor) added precision directive to all values I found --- code/DAMASK_spectral.f90 | 950 +++++++++++++++++++++------------------ code/IO.f90 | 2 + code/config/debug.config | 4 + code/debug.f90 | 14 +- code/numerics.f90 | 108 +++-- 5 files changed, 583 insertions(+), 495 deletions(-) diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index bd010a9a6..0c39f7cd3 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -39,7 +39,8 @@ ! -l (--load, --loadcase) PathToLoadFile/NameOfLoadFile.load ! - PathToGeomFile will be the working directory ! - make sure the file "material.config" exists in the working -! directory. For further configuration use "numerics.config" +! directory. For further configuration use "numerics.config" and +! "numerics.config" !******************************************************************** program DAMASK_spectral !******************************************************************** @@ -47,20 +48,21 @@ program DAMASK_spectral use DAMASK_interface use prec, only: pInt, pReal use IO - use debug, only: debug_Verbosity + use debug, only: debug_verbosity, spectral_debug_verbosity use math use mesh, only: mesh_ipCenterOfGravity use CPFEM, only: CPFEM_general, CPFEM_initAll use FEsolving, only: restartWrite, restartReadSpectral, restartReadStep use numerics, only: err_div_tol, err_stress_tol, err_stress_tolrel , rotation_tol,& - itmax, memory_efficient, DAMASK_NumThreadsInt,& + itmax, memory_efficient, DAMASK_NumThreadsInt, divergence_correction, & fftw_planner_flag, fftw_timelimit use homogenization, only: materialpoint_sizeResults, materialpoint_results !$ use OMP_LIB ! the openMP function library implicit none ! variables to read from loadcase and geom file - real(pReal), dimension(9) :: valueVector ! stores information temporarily from loadcase file + real(pReal), dimension(9) :: temp_valueVector ! stores information temporarily from loadcase file + logical, dimension(9) :: temp_maskVector integer(pInt), parameter :: maxNchunksLoadcase = & (1_pInt + 9_pInt)*3_pInt + & ! deformation, rotation, and stress (1_pInt + 1_pInt)*5_pInt + & ! time, (log)incs, temp, restartfrequency, and outputfrequency @@ -68,27 +70,30 @@ program DAMASK_spectral integer(pInt), dimension (1 + maxNchunksLoadcase*2) :: posLoadcase integer(pInt), parameter :: maxNchunksGeom = 7_pInt ! 4 identifiers, 3 values integer(pInt), dimension (1 + maxNchunksGeom*2) :: posGeom - integer(pInt) :: headerLength,N_l=0_pInt, N_t=0_pInt, N_n=0_pInt, N_Fdot=0_pInt + integer(pInt) :: headerLength, N_l=0_pInt, N_t=0_pInt, N_n=0_pInt, N_Fdot=0_pInt integer(pInt), parameter :: myUnit = 234_pInt character(len=1024) :: path, line, keyword logical :: gotResolution =.false., gotDimension =.false., gotHomogenization = .false. -! variables storing information from loadcase file -!ToDo: create Data Structure loadcase - real(pReal), dimension (:,:,:), allocatable :: bc_deformation, & ! applied velocity gradient or time derivative of deformation gradient - bc_stress, & ! stress BC (if applicable) - bc_rotation ! rotation of BC (if applicable) - real(pReal), dimension(:), allocatable :: bc_timeIncrement, & ! length of increment - bc_temperature ! isothermal starting conditions - integer(pInt), dimension(:), allocatable :: bc_steps, & ! number of steps - bc_outputfrequency, & ! frequency of result writes - bc_restartfrequency, & ! frequency of result writes - bc_logscale ! linear/logaritmic time step flag - logical, dimension(:), allocatable :: bc_followFormerTrajectory,& ! follow trajectory of former loadcase - bc_velGradApplied ! decide wether velocity gradient or fdot is given - logical, dimension(:,:,:,:), allocatable :: bc_mask ! mask of boundary conditions - logical, dimension(:,:,:), allocatable :: bc_maskvector ! linear mask of boundary conditions - character(len=3) :: loadcase_string + 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 ! variables storing information from geom file real(pReal) :: wgt @@ -110,147 +115,154 @@ program DAMASK_spectral integer(pInt) :: size_reduced = 0.0_pReal ! number of stress BCs ! pointwise data - real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold - real(pReal), dimension(:,:,:,:), allocatable :: coordinates - real(pReal), dimension(:,:,:), allocatable :: temperature + real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold + real(pReal), dimension(:,:,:,:), allocatable :: coordinates + real(pReal), dimension(:,:,:), allocatable :: temperature - ! variables storing information for spectral method and FFTW real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field integer(pInt), dimension(3) :: k_s - integer*8, dimension(2) :: fftw_plan ! plans for fftw (forward and backward) + integer*8, dimension(3) :: fftw_plan ! plans for fftw (forward and backward) integer*8 :: fftw_flag ! planner flag for fftw ! loop variables, convergence etc. real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc ! elapsed time, begin of interval, time interval real(pReal) :: guessmode, err_div, err_stress, p_hat_avg - 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 + 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 integer(pInt) :: i, j, k, l, m, n, p - integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode,& + integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode, stepZero=1_pInt, & ierr, notConvergedCounter = 0_pInt, totalStepsCounter = 0_pInt - logical errmatinv + logical :: errmatinv, regrid = .false. real(pReal) :: defgradDet, defgradDetMax, defgradDetMin + 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. !Initializing -!$ 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) ! check for correct number of given arguments + !$ 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 call DAMASK_interface_init() + 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. + !$OMP CRITICAL (write2out) print '(a)', '' - print '(a,a)', '<<<+- DAMASK_spectral init -+>>>' - print '(a,a)', '$Id$' + print '(a,a)', ' <<<+- DAMASK_spectral init -+>>>' + print '(a,a)', ' $Id$' print '(a)', '' - print '(a,a)', 'Working Directory: ',trim(getSolverWorkingDirectoryName()) - print '(a,a)', 'Solver Job Name: ',trim(getSolverJobName()) + print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName()) + print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName()) print '(a)', '' !$OMP END CRITICAL (write2out) ! Reading the loadcase file and allocate variables path = getLoadcaseName() - if (.not. IO_open_file(myUnit,path)) call IO_error(error_ID=30,ext_msg = trim(path)) + if (.not. IO_open_file(myUnit,path)) call IO_error(error_ID=30_pInt,ext_msg = trim(path)) rewind(myUnit) do read(myUnit,'(a1024)',END = 100) line if (IO_isBlank(line)) cycle ! skip empty lines posLoadcase = IO_stringPos(line,maxNchunksLoadcase) - do i = 1, maxNchunksLoadcase, 1 ! reading compulsory parameters for loadcase + do i = 1_pInt, maxNchunksLoadcase, 1_pInt ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,posLoadcase,i))) case('l', 'velocitygrad', 'velgrad','velocitygradient') - N_l = N_l+1 + N_l = N_l + 1_pInt case('fdot') - N_Fdot = N_Fdot+1 + N_Fdot = N_Fdot + 1_pInt case('t', 'time', 'delta') - N_t = N_t+1 + N_t = N_t + 1_pInt case('n', 'incs', 'increments', 'steps', 'logincs', 'logsteps') - N_n = N_n+1 + N_n = N_n + 1_pInt end select enddo ! count all identifiers to allocate memory and do sanity check enddo 100 N_Loadcases = N_n if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check - call IO_error(error_ID=37,ext_msg = trim(path)) ! error message for incomplete loadcase + call IO_error(error_ID=37_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase - allocate (bc_deformation(3,3,N_Loadcases)); bc_deformation = 0.0_pReal - allocate (bc_stress(3,3,N_Loadcases)); bc_stress = 0.0_pReal - allocate (bc_mask(3,3,2,N_Loadcases)); bc_mask = .false. - allocate (bc_maskvector(9,2,N_Loadcases)); bc_maskvector = .false. - allocate (bc_velGradApplied(N_Loadcases)); bc_velGradApplied = .false. - allocate (bc_timeIncrement(N_Loadcases)); bc_timeIncrement = 0.0_pReal - allocate (bc_temperature(N_Loadcases)); bc_temperature = 300.0_pReal - allocate (bc_steps(N_Loadcases)); bc_steps = 0_pInt - allocate (bc_logscale(N_Loadcases)); bc_logscale = 0_pInt - allocate (bc_outputfrequency(N_Loadcases)); bc_outputfrequency = 1_pInt - allocate (bc_restartfrequency(N_Loadcases)); bc_restartfrequency = 1_pInt - allocate (bc_followFormerTrajectory(N_Loadcases)); bc_followFormerTrajectory = .true. - allocate (bc_rotation(3,3,N_Loadcases)); bc_rotation = 0.0_pReal + allocate (bc(N_Loadcases)) rewind(myUnit) loadcase = 0_pInt do read(myUnit,'(a1024)',END = 101) line - if (IO_isBlank(line)) cycle ! skip empty lines - loadcase = loadcase + 1 - bc_rotation(:,:,loadcase) = math_I3 ! assume no rotation, overwrite later in case rotation of loadcase is given + 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 = 300.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 posLoadcase = IO_stringPos(line,maxNchunksLoadcase) - do j = 1,maxNchunksLoadcase + do j = 1_pInt,maxNchunksLoadcase select case (IO_lc(IO_stringValue(line,posLoadcase,j))) - case('fdot','l','velocitygrad','velgrad','velocitygradient') ! assign values for the deformation BC matrix - bc_velGradApplied(loadcase) = (IO_lc(IO_stringValue(line,posLoadcase,j)) == 'l' .or. & - IO_lc(IO_stringValue(line,posLoadcase,j)) == 'velocitygrad') ! in case of given L, set flag to true - valueVector = 0.0_pReal - forall (k = 1:9) bc_maskvector(k,1,loadcase) = IO_stringValue(line,posLoadcase,j+k) /= '*' - do k = 1,9 - if (bc_maskvector(k,1,loadcase)) valueVector(k) = IO_floatValue(line,posLoadcase,j+k) + 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) enddo - bc_mask(:,:,1,loadcase) = transpose(reshape(bc_maskvector(1:9,1,loadcase),(/3,3/))) - bc_deformation(:,:,loadcase) = math_plain9to33(valueVector) + bc(loadcase)%maskDeformation = transpose(reshape(temp_maskVector,(/3,3/))) + bc(loadcase)%deformation = math_plain9to33(temp_valueVector) case('p', 'pk1', 'piolakirchhoff', 'stress') - valueVector = 0.0_pReal - forall (k = 1:9) bc_maskvector(k,2,loadcase) = IO_stringValue(line,posLoadcase,j+k) /= '*' - do k = 1,9 - if (bc_maskvector(k,2,loadcase)) valueVector(k) = IO_floatValue(line,posLoadcase,j+k) ! assign values for the bc_stress matrix + 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 enddo - bc_mask(:,:,2,loadcase) = transpose(reshape(bc_maskvector(1:9,2,loadcase),(/3,3/))) - bc_stress(:,:,loadcase) = math_plain9to33(valueVector) - case('t','time','delta') ! increment time - bc_timeIncrement(loadcase) = IO_floatValue(line,posLoadcase,j+1) - case('temp','temperature') ! starting temperature - bc_temperature(loadcase) = IO_floatValue(line,posLoadcase,j+1) - case('n','incs','increments','steps') ! bc_steps - bc_steps(loadcase) = IO_intValue(line,posLoadcase,j+1) - case('logincs','logsteps') ! true, if log scale - bc_steps(loadcase) = IO_intValue(line,posLoadcase,j+1) - bc_logscale(loadcase) = 1_pInt - case('f','freq','frequency','outputfreq') ! frequency of result writings - bc_outputfrequency(loadcase) = IO_intValue(line,posLoadcase,j+1) - case('r','restart','restartwrite') ! frequency of writing restart information - bc_restartfrequency(loadcase) = IO_intValue(line,posLoadcase,j+1) + 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) case('guessreset','dropguessing') - bc_followFormerTrajectory(loadcase) = .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))) + 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))) case('deg','degree') - p = 1_pInt ! for conversion from degree to radian + p = 1_pInt ! for conversion from degree to radian case('rad','radian') case default - l = 0_pInt ! imediately reading in angles, assuming radians + l = 0_pInt ! imediately reading in angles, assuming radians end select - forall(k = 1:3) temp33_Real(k,1) = IO_floatValue(line,posLoadcase,j +l +k) * real(p,pReal) * inRad - bc_rotation(:,:,loadcase) = math_EulerToR(temp33_Real(:,1)) - case('rotation','rot') ! assign values for the rotation of loadcase matrix - valueVector = 0.0_pReal - forall (k = 1:9) valueVector(k) = IO_floatValue(line,posLoadcase,j+k) - bc_rotation(:,:,loadcase) = math_plain9to33(valueVector) + 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) end select enddo; enddo @@ -260,47 +272,47 @@ program DAMASK_spectral path = getModelName() if (.not. IO_open_file(myUnit,trim(path)//InputFileExtension))& - call IO_error(error_ID=101,ext_msg = trim(path)//InputFileExtension) + call IO_error(error_ID=101_pInt,ext_msg = trim(path)//InputFileExtension) rewind(myUnit) read(myUnit,'(a1024)') line - posGeom = IO_stringPos(line,2) - keyword = IO_lc(IO_StringValue(line,posGeom,2)) + posGeom = IO_stringPos(line,2_pInt) + keyword = IO_lc(IO_StringValue(line,posGeom,2_pInt)) if (keyword(1:4) == 'head') then - headerLength = IO_intValue(line,posGeom,1) + 1_pInt + headerLength = IO_intValue(line,posGeom,1_pInt) + 1_pInt else - call IO_error(error_ID=42) + call IO_error(error_ID=42_pInt) endif rewind(myUnit) - do i = 1, headerLength + do i = 1_pInt, headerLength read(myUnit,'(a1024)') line posGeom = IO_stringPos(line,maxNchunksGeom) select case ( IO_lc(IO_StringValue(line,posGeom,1)) ) case ('dimension') gotDimension = .true. - do j = 2,6,2 + do j = 2_pInt,6_pInt,2_pInt select case (IO_lc(IO_stringValue(line,posGeom,j))) case('x') - geomdimension(1) = IO_floatValue(line,posGeom,j+1) + geomdimension(1) = IO_floatValue(line,posGeom,j+1_pInt) case('y') - geomdimension(2) = IO_floatValue(line,posGeom,j+1) + geomdimension(2) = IO_floatValue(line,posGeom,j+1_pInt) case('z') - geomdimension(3) = IO_floatValue(line,posGeom,j+1) + geomdimension(3) = IO_floatValue(line,posGeom,j+1_pInt) end select enddo case ('homogenization') gotHomogenization = .true. - homog = IO_intValue(line,posGeom,2) + homog = IO_intValue(line,posGeom,2_pInt) case ('resolution') gotResolution = .true. - do j = 2,6,2 + do j = 2_pInt,6_pInt,2_pInt select case (IO_lc(IO_stringValue(line,posGeom,j))) case('a') - resolution(1) = IO_intValue(line,posGeom,j+1) + resolution(1) = IO_intValue(line,posGeom,j+1_pInt) case('b') - resolution(2) = IO_intValue(line,posGeom,j+1) + resolution(2) = IO_intValue(line,posGeom,j+1_pInt) case('c') - resolution(3) = IO_intValue(line,posGeom,j+1) + resolution(3) = IO_intValue(line,posGeom,j+1_pInt) end select enddo case ('picture') @@ -308,226 +320,109 @@ program DAMASK_spectral end select enddo close(myUnit) - if (.not.(gotDimension .and. gotHomogenization .and. gotResolution)) call IO_error(error_ID=45) + if (.not.(gotDimension .and. gotHomogenization .and. gotResolution)) call IO_error(error_ID=45_pInt) if(mod(resolution(1),2_pInt)/=0_pInt .or.& mod(resolution(2),2_pInt)/=0_pInt .or.& - (mod(resolution(3),2_pInt)/=0_pInt .and. resolution(3)/= 1_pInt)) call IO_error(error_ID=103) + (mod(resolution(3),2_pInt)/=0_pInt .and. resolution(3)/= 1_pInt)) call IO_error(error_ID=103_pInt) - 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_temperature(1) ! start out isothermally - allocate (xi (3,resolution(1)/2+1,resolution(2),resolution(3))); xi =0.0_pReal - - wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal) - ! Initialization of CPFEM_general (= constitutive law) and of deformation gradient field - call CPFEM_initAll(bc_temperature(1),1_pInt,1_pInt) + call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) !Output of geom file !$OMP CRITICAL (write2out) print '(a)', '' - print '(a)', '*************************************************************' + print '(a)', '#############################################################' print '(a)', 'DAMASK spectral:' print '(a)', 'The spectral method boundary value problem solver for' print '(a)', 'the Duesseldorf Advanced Material Simulation Kit' - print '(a)', '*************************************************************' + print '(a)', '#############################################################' print '(a,a)', 'Geom File Name: ',trim(path)//'.geom' - 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 + 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 print '(a,i5)','homogenization: ',homog print '(a,L)','spectralPictureMode: ',spectralPictureMode - print '(a)', '************************************************************' + print '(a)', '#############################################################' print '(a,a)','Loadcase File Name: ',trim(getLoadcaseName()) !$OMP END CRITICAL (write2out) - if (bc_followFormerTrajectory(1)) then + if (bc(1)%followFormerTrajectory) then call IO_warning(warning_ID=33_pInt) ! cannot guess along trajectory for first step of first loadcase - bc_followFormerTrajectory(1) = .false. + bc(1)%followFormerTrajectory = .false. endif ! consistency checks and output of loadcase do loadcase = 1_pInt, N_Loadcases !$OMP CRITICAL (write2out) - print '(a)', '-------------------------------------------------------------' + print '(a)', '=============================================================' print '(a,i5)', 'Loadcase: ', loadcase write (loadcase_string, '(i3)' ) loadcase - if (.not. bc_followFormerTrajectory(loadcase)) & + if (.not. bc(loadcase)%followFormerTrajectory) & print '(a)', 'Drop Guessing Along Trajectory' !$OMP END CRITICAL (write2out) - if (any(bc_mask(:,:,1,loadcase) .eqv. bc_mask(1:3,1:3,2,loadcase)))& ! exclusive or masking only - call IO_error(error_ID=31,ext_msg=loadcase_string) - if (any(bc_mask(1:3,1:3,2,loadcase).and.transpose(bc_mask(1:3,1:3,2,loadcase)).and.& !checking if no rotation is allowed by stress BC + 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 reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/))))& - call IO_error(error_ID=38,ext_msg=loadcase_string) - if (bc_velGradApplied(loadcase)) then - do j = 1, 3 - if (any(bc_mask(j,1:3,1,loadcase) .eqv. .true.) .and.& - any(bc_mask(j,1:3,1,loadcase) .eqv. .false.)) call IO_error(error_ID=32,ext_msg=loadcase_string) ! each line should be either fully or not at all defined + 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 enddo !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.6,x)/))','Velocity Gradient:', merge(math_transpose3x3(bc_deformation(1:3,1:3,loadcase)),& - reshape(spread(DAMASK_NaN,1,9),(/3,3/)),& - transpose(bc_mask(1:3,1:3,1,loadcase))) + print '(a)','Velocity Gradient:' !$OMP END CRITICAL (write2out) else !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.6,x)/))','Change of Deformation Gradient:', merge(math_transpose3x3(bc_deformation(1:3,1:3,loadcase)),& - reshape(spread(DAMASK_NaN,1,9),(/3,3/)),& - transpose(bc_mask(1:3,1:3,1,loadcase))) + print '(a)','Change of Deformation Gradient:' !$OMP END CRITICAL (write2out) endif !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.6,x)/))','Stress Boundary Condition/MPa:',merge(math_transpose3x3(bc_stress(1:3,1:3,loadcase)),& - reshape(spread(DAMASK_NaN,1,9),(/3,3/)),& - transpose(bc_mask(:,:,2,loadcase)))*1e-6 + 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 !$OMP END CRITICAL (write2out) - if (any(abs(math_mul33x33(bc_rotation(1:3,1:3,loadcase),math_transpose3x3(bc_rotation(1:3,1:3,loadcase)))-math_I3)& + if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose3x3(bc(loadcase)%rotation))-math_I3)& ! given rotation matrix contains strain >reshape(spread(rotation_tol,1,9),(/3,3/)))& - .or. abs(math_det3x3(bc_rotation(1:3,1:3,loadcase)))>1.0_pReal + rotation_tol) call IO_error(error_ID=46,ext_msg=loadcase_string) + .or. abs(math_det3x3(bc(loadcase)%rotation))>1.0_pReal + rotation_tol) call IO_error(error_ID=46_pInt,ext_msg=loadcase_string) !$OMP CRITICAL (write2out) - if (any(bc_rotation(1:3,1:3,loadcase)/=math_I3)) & - print '(a,/,3(3(f12.6,x)/))','Rotation of BCs:',math_transpose3x3(bc_rotation(1:3,1:3,loadcase)) + if (any(bc(loadcase)%rotation/=math_I3)) & + print '(a,3(3(f12.6,x)/)\)','Rotation of BCs:',math_transpose3x3(bc(loadcase)%rotation) !$OMP END CRITICAL (write2out) - if (bc_timeIncrement(loadcase) < 0.0_pReal) call IO_error(error_ID=34,ext_msg=loadcase_string) ! negative time increment + if (bc(loadcase)%timeIncrement < 0.0_pReal) call IO_error(error_ID=34_pInt,ext_msg=loadcase_string) ! negative time increment !$OMP CRITICAL (write2out) - print '(a,f12.6)','Temperature: ',bc_temperature(loadcase) - print '(a,f12.6)','Time: ',bc_timeIncrement(loadcase) + print '(a,f12.6)','Temperature:',bc(loadcase)%temperature + print '(a,f12.6)','Time: ',bc(loadcase)%timeIncrement !$OMP END CRITICAL (write2out) - if (bc_steps(loadcase) < 1_pInt) call IO_error(error_ID=35,ext_msg=loadcase_string) ! non-positive increment count + if (bc(loadcase)%steps < 1_pInt) call IO_error(error_ID=35_pInt,ext_msg=loadcase_string) ! non-positive increment count !$OMP CRITICAL (write2out) - print '(a,i5)','Increments: ',bc_steps(loadcase) + print '(a,i5)','Steps: ',bc(loadcase)%steps !$OMP END CRITICAL (write2out) - if (bc_outputfrequency(loadcase) < 1_pInt) call IO_error(error_ID=36,ext_msg=loadcase_string) ! non-positive result frequency + if (bc(loadcase)%outputfrequency < 1_pInt) call IO_error(error_ID=36_pInt,ext_msg=loadcase_string) ! non-positive result frequency !$OMP CRITICAL (write2out) - print '(a,i5)','Freq. of Restults Output: ',bc_outputfrequency(loadcase) + print '(a,i5)','Freq. of Results Output: ',bc(loadcase)%outputfrequency !$OMP END CRITICAL (write2out) - if (bc_restartfrequency(loadcase) < 1_pInt) call IO_error(error_ID=39,ext_msg=loadcase_string) ! non-positive restart frequency + if (bc(loadcase)%restartfrequency < 1_pInt) call IO_error(error_ID=39_pInt,ext_msg=loadcase_string) ! non-positive restart frequency !$OMP CRITICAL (write2out) - print '(a,i5)','Freq. of Restart Information Output: ',bc_restartfrequency(loadcase) + print '(a,i5)','Freq. of Restart Information Output: ',bc(loadcase)%restartfrequency !$OMP END CRITICAL (write2out) enddo - if (.not. restartReadSpectral) then ! no deformation at the beginning - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, 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 - loadcase = 1_pInt - step = 1_pInt - else ! using old values - if (IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',trim(getModelName()),size(defgrad))) then - read (777,rec=1) defgrad - close (777) - endif - defgradold = defgrad - defgradAim = 0.0_pReal - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) - defgradAim = defgradAim + defgrad(i,j,k,1:3,1:3) - enddo; enddo; enddo - defgradAim = defgradAim * wgt - defgradAimOld = defgradAim - do i = 1_pInt, N_loadcases ! looping over ALL loadcase - time0 = time ! loadcase start time - timeinc = bc_timeIncrement(i)/bc_steps(i) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used - do j = 1_pInt, bc_steps(i) ! looping over ALL steps in current loadcase - if (totalStepsCounter <= restartReadStep) then ! forwarding to restart step - totalStepsCounter = totalStepsCounter + 1_pInt - if (bc_logscale(i) == 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_timeIncrement(1)*(2.0**(1 - bc_steps(1))) ! assume 1st step is equal to 2nd - else ! not-1st step of 1st loadcase of loglinear scale - timeinc = bc_timeIncrement(1)*(2.0**(j - (1 + bc_steps(1)))) - endif - else ! not-1st loadcase of loglinear scale - timeinc = time0 * ( ((1.0_pReal+bc_timeIncrement(i)/time0)**(float( j )/(bc_steps(i)))) & - - ((1.0_pReal+bc_timeIncrement(i)/time0)**(float((j-1))/(bc_steps(i)))) ) - endif - endif - time = time + timeinc - step = j - loadcase = i - endif - enddo - enddo - totalStepsCounter = totalStepsCounter - 1_pInt - endif - ielem = 0_pInt - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) - ielem = ielem +1 - 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!!! - call CPFEM_general(2,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_rotation(1:3,1:3,loadcase)) ! rotate_forward: lab -> load system - - if (debug_verbosity > 1) then - !$OMP CRITICAL (write2out) - write (6,*) 'First Call to CPFEM_general finished' - !$OMP END CRITICAL (write2out) - endif - - do k = 1, resolution(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - k_s(3) = k-1 - if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) - do j = 1, resolution(2) - k_s(2) = j-1 - if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) - do i = 1, resolution(1)/2+1 - k_s(1) = i-1 - xi(3,i,j,k) = 0.0_pReal ! 2D case - if(resolution(3) > 1) 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,resolution(3); do j = 1,resolution(2); do i = 1,resolution(1)/2+1 - if(k==resolution(3)/2+1) xi(3,i,j,k)= 0.0_pReal - if(j==resolution(2)/2+1) xi(2,i,j,k)= 0.0_pReal - if(i==resolution(1)/2+1) 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+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 - if (any(xi(:,i,j,k) /= 0.0_pReal)) then - do l = 1,3; do m = 1,3 - 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,3; do m=1,3; do n=1,3; do p=1,3 - 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 - - allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal - ! Initialization of fftw (see manual on fftw.org for more details) #ifdef _OPENMP if(DAMASK_NumThreadsInt>0_pInt) then call dfftw_init_threads(ierr) - if(ierr == 0_pInt) call IO_error(error_ID=104) + if(ierr == 0_pInt) call IO_error(error_ID=104_pInt) call dfftw_plan_with_nthreads(DAMASK_NumThreadsInt) endif #endif - !is not working, have to find out how it is working in FORTRAN - !call dfftw_timelimit(fftw_timelimit) + !call dfftw_timelimit(fftw_timelimit) !is not working, have to fix it in FFTW source file - ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f - ! ordered from slow execution (but fast plan creation) to fast execution - select case(IO_lc(fftw_planner_flag)) - case('estimate','fftw_estimate') + 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 fftw_flag = 64 case('measure','fftw_measure') fftw_flag = 0 @@ -536,94 +431,243 @@ program DAMASK_spectral case('exhaustive','fftw_exhaustive') fftw_flag = 8 case default - !$OMP CRITICAL (write2out) - write (6,*) 'No valid parameter for FFTW given, using FFTW_PATIENT' - !$OMP END CRITICAL (write2out) + call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_planner_flag))) fftw_flag = 32 end select - - call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/resolution(1),resolution(2),resolution(3)/),9,& - workfft,(/resolution(1) +2,resolution(2),resolution(3)/),1,(resolution(1) +2)*resolution(2)*resolution(3),& - workfft,(/resolution(1)/2+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*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+1,resolution(2),resolution(3)/),1,(resolution(1)/2+1)*resolution(2)*resolution(3),& - workfft,(/resolution(1) +2,resolution(2),resolution(3)/),1,(resolution(1) +2)*resolution(2)*resolution(3),fftw_flag) - !$OMP CRITICAL (write2out) - if (debug_verbosity > 1) then - write (6,*) 'FFTW initialized' - endif + + 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 -! write header of output file - open(538,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())& - //'.spectralOut',form='UNFORMATTED',status='REPLACE')!,access='DIRECT') - 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_logscale ! one entry per loadcase (0: linear, 1: log) - write(538), 'frequencies', bc_outputfrequency ! one entry per loadcase - write(538), 'times', bc_timeIncrement ! one entry per loadcase - bc_steps(1)= bc_steps(1) + 1_pInt - write(538), 'increments', bc_steps ! one entry per loadcase ToDo: rename keyword to steps - bc_steps(1)= bc_steps(1) - 1_pInt - write(538), 'startingIncrement', totalStepsCounter - write(538), 'eoh' ! end of header - write(538), materialpoint_results(:,1,:) ! initial (non-deformed) results -!$OMP END CRITICAL (write2out) -! Initialization done - + endif +print*, totalStepsCounter +print*, loadcase +print*, step +pause !************************************************************* ! Loop over loadcases defined in the loadcase file do loadcase = loadcase, N_Loadcases !************************************************************* time0 = time ! loadcase start time - if (bc_followFormerTrajectory(loadcase)) then ! continue to guess along former trajectory where applicable + if (bc(loadcase)%followFormerTrajectory) 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 step endif - mask_defgrad = merge(ones,zeroes,bc_mask(:,:,1,loadcase)) - mask_stress = merge(ones,zeroes,bc_mask(:,:,2,loadcase)) - size_reduced = count(bc_maskvector(1:9,2,loadcase)) + mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation) + mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress) + size_reduced = count(bc(loadcase)%maskStressVector) allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal - timeinc = bc_timeIncrement(loadcase)/bc_steps(loadcase) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used - fDot = bc_deformation(:,:,loadcase) ! only valid for given fDot. will be overwritten later in case L is given + 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 !************************************************************* ! loop oper steps defined in input file for current loadcase - do step = step, bc_steps(loadcase) + do step = step, bc(loadcase)%steps !************************************************************* - if (mod(step - 1_pInt,bc_restartFrequency(loadcase))==0_pInt) then ! setting restart parameter for FEsolving - restartWrite = .true. - if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then - write (777,rec=1) defgrad +!************************************************************* +! 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 close (777) endif - else - restartWrite = .false. + 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 endif - if (bc_logscale(loadcase) == 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_timeIncrement(1)*(2.0**(1 - bc_steps(1))) ! assume 1st step is equal to 2nd - else ! not-1st step of 1st loadcase of loglinear scale - timeinc = bc_timeIncrement(1)*(2.0**(step - (1 + bc_steps(1)))) - endif - else ! not-1st loadcase of loglinear scale - timeinc = time0 * ( ((1.0_pReal+bc_timeIncrement(loadcase)/time0)**(float( step )/(bc_steps(loadcase)))) & - - ((1.0_pReal+bc_timeIncrement(loadcase)/time0)**(float((step-1))/(bc_steps(loadcase)))) ) - endif + + 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 - time = time + timeinc + + 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) ) + endif + endif + time = time + timeinc - if (bc_velGradApplied(loadcase)) & ! calculate fDot from given L and current F - fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase), defgradAim) + if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F + fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim) !winding forward of deformation aim in loadcase system temp33_Real = defgradAim @@ -633,23 +677,23 @@ program DAMASK_spectral defgradAimOld = temp33_Real ! update local deformation gradient - if (any(bc_rotation(1:3,1:3,loadcase)/=math_I3)) then ! lab and loadcase coordinate system are NOT the same - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + 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) temp33_Real = defgrad(i,j,k,1:3,1:3) - if (bc_velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing) - fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase),& - math_rotate_forward3x3(defgradold(i,j,k,1:3,1:3),bc_rotation(1:3,1:3,loadcase))) + if (bc(loadcase)%velGradApplied) & ! use velocity gradient to calculate new deformation gradient (if not guessing) + fDot = math_mul33x33(bc(loadcase)%deformation,& + math_rotate_forward3x3(defgradold(i,j,k,1:3,1:3),bc(loadcase)%rotation)) defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon + guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing... + math_rotate_backward3x3((1.0_pReal-guessmode) * mask_defgrad * fDot,& - bc_rotation(1:3,1:3,loadcase)) *timeinc ! apply the prescribed value where deformation is given if not guessing + bc(loadcase)%rotation) *timeinc ! apply the prescribed value where deformation is given if not guessing defgradold(i,j,k,1:3,1:3) = temp33_Real enddo; enddo; enddo - else ! one coordinate system for lab and loadcase, save some multiplication - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + 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) temp33_Real = defgrad(i,j,k,1:3,1:3) - if (bc_velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing) - fDot = math_mul33x33(bc_deformation(1:3,1:3,loadcase),defgradold(i,j,k,1:3,1:3)) + if (bc(loadcase)%velGradApplied) & ! use velocity gradient to calculate new deformation gradient (if not guessing) + fDot = math_mul33x33(bc(loadcase)%deformation,defgradold(i,j,k,1:3,1:3)) defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon + guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing... + (1.0_pReal-guessmode) * mask_defgrad * fDot * timeinc ! apply the prescribed value where deformation is given if not guessing @@ -665,12 +709,12 @@ program DAMASK_spectral if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied c_prev99 = math_Plain3333to99(c_prev) k = 0_pInt ! build reduced stiffness - do n = 1,9 - if(bc_maskvector(n,2,loadcase)) then + do n = 1_pInt,9_pInt + if(bc(loadcase)%maskStressVector(n)) then k = k + 1_pInt j = 0_pInt - do m = 1,9 - if(bc_maskvector(m,2,loadcase)) then + do m = 1_pInt,9_pInt + if(bc(loadcase)%maskStressVector(m)) then j = j + 1_pInt c_reduced(k,j) = c_prev99(n,m) endif; enddo; endif; enddo @@ -678,112 +722,123 @@ program DAMASK_spectral if(errmatinv) call IO_error(error_ID=800) s_prev99 = 0.0_pReal ! build full compliance k = 0_pInt - do n = 1,9 - if(bc_maskvector(n,2,loadcase)) then + do n = 1_pInt,9_pInt + if(bc(loadcase)%maskStressVector(n)) then k = k + 1_pInt j = 0_pInt - do m = 1,9 - if(bc_maskvector(m,2,loadcase)) then + do m = 1_pInt,9_pInt + if(bc(loadcase)%maskStressVector(m)) then j = j + 1_pInt s_prev99(n,m) = s_reduced(k,j) endif; enddo; endif; enddo s_prev = (math_Plain99to3333(s_prev99)) endif + + !$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) + !************************************************************* ! convergence loop do while(iter < itmax .and. & (err_div > err_div_tol .or. & err_stress > err_stress_tol)) iter = iter + 1_pInt - !$OMP CRITICAL (write2out) - print '(A)', '************************************************************' - print '(3(A,I5.5,tr2)A)', '**** Loadcase = ',loadcase, 'Step = ',step, 'Iteration = ',iter,'****' - print '(A)', '************************************************************' - !$OMP END CRITICAL (write2out) - if (restartWrite .eq. .true. .and. iter == 1_pInt) then - !$OMP CRITICAL (write2out) - print '(A)', 'Writing converged Results of previous Step for Restart' - print '(A)', '************************************************************' - !$OMP END CRITICAL (write2out) - endif - workfft = 0.0_pReal ! needed because of the padding for FFTW !************************************************************* - do n = 1,3; do m = 1,3 - defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt + 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 enddo; enddo !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.7,x)/))', 'Deformation Gradient:',math_transpose3x3(defgrad_av) - - print '(A,/)', '== Update Stress Field (Constitutive Evaluation P(F)) ======' + print '(a,/,3(3(f12.7,x)/)\)', 'Deformation Gradient:',math_transpose3x3(defgrad_av) + print '(A)', '... Update Stress Field (Constitutive Evaluation P(F)) ......' !$OMP END CRITICAL (write2out) ielem = 0_pInt - defgradDetMax = 0.0_pReal - defgradDetMin = 999.0_pReal - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + defgradDetMax = -999.0_pReal + defgradDetMin = 999.0_pReal + do k = 1_pInt, resolution(3); do j = 1_pInt, resolution(2); do i = 1_pInt, resolution(1) defgradDet = math_det3x3(defgrad(i,j,k,1:3,1:3)) defgradDetMax = max(defgradDetMax,defgradDet) defgradDetMin = min(defgradDetMin,defgradDet) - ielem = ielem + 1 - call CPFEM_general(3,& ! collect cycle - coordinates(1:3,i,j,k), defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& + ielem = ielem + 1_pInt + call CPFEM_general(3_pInt,& ! collect cycle + coordinates(1:3,i,j,k), defgradold(i,j,k,1:3,1:3), defgrad(i,j,k,1:3,1:3),& temperature(i,j,k),timeinc,ielem,1_pInt,& cstress,dsde, pstress, dPdF) enddo; enddo; enddo - + + 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 c_current = 0.0_pReal ielem = 0_pInt - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) + do k = 1_pInt, resolution(3); do j = 1_pInt, resolution(2); do i = 1_pInt, resolution(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, coordinates(1:3,i,j,k),& - defgradold(i,j,k,:,:), defgrad(i,j,k,:,:),& ! others get 2 (saves winding forward effort) + defgradold(i,j,k,1:3,1:3), defgrad(i,j,k,1:3,1:3),& ! others get 2 (saves winding forward effort) temperature(i,j,k),timeinc,ielem,1_pInt,& cstress,dsde, pstress,dPdF) CPFEM_mode = 2_pInt - workfft(i,j,k,:,:) = pstress ! build up average P-K stress + workfft(i,j,k,1:3,1:3) = pstress ! build up average P-K stress c_current = c_current + dPdF enddo; enddo; enddo - restartWrite = .false. ! ToDo: don't know if we need it - do n = 1,3; do m = 1,3 - pstress_av(m,n) = sum(workfft(1:resolution(1),:,:,m,n)) * wgt + 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 enddo; enddo !$OMP CRITICAL (write2out) - print '(a,/,3(3(f12.7,x)/))', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 + print '(a,/,3(3(f12.7,x)/)\)', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 err_stress_tol = 0.0_pReal - pstress_av_load = math_rotate_forward3x3(pstress_av,bc_rotation(1:3,1:3,loadcase)) - if(size_reduced > 0_pInt) then ! calculate stress BC if applied - err_stress = maxval(abs(mask_stress * (pstress_av_load - bc_stress(1:3,1:3,loadcase)))) ! 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,E10.5)/)', 'Error Stress = ',err_stress, ', Tol. = ', err_stress_tol - defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av_load - bc_stress(1:3,1:3,loadcase)))) ! residual on given stress components + 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 defgradAim = defgradAim + defgradAimCorr - print '(a,/,3(3(f12.7,x)/))', 'Deformation Aim: ',math_transpose3x3(math_rotate_backward3x3(& - defgradAim,bc_rotation(1:3,1:3,loadcase))) - print '(a,x,f12.7,/)' , 'Determinant of Deformation Aim: ', math_det3x3(defgradAim) - print '(a,x,f12.7,/)' , 'Volume Change Max: ', defgradDetMax - print '(a,x,f12.7,/)' , 'Volume Change Min: ', defgradDetMin + 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) endif - print '(A,/)', '== Calculating Equilibrium Using Spectral Method ===========' + print '(A)', '... Calculating Equilibrium Using Spectral Method ...........' !$OMP END CRITICAL (write2out) - call dfftw_execute_dft_r2c(fftw_plan(1),workfft,workfft) ! FFT of pstress + call dfftw_execute_dft_r2c(fftw_plan(1),workfft,workfft) ! FFT of pstress p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(workfft(1,1,1,1:3,1:3),& ! L_2 norm of average stress in fourier space, math_transpose3x3(workfft(1,1,1,1:3,1:3)))))) ! ignore imaginary part as it is always zero for real only input)) err_div = 0.0_pReal - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 - err_div = err_div + sqrt(sum((math_mul33x3_complex(workfft(i*2-1,j,k,1:3,1:3)+& ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain) - workfft(i*2, j,k,1:3,1:3)*img,xi(1:3,i,j,k)))**2.0)) + 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)))) enddo; enddo; enddo - - err_div = err_div*wgt/p_hat_avg*(minval(geomdimension)*wgt**(-1/4)) ! weigthting, multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency + 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 if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat - do k = 1, resolution(3); do j = 1, resolution(2) ;do i = 1, resolution(1)/2+1 + 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,3; do m = 1,3 + 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)) @@ -791,72 +846,89 @@ program DAMASK_spectral xiDyad = 0.0_pReal temp33_Real = 0.0_pReal endif - do l=1,3; do m=1,3; do n=1,3; do p=1,3 + do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt; do p=1_pInt,3_pInt gamma_hat(1,1,1, l,m,n,p) = - 0.25_pReal*(temp33_Real(l,n)+temp33_Real(n,l))*& (xiDyad(m,p) +xiDyad(p,m)) enddo; enddo; enddo; enddo - do m = 1,3; do n = 1,3 - temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,:,:) *(workfft(i*2-1,j,k,:,:)& - +workfft(i*2 ,j,k,:,:)*img)) + do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt + temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,1:3,1:3) *(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3)& + +workfft(i*2_pInt ,j,k,1:3,1:3)*img)) enddo; enddo - workfft(i*2-1,j,k,:,:) = real (temp33_Complex) - workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) + workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) = real (temp33_Complex) + workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex) enddo; enddo; enddo else ! use precalculated gamma-operator - do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 - do m = 1,3; do n = 1,3 - temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,:,:) *(workfft(i*2-1,j,k,:,:)& - + workfft(i*2 ,j,k,:,:)*img)) + 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)) enddo; enddo - workfft(i*2-1,j,k,:,:) = real (temp33_Complex) - workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) + workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) = real (temp33_Complex) + workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex) enddo; enddo; enddo endif + if(debugDivergence) then + 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 ! average strain - workfft(1,1,1,:,:) = defgrad_av - math_I3 ! zero frequency (real part) - workfft(2,1,1,:,:) = 0.0_pReal ! zero frequency (imaginary part) + workfft(1,1,1,1:3,1:3) = defgrad_av - math_I3 ! zero frequency (real part) + workfft(2,1,1,1:3,1:3) = 0.0_pReal ! zero frequency (imaginary part) call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft) defgrad = defgrad + workfft(1:resolution(1),:,:,:,:)*wgt do m = 1,3; do n = 1,3 defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt enddo; enddo - defgradAim_lab = math_rotate_backward3x3(defgradAim,bc_rotation(1:3,1:3,loadcase)) + defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation) 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 enddo; enddo !$OMP CRITICAL (write2out) - print '(2(a,E10.5)/)', 'Error Divergence = ',err_div, ', Tol. = ', err_div_tol + print '(2(a,es10.4))', 'Error Divergence = ',err_div, ', Tol. = ', err_div_tol !$OMP END CRITICAL (write2out) enddo ! end looping when convergency is achieved - c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc_rotation(1:3,1:3,loadcase)) ! calculate stiffness for next step + c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness for next step !ToDo: Incfluence for next loadcase - if (mod(totalStepsCounter,bc_outputfrequency(loadcase)) == 0_pInt) then ! at output frequency - write(538), materialpoint_results(:,1,:) ! write result to file - endif - totalStepsCounter = totalStepsCounter + 1_pInt !$OMP CRITICAL (write2out) + print '(a)', '=============================================================' if(err_div<=err_div_tol .and. err_stress<=err_stress_tol) then - print '(3(A,I5.5),A,/)', '== Step ',step, ' of Loadcase ',loadcase,' (Total ', totalStepsCounter,') Converged ====' + print '(A,I5.5,A)', 'Increment ', totalStepsCounter, ' Converged' else - print '(3(A,I5.5),A,/)', '== Step ',step, ' of Loadcase ',loadcase,' (Total ', totalStepsCounter,') NOT Converged ' + print '(A,I5.5,A)', 'Increment ', totalStepsCounter, ' NOT Converged' notConvergedCounter = notConvergedCounter + 1 endif + 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 !$OMP END CRITICAL (write2out) enddo ! end looping over steps in current loadcase deallocate(c_reduced) deallocate(s_reduced) + step = 1_pInt ! Reset Step Counter enddo ! end looping over loadcases !$OMP CRITICAL (write2out) - print '(A,/)', '############################################################' - print '(a,i5.5,a,i5.5,a)', 'Of ', totalStepsCounter -restartReadStep, ' Total Steps, ', notConvergedCounter, ' Steps did not Converge!' + print '(a)', '#############################################################' + print '(a,i5.5,a,i5.5,a)', 'Of ', totalStepsCounter -restartReadStep, ' Calculated Steps, ', notConvergedCounter, ' Steps did not Converge!' !$OMP END CRITICAL (write2out) close(538) call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2)) - + if(debugDivergence) call dfftw_destroy_plan(fftw_plan(3)) end program DAMASK_spectral !******************************************************************** diff --git a/code/IO.f90 b/code/IO.f90 index e59afb2e3..2259ad2a5 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1416,6 +1416,8 @@ endfunction select case (warning_ID) case (33_pInt) msg = 'cannot guess along trajectory for first step of first loadcase' + case (47_pInt) + msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' case (101_pInt) msg = '+ crystallite debugging off... +' case (600_pInt) diff --git a/code/config/debug.config b/code/config/debug.config index 8f0a28b03..87698a6a2 100644 --- a/code/config/debug.config +++ b/code/config/debug.config @@ -15,3 +15,7 @@ selective 1 # >0 true to switch on e,i,g selective deb element 1 # selected element for debugging (synonymous: "el", "e") ip 1 # selected integration point for debugging (synonymous: "integrationpoint", "i") grain 1 # selected grain at ip for debugging (synonymous: "gr", "g") + +### spectral solver debugging parameters ### +generalDebugSpectral 0 # > 0: general (algorithmical) debug outputs +divergenceDebugSpectral 0 # > 0: calculate more divergence measures and print them out \ No newline at end of file diff --git a/code/debug.f90 b/code/debug.f90 index be7b2b68a..49019f469 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -51,6 +51,7 @@ real(pReal) :: debug_jacobianMax real(pReal) :: debug_jacobianMin logical :: debug_selectiveDebugger = .true. integer(pInt) :: debug_verbosity = 1_pInt +integer(pInt) :: spectral_debug_verbosity = 0_pInt CONTAINS @@ -122,6 +123,12 @@ subroutine debug_init() debug_selectiveDebugger = IO_intValue(line,positions,2) > 0_pInt case ('verbosity') debug_verbosity = IO_intValue(line,positions,2) + case ('generaldebugspectral') ! use bitwise logical and, continue with +8_pInt + if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 1_pInt + case ('divergencedebugspectral') + if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 2_pInt + case ('restartdebugspectral') + if(IO_intValue(line,positions,2)) spectral_debug_verbosity = spectral_debug_verbosity + 4_pInt endselect enddo 100 close(fileunit) @@ -164,7 +171,11 @@ subroutine debug_init() debug_i = 0_pInt debug_g = 0_pInt endif - + !$OMP CRITICAL (write2out) ! bitwise coded + if (iand(spectral_debug_verbosity,1_pInt)==1_pInt) write(6,'(a)') ' Spectral General Debugging' + if (iand(spectral_debug_verbosity,2_pInt)==2_pInt) write(6,'(a)') ' Spectral Divergence Debugging' + if (iand(spectral_debug_verbosity,4_pInt)==4_pInt) write(6,'(a)') ' Spectral Restart Debugging' + !$OMP END CRITICAL (write2out) endsubroutine @@ -196,6 +207,7 @@ subroutine debug_reset() debug_stressMin = huge(1.0_pReal) debug_jacobianMax = -huge(1.0_pReal) debug_jacobianMin = huge(1.0_pReal) + spectral_debug_verbosity = 0.0_pReal endsubroutine diff --git a/code/numerics.f90 b/code/numerics.f90 index 71502b306..1d758e86c 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -70,9 +70,10 @@ real(pReal) :: relevantStrain, & ! strain err_stress_tolrel, & ! factor to multiply with highest stress to get err_stress_tol fftw_timelimit, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org rotation_tol ! tolerance of rotation specified in loadcase -character(len=64) fftw_planner_flag ! sets the planig-rigor flag, see manual on www.fftw.org -logical memory_efficient ! for fast execution (pre calculation of gamma_hat) -integer(pInt) itmax , & ! maximum number of iterations +character(len=64) :: fftw_planner_flag ! sets the planig-rigor flag, see manual on www.fftw.org +logical :: memory_efficient,& ! for fast execution (pre calculation of gamma_hat) + divergence_correction ! correct divergence calculation in fourier space +integer(pInt) :: itmax , & ! maximum number of iterations !* Random seeding parameters @@ -172,8 +173,8 @@ subroutine numerics_init() fftw_timelimit = -1.0_pReal ! no timelimit of plan creation for FFTW fftw_planner_flag ='FFTW_PATIENT' rotation_tol = 1.0e-12 - -!* Random seeding parameters: added <<>> + divergence_correction = .true. +!* Random seeding parameters fixedSeed = 0_pInt @@ -290,6 +291,8 @@ subroutine numerics_init() fftw_planner_flag = IO_stringValue(line,positions,2) case ('rotation_tol') rotation_tol = IO_floatValue(line,positions,2) + case ('divergence_correction') + divergence_correction = IO_intValue(line,positions,2) > 0_pInt !* Random seeding parameters case ('fixed_seed') @@ -310,69 +313,64 @@ subroutine numerics_init() ! writing parameters to output file !$OMP CRITICAL (write2out) - write(6,'(a24,x,e8.1)') 'relevantStrain: ',relevantStrain - write(6,'(a24,x,e8.1)') 'defgradTolerance: ',defgradTolerance - write(6,'(a24,x,i8)') 'iJacoStiffness: ',iJacoStiffness - write(6,'(a24,x,i8)') 'iJacoLpresiduum: ',iJacoLpresiduum - write(6,'(a24,x,e8.1)') 'pert_Fg: ',pert_Fg - write(6,'(a24,x,i8)') 'pert_method: ',pert_method - write(6,'(a24,x,i8)') 'nCryst: ',nCryst - write(6,'(a24,x,e8.1)') 'subStepMinCryst: ',subStepMinCryst - write(6,'(a24,x,e8.1)') 'subStepSizeCryst: ',subStepSizeCryst - write(6,'(a24,x,e8.1)') 'stepIncreaseCryst: ',stepIncreaseCryst - write(6,'(a24,x,i8)') 'nState: ',nState - write(6,'(a24,x,i8)') 'nStress: ',nStress - write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState - write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature - write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress - write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress - write(6,'(a24,2(x,i8))')'integrator: ',numerics_integrator - write(6,*) + write(6,'(a24,x,e8.1)') ' relevantStrain: ',relevantStrain + write(6,'(a24,x,e8.1)') ' defgradTolerance: ',defgradTolerance + write(6,'(a24,x,i8)') ' iJacoStiffness: ',iJacoStiffness + write(6,'(a24,x,i8)') ' iJacoLpresiduum: ',iJacoLpresiduum + write(6,'(a24,x,e8.1)') ' pert_Fg: ',pert_Fg + write(6,'(a24,x,i8)') ' pert_method: ',pert_method + write(6,'(a24,x,i8)') ' nCryst: ',nCryst + write(6,'(a24,x,e8.1)') ' subStepMinCryst: ',subStepMinCryst + write(6,'(a24,x,e8.1)') ' subStepSizeCryst: ',subStepSizeCryst + write(6,'(a24,x,e8.1)') ' stepIncreaseCryst: ',stepIncreaseCryst + write(6,'(a24,x,i8)') ' nState: ',nState + write(6,'(a24,x,i8)') ' nStress: ',nStress + write(6,'(a24,x,e8.1)') ' rTol_crystalliteState: ',rTol_crystalliteState + write(6,'(a24,x,e8.1)') ' rTol_crystalliteTemp: ',rTol_crystalliteTemperature + write(6,'(a24,x,e8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress + write(6,'(a24,x,e8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress + write(6,'(a24,2(x,i8),/)')' integrator: ',numerics_integrator - write(6,'(a24,x,i8)') 'nHomog: ',nHomog - write(6,'(a24,x,e8.1)') 'subStepMinHomog: ',subStepMinHomog - write(6,'(a24,x,e8.1)') 'subStepSizeHomog: ',subStepSizeHomog - write(6,'(a24,x,e8.1)') 'stepIncreaseHomog: ',stepIncreaseHomog - write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate - write(6,*) + write(6,'(a24,x,i8)') ' nHomog: ',nHomog + write(6,'(a24,x,e8.1)') ' subStepMinHomog: ',subStepMinHomog + write(6,'(a24,x,e8.1)') ' subStepSizeHomog: ',subStepSizeHomog + write(6,'(a24,x,e8.1)') ' stepIncreaseHomog: ',stepIncreaseHomog + write(6,'(a24,x,i8,/)') ' nMPstate: ',nMPstate !* RGC parameters - write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC - write(6,'(a24,x,e8.1)') 'rTol_RGC: ',relTol_RGC - write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_RGC - write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC - write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC - write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC - write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',viscPower_RGC - write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModus_RGC - write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC - write(6,'(a24,x,e8.1)') 'maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC - write(6,'(a24,x,e8.1)') 'volDiscrepancyMod_RGC: ',volDiscrMod_RGC - write(6,'(a24,x,e8.1)') 'discrepancyPower_RGC: ',volDiscrPow_RGC - write(6,*) + write(6,'(a24,x,e8.1)') ' aTol_RGC: ',absTol_RGC + write(6,'(a24,x,e8.1)') ' rTol_RGC: ',relTol_RGC + write(6,'(a24,x,e8.1)') ' aMax_RGC: ',absMax_RGC + write(6,'(a24,x,e8.1)') ' rMax_RGC: ',relMax_RGC + write(6,'(a24,x,e8.1)') ' perturbPenalty_RGC: ',pPert_RGC + write(6,'(a24,x,e8.1)') ' relevantMismatch_RGC: ',xSmoo_RGC + write(6,'(a24,x,e8.1)') ' viscosityrate_RGC: ',viscPower_RGC + write(6,'(a24,x,e8.1)') ' viscositymodulus_RGC: ',viscModus_RGC + write(6,'(a24,x,e8.1)') ' maxrelaxation_RGC: ',maxdRelax_RGC + write(6,'(a24,x,e8.1)') ' maxVolDiscrepancy_RGC: ',maxVolDiscr_RGC + write(6,'(a24,x,e8.1)') ' volDiscrepancyMod_RGC: ',volDiscrMod_RGC + write(6,'(a24,x,e8.1,/)') ' discrepancyPower_RGC: ',volDiscrPow_RGC !* spectral parameters - write(6,'(a24,x,e8.1)') 'err_div_tol: ',err_div_tol - write(6,'(a24,x,e8.1)') 'err_stress_tolrel: ',err_stress_tolrel - write(6,'(a24,x,i8)') 'itmax: ',itmax - write(6,'(a24,x,L8)') 'memory_efficient: ',memory_efficient + write(6,'(a24,x,e8.1)') ' err_div_tol: ',err_div_tol + write(6,'(a24,x,e8.1)') ' err_stress_tolrel: ',err_stress_tolrel + write(6,'(a24,x,i8)') ' itmax: ',itmax + write(6,'(a24,x,L8)') ' memory_efficient: ',memory_efficient if(fftw_timelimit<0) then - write(6,'(a24,x,L8)') 'fftw_timelimit: ',.false. + write(6,'(a24,x,L8)') ' fftw_timelimit: ',.false. else - write(6,'(a24,x,e8.1)') 'fftw_timelimit: ',fftw_timelimit + write(6,'(a24,x,e8.1)') ' fftw_timelimit: ',fftw_timelimit endif - write(6,'(a24,x,a)') 'fftw_planner_flag: ',trim(fftw_planner_flag) - write(6,'(a24,x,e8.1)') 'rotation_tol: ',rotation_tol - write(6,*) + write(6,'(a24,x,a)') ' fftw_planner_flag: ',trim(fftw_planner_flag) + write(6,'(a24,x,e8.1)') ' rotation_tol: ',rotation_tol + write(6,'(a24,x,L8,/)') ' divergence_correction: ',divergence_correction !* Random seeding parameters - write(6,'(a24,x,i16)') 'fixed_seed: ',fixedSeed - write(6,*) + write(6,'(a24,x,i16,/)') ' fixed_seed: ',fixedSeed !$OMP END CRITICAL (write2out) !* openMP parameter -!$ write(6,'(a24,x,i8)') 'number of threads: ',DAMASK_NumThreadsInt -!$ write(6,*) +!$ write(6,'(a24,x,i8,/)') ' number of threads: ',DAMASK_NumThreadsInt ! sanity check if (relevantStrain <= 0.0_pReal) call IO_error(260)