From 13b55275b1583e7d6d9b4168bd9b082dc762437e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 24 Oct 2012 11:31:40 +0000 Subject: [PATCH] documented utilities and structured, worked on the restart capabilities of the new basic solver --- code/CPFEM.f90 | 2 +- code/DAMASK_spectral.f90 | 18 +- code/DAMASK_spectral_driver.f90 | 133 ++-- code/DAMASK_spectral_solverAL.f90 | 4 +- code/DAMASK_spectral_solverBasic.f90 | 90 ++- code/DAMASK_spectral_solverBasicPETSc.f90 | 4 +- code/DAMASK_spectral_utilities.f90 | 803 ++++++++++++---------- 7 files changed, 592 insertions(+), 462 deletions(-) diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index f1ed4e9dc..5ca5f187c 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -209,7 +209,7 @@ subroutine CPFEM_init !$OMP CRITICAL (write2out) write(6,*) - write(6,*) '<<<+- cpfem init -+>>>' + write(6,*) '<<<+- CPFEM init -+>>>' write(6,*) '$Id$' #include "compilation_info.f90" if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 4bf413311..cd5346251 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -238,15 +238,9 @@ program DAMASK_spectral complex(pReal), dimension(:,:,:), pointer :: scalarField_real complex(pReal), dimension(:,:,:), pointer :: scalarField_fourier integer(pInt) :: row, column - !################################################################################################## ! reading of information from load case file and geometry file !################################################################################################## -#ifdef PETSC - integer :: ierr_psc - call PetscInitialize(PETSC_NULL_CHARACTER, ierr_psc) -#endif - !-------------------------------------------------------------------------------------------------- ! initialization of all related DAMASK modules (e.g. mesh.f90 reads in geometry) call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt) @@ -393,13 +387,13 @@ program DAMASK_spectral else write(6,'(a)')'deformation gradient rate:' endif - write (6,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& + write(6,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(loadcase)%deformation),& reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskDeformation)) - write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',& + write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',& 1e-9_pReal*merge(math_transpose33(bc(loadcase)%stress),& reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(loadcase)%maskStress)) if (any(bc(loadcase)%rotation /= math_I3)) & - write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& + write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& math_transpose33(bc(loadcase)%rotation) write(6,'(a,f12.6)') 'temperature:', bc(loadcase)%temperature write(6,'(a,f12.6)') 'time: ', bc(loadcase)%time @@ -539,7 +533,7 @@ program DAMASK_spectral enddo; enddo; enddo if (debugGeneral) write(6,'(a)') 'First call to CPFEM finished' C = C * wgt - + !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) @@ -805,7 +799,7 @@ program DAMASK_spectral P_av_lab = real(P_fourier(1,1,1,1:3,1:3),pReal)*wgt P_av = math_rotate_forward33(P_av_lab,bc(loadcase)%rotation) - write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& + write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& math_transpose33(P_av)/1.e6_pReal !-------------------------------------------------------------------------------------------------- @@ -1051,7 +1045,7 @@ program DAMASK_spectral if (mod(inc,bc(loadcase)%outputFrequency) == 0_pInt) then ! at output frequency write(6,'(a)') '' write(6,'(a)') '... writing results to file ......................................' - write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! write result to file + write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! write result to file flush(538) endif diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index d0e245f60..6992bb5f8 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -1,5 +1,5 @@ !-------------------------------------------------------------------------------------------------- -!* $Id$ +! $Id$ !-------------------------------------------------------------------------------------------------- !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH @@ -15,12 +15,10 @@ program DAMASK_spectral_Driver getSolverWorkingDirectoryName, & getSolverJobName, & appendToOutFile - use prec, only: & pInt, & pReal, & DAMASK_NaN - use IO, only: & IO_isBlank, & IO_open_file, & @@ -33,36 +31,28 @@ program DAMASK_spectral_Driver IO_read_jobBinaryFile, & IO_write_jobBinaryFile, & IO_intOut - - use math - + use math ! need to include the whole module for FFTW use mesh, only : & res, & geomdim, & mesh_NcpElems - use CPFEM, only: & CPFEM_initAll - use FEsolving, only: & restartWrite, & restartInc - use numerics, only: & maxCutBack, & rotation_tol, & mySpectralSolver - use homogenization, only: & materialpoint_sizeResults, & materialpoint_results - use DAMASK_spectral_Utilities, only: & tBoundaryCondition, & tSolutionState, & debugGeneral, & cutBack - use DAMASK_spectral_SolverBasic #ifdef PETSc use DAMASK_spectral_SolverBasicPETSC @@ -96,12 +86,11 @@ program DAMASK_spectral_Driver N_l = 0_pInt, & N_t = 0_pInt, & N_n = 0_pInt, & - N_Fdot = 0_pInt, & ! number of Fourier points + N_Fdot = 0_pInt, & ! myUnit = 234_pInt character(len=1024) :: & line - !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeros = 0.0_pReal @@ -120,11 +109,11 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt) - write(6,'(a)') '' write(6,'(a)') ' <<<+- DAMASK_spectral_driver init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" + !-------------------------------------------------------------------------------------------------- ! reading basic information from load case file and allocate data structure containing load cases call IO_open_file(myUnit,trim(loadCaseFile)) @@ -147,12 +136,12 @@ program DAMASK_spectral_Driver enddo ! count all identifiers to allocate memory and do sanity check enddo -100 if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check - call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase +100 if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check + call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase allocate (loadCases(N_n)) loadCases%P%myType='p' - !-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure rewind(myUnit) do @@ -163,7 +152,7 @@ program DAMASK_spectral_Driver do i = 1_pInt,maxNchunksLoadcase select case (IO_lc(IO_stringValue(line,positions,i))) case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient') ! assign values for the deformation BC matrix - if (IO_lc(IO_stringValue(line,positions,i)) == 'l'.or. & ! in case of given L, set flag to true + if (IO_lc(IO_stringValue(line,positions,i)) == 'l'.or. & ! in case of given L, set flag to true IO_lc(IO_stringValue(line,positions,i)) == 'velocitygrad'.or.& IO_lc(IO_stringValue(line,positions,i)) == 'velgrad'.or.& IO_lc(IO_stringValue(line,positions,i)) == 'velocitygradient') then @@ -203,7 +192,7 @@ program DAMASK_spectral_Driver case('r','restart','restartwrite') ! frequency of writing restart information loadCases(currentLoadCase)%restartfrequency = max(0_pInt,IO_intValue(line,positions,i+1_pInt)) case('guessreset','dropguessing') - loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory + loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory case('euler') ! rotation of currentLoadCase given in euler angles l = 0_pInt ! assuming values given in degrees k = 0_pInt ! assuming keyword indicating degree/radians @@ -228,58 +217,65 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case - loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase + loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase errorID = 0_pInt checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases) write (loadcase_string, '(i6)' ) currentLoadCase write(6,'(1x,a,i6)') 'load case: ', currentLoadCase - if (.not. loadCases(currentLoadCase)%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory' + if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & + write(6,'(2x,a)') 'drop guessing along trajectory' if (loadCases(currentLoadCase)%deformation%myType=='l') then do j = 1_pInt, 3_pInt if (any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & - any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined + any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) & + errorID = 832_pInt ! each row should be either fully or not at all defined enddo write(6,'(2x,a)') 'velocity gradient:' else write(6,'(2x,a)') 'deformation gradient rate:' endif - write(6,'(3(3(3x,f12.7,1x)/))',advance='no') merge(math_transpose33(loadCases(currentLoadCase)%deformation%values),& - reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%deformation%maskLogical)) + write(6,'(3(3(3x,f12.7,1x)/))',advance='no') & + merge(math_transpose33(loadCases(currentLoadCase)%deformation%values), & + reshape(spread(huge(1.0_pReal),1,9),[ 3,3]), & + transpose(loadCases(currentLoadCase)%deformation%maskLogical)) + if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. & + loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only + if (any(loadCases(currentLoadCase)%P%maskLogical .and. & + transpose(loadCases(currentLoadCase)%P%maskLogical) .and. & + reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & + errorID = 838_pInt ! no rotation is allowed by stress BC write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'stress / GPa:',& - 1e-9_pReal*merge(math_transpose33(loadCases(currentLoadCase)%P%values),& - reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%P%maskLogical)) + 1e-9_pReal*merge(math_transpose33(loadCases(currentLoadCase)%P%values),& + reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),& + transpose(loadCases(currentLoadCase)%P%maskLogical)) + if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation, & + math_transpose33(loadCases(currentLoadCase)%rotation))-math_I3) >& + reshape(spread(rotation_tol,1,9),[ 3,3]))& + .or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > & + 1.0_pReal + rotation_tol) errorID = 846_pInt ! given rotation matrix contains strain if (any(loadCases(currentLoadCase)%rotation /= math_I3)) & write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& - math_transpose33(loadCases(currentLoadCase)%rotation) + math_transpose33(loadCases(currentLoadCase)%rotation) write(6,'(2x,a,f12.6)') 'temperature:', loadCases(currentLoadCase)%temperature + if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time + if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs - write(6,'(2x,a,i5)') 'output frequency: ', loadCases(currentLoadCase)%outputfrequency - write(6,'(2x,a,i5,/)') 'restart frequency: ', loadCases(currentLoadCase)%restartfrequency - if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only - if (any(loadCases(currentLoadCase)%P%maskLogical .and. transpose(loadCases(currentLoadCase)%P%maskLogical) .and. & - reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & - errorID = 838_pInt ! no rotation is allowed by stress BC - if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation,math_transpose33(loadCases(currentLoadCase)%rotation))& - -math_I3) > reshape(spread(rotation_tol,1,9),[ 3,3]))& - .or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > 1.0_pReal + rotation_tol)& - errorID = 846_pInt ! given rotation matrix contains strain - if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment - if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count - if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency + if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency + write(6,'(2x,a,i5)') 'output frequency: ', & + loadCases(currentLoadCase)%outputfrequency + write(6,'(2x,a,i5,/)') 'restart frequency: ', & + loadCases(currentLoadCase)%restartfrequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) enddo checkLoadcases select case (myspectralsolver) - case (DAMASK_spectral_SolverBasic_label) call basic_init() - #ifdef PETSc case (DAMASK_spectral_SolverBasicPETSc_label) call basicPETSc_init() - case (DAMASK_spectral_SolverAL_label) call AL_init() #endif @@ -291,14 +287,12 @@ program DAMASK_spectral_Driver ! write header of output file if (appendToOutFile) then open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD') + '.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD') open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.sta',form='FORMATTED', position='APPEND', status='OLD') + '.sta',form='FORMATTED', position='APPEND', status='OLD') else open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.spectralOut',form='UNFORMATTED',status='REPLACE') - open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& - '.sta',form='FORMATTED',status='REPLACE') + '.spectralOut',form='UNFORMATTED',status='REPLACE') write(resUnit) 'load', trim(loadCaseFile) write(resUnit) 'workingdir', trim(getSolverWorkingDirectoryName()) write(resUnit) 'geometry', trim(geometryFile) @@ -306,16 +300,18 @@ program DAMASK_spectral_Driver write(resUnit) 'dimension', geomdim write(resUnit) 'materialpoint_sizeResults', materialpoint_sizeResults write(resUnit) 'loadcases', size(loadCases) - write(resUnit) 'frequencies', loadCases%outputfrequency ! one entry per currentLoadCase - write(resUnit) 'times', loadCases%time ! one entry per currentLoadCase + write(resUnit) 'frequencies', loadCases%outputfrequency ! one entry per currentLoadCase + write(resUnit) 'times', loadCases%time ! one entry per currentLoadCase write(resUnit) 'logscales', loadCases%logscale - write(resUnit) 'increments', loadCases%incs ! one entry per currentLoadCase - write(resUnit) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc - write(resUnit) 'eoh' ! end of header - write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results + write(resUnit) 'increments', loadCases%incs ! one entry per currentLoadCase + write(resUnit) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc + write(resUnit) 'eoh' ! end of header + write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results + open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + '.sta',form='FORMATTED',status='REPLACE') + write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' if (debugGeneral) write(6,'(a)') 'Header of result file written out' endif - !-------------------------------------------------------------------------------------------------- ! loopping over loadcases loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) @@ -344,10 +340,10 @@ program DAMASK_spectral_Driver timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) endif else ! not-1st currentLoadCase of logarithmic scale - timeinc = time0 *( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal))& - -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& - real(loadCases(currentLoadCase)%incs ,pReal)) ) + timeinc = time0 *( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& + real(loadCases(currentLoadCase)%incs ,pReal))& + -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& + real(loadCases(currentLoadCase)%incs ,pReal)) ) endif endif timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) @@ -364,12 +360,12 @@ program DAMASK_spectral_Driver ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//& ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & - 'Time', time, & - 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& - '-', stepFraction, '/', 2_pInt**cutBackLevel,& - ' of load case ', currentLoadCase,'/',size(loadCases) + 'Time', time, & + 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& + '-', stepFraction, '/', 2_pInt**cutBackLevel,& + ' of load case ', currentLoadCase,'/',size(loadCases) write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases(:)%incs))//& - ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//')') & + ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//')') & 'Inc. ',totalIncsCounter,'/',sum(loadCases(:)%incs),& '-',stepFraction, '/', 2_pInt**cutBackLevel select case(myspectralsolver) @@ -435,8 +431,13 @@ program DAMASK_spectral_Driver write(6,'(1/,a)') '... writing results to file ......................................' write(resUnit) materialpoint_results ! write result to file endif + if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & + mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! at frequency of writing restart information set restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) + restartWrite = .true. + endif else !just time forwarding - time = time + timeinc + time = time + timeinc + guessmode = 1.0_pReal endif ! end calculation/forwarding enddo incLooping diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 65bc705d1..253321784 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -208,7 +208,7 @@ subroutine AL_init() close (777) endif - call Utilities_updateGamma(C) + call Utilities_updateGamma(C,.True.) C_scale = C S_scale = math_invSym3333(C) @@ -328,7 +328,7 @@ else !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C) + if (update_gamma) call Utilities_updateGamma(C,restartWrite) ForwardData = .True. mask_stress = P_BC%maskFloat diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index 4555425ff..ab33553dd 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -33,7 +33,8 @@ module DAMASK_spectral_SolverBasic ! stress, stiffness and compliance average etc. real(pReal), private, dimension(3,3) :: & F_aim = math_I3, & - F_aim_lastInc = math_I3 + F_aim_lastInc = math_I3, & + F_aimDot = 0.0_pReal real(pReal), private,dimension(3,3,3,3) :: & C = 0.0_pReal, C_lastInc = 0.0_pReal @@ -69,8 +70,9 @@ subroutine basic_init() implicit none integer(pInt) :: i,j,k real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3,3,3) :: temp3333_Real - + call Utilities_Init() write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(a)') ' $Id$' @@ -81,8 +83,8 @@ subroutine basic_init() allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (P ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) - allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) - allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) + allocate (coordinates( res(1), res(2),res(3),3),source = 0.0_pReal) + allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! init fields @@ -93,6 +95,7 @@ subroutine basic_init() coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - geomdim/real(2_pInt*res,pReal) enddo; enddo; enddo + elseif (restartInc > 1_pInt) then ! using old values from file if (debugRestart) write(6,'(a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'Reading values of increment', restartInc - 1_pInt, 'from file' @@ -110,26 +113,32 @@ subroutine basic_init() call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) read (777,rec=1) F_aim_lastInc close (777) - + call IO_read_jobBinaryFile(777,'C_lastInc',trim(getSolverJobName()),size(C_lastInc)) + read (777,rec=1) C_lastInc + close (777) + call IO_read_jobBinaryFile(777,'C',trim(getSolverJobName()),size(C)) + read (777,rec=1) C + close (777) + call IO_read_jobBinaryFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) + read (777,rec=1) f_aimDot + close (777) + call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real)) + read (777,rec=1) temp3333_Real + close (777) coordinates = 0.0 ! change it later!!! endif - !no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) - call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& + call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,& P,C,temp33_Real,.false.,math_I3) + !no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) + !-------------------------------------------------------------------------------------------------- ! reference stiffness - if (restartInc == 1_pInt) then - call IO_write_jobBinaryFile(777,'C_ref',size(C)) - write (777,rec=1) C - close(777) - elseif (restartInc > 1_pInt) then - call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) - read (777,rec=1) C - close (777) - endif + if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness + temp3333_Real = C + endif - call Utilities_updateGamma(C) + call Utilities_updateGamma(temp3333_Real,.True.) end subroutine basic_init @@ -173,8 +182,11 @@ type(tSolutionState) function & use FEsolving, only: & restartWrite, & + restartRead, & terminallyIll - use DAMASK_spectral_Utilities, only: cutBack + use DAMASK_spectral_Utilities, only: & + cutBack + implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution @@ -184,8 +196,7 @@ type(tSolutionState) function & real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3,3,3) :: S - real(pReal), dimension(3,3), save :: f_aimDot = 0.0_pReal - real(pReal), dimension(3,3) :: F_aim_lab, & + real(pReal), dimension(3,3) :: F_aim_lab, & F_aim_lab_lastIter, & P_av !-------------------------------------------------------------------------------------------------- @@ -200,21 +211,35 @@ type(tSolutionState) function & ! restart information for spectral solver if (restartWrite) then write(6,'(a)') 'writing converged results for restart' - call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) - write (777,rec=1) F_LastInc + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file + write (777,rec=1) F close (777) + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',size(F_lastInc)) ! writing F_lastInc field to file + write (777,rec=1) F_lastInc + close (777) + call IO_write_jobBinaryFile(777,'F_aim',size(F_aim)) + write (777,rec=1) F_aim + close(777) + call IO_write_jobBinaryFile(777,'F_aim_lastInc',size(F_aim_lastInc)) + write (777,rec=1) F_aim_lastInc + close(777) call IO_write_jobBinaryFile(777,'C',size(C)) write (777,rec=1) C close(777) + call IO_write_jobBinaryFile(777,'C_lastInc',size(C_lastInc)) + write (777,rec=1) C_lastInc + close(777) + call IO_write_jobBinaryFile(777,'F_aimDot',size(f_aimDot)) + write (777,rec=1) f_aimDot + close(777) endif - if ( cutBack) then - F_aim = F_aim_lastInc - F = F_lastInc - C = C_lastInc -else - C_lastInc = C + F_aim = F_aim_lastInc + F = F_lastInc + C = C_lastInc + else + C_lastInc = C !-------------------------------------------------------------------------------------------------- ! calculate rate for aim if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F @@ -225,7 +250,8 @@ else f_aimDot = f_aimDot & + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim - + print*, 'F_aimDot', f_aimDot + print*, 'guessmode', guessmode !-------------------------------------------------------------------------------------------------- ! update coordinates and rate and forward last inc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & @@ -235,15 +261,15 @@ else F_lastInc = F endif F_aim = F_aim + f_aimDot * timeinc - F = Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot) !I think F aim should be rotated here + F = Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot) !I think F aim should be rotated here !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C) + if (update_gamma) call Utilities_updateGamma(C,restartWrite) - ForwardData = .True. + if (.not. restartRead) ForwardData = .True. iter = 0_pInt convergenceLoop: do while(iter < itmax) diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 1b297cd0c..a01298c3a 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -191,7 +191,7 @@ subroutine BasicPETSC_init() close (777) endif - call Utilities_updateGamma(C) + call Utilities_updateGamma(C,.True.) end subroutine BasicPETSC_init @@ -294,7 +294,7 @@ else !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C) + if (update_gamma) call Utilities_updateGamma(C,restartWrite) ForwardData = .True. mask_stress = P_BC%maskFloat diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 06ba4b5f7..a0648915e 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -1,70 +1,81 @@ !-------------------------------------------------------------------------------------------------- -!* $Id$ +! $Id$ !-------------------------------------------------------------------------------------------------- !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Utilities used by the different spectral solver variants !-------------------------------------------------------------------------------------------------- -module DAMASK_spectral_Utilities +module DAMASK_spectral_utilities + use, intrinsic :: iso_c_binding use prec, only: & pReal, & pInt - use mesh, only : & - res, & - res1_red, & - geomdim, & - mesh_NcpElems, & - wgt - - use math - - use IO, only: & - IO_error - + implicit none + logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill + !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - type(C_PTR), private :: plan_forward, plan_backward ! plans for fftw - real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method - real(pReal), private, dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator - complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier - real(pReal), private, dimension(3,3,3,3) :: C_ref - - real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real - logical, public :: cutBack =.false. + real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress of deformation) of field_fourier + type(C_PTR), private :: plan_forward, plan_backward !< plans for FFTW + real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method + real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator + complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates + real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness + !-------------------------------------------------------------------------------------------------- ! debug fftw - type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back - complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real - complex(pReal),private, dimension(:,:,:), pointer :: scalarField_fourier + type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back !< plans for FFTW in case of debugging the Fourier transform + complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for debug of FFTW + complex(pReal),private, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field complex representation for debug of FFTW !-------------------------------------------------------------------------------------------------- ! debug divergence - type(C_PTR), private :: plan_divergence - real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real - complex(pReal), private, dimension(:,:,:,:), pointer :: divergence_fourier - real(pReal), private, dimension(:,:,:,:), allocatable :: divergence_post + type(C_PTR), private :: plan_divergence !< plan for FFTW in case of debugging divergence calculation + real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real !< scalar field real representation for debugging divergence calculation + complex(pReal),private, dimension(:,:,:,:), pointer :: divergence_fourier !< scalar field real representation for debugging divergence calculation + real(pReal), private, dimension(:,:,:,:), allocatable :: divergence_post !< data of divergence calculation using function from core modules (serves as a reference) !-------------------------------------------------------------------------------------------------- -!variables controlling debugging - logical,public :: debugGeneral, debugDivergence, debugRestart, debugFFTW +! variables controlling debugging + logical, public :: & + debugGeneral, & !< general debugging of spectral solver + debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing) + debugRestart, & !< debbuging of restart features + debugFFTW !< doing additional FFT on scalar field and compare to results of strided 3D FFT !-------------------------------------------------------------------------------------------------- ! derived types - type tSolutionState - logical :: converged = .true. - logical :: regrid = .false. - logical :: termIll = .false. - integer(pInt) :: iterationsNeeded = 0_pInt + type tSolutionState !< return type of solution from spectral solver variants + logical :: converged = .true. + logical :: regrid = .false. + logical :: termIll = .false. + integer(pInt) :: iterationsNeeded = 0_pInt end type tSolutionState - type tBoundaryCondition - real(pReal), dimension(3,3) :: values = 0.0_pReal - real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal + type tBoundaryCondition !< set of parameters defining a boundary condition + real(pReal), dimension(3,3) :: values = 0.0_pReal + real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal logical, dimension(3,3) :: maskLogical = .false. - character(len=64) :: myType = 'None' + character(len=64) :: myType = 'None' end type tBoundaryCondition + + public :: & + utilities_init, & + utilities_updateGamma, & + utilities_FFTforward, & + utilities_FFTbackward, & + utilities_fourierConvolution, & + utilities_divergenceRMS, & + utilities_maskedCompliance, & + utilities_constitutiveResponse, & + utilities_calculateRate, & + utilities_forwardField, & + utilities_destroy + + private :: & + utilities_getFilter contains @@ -76,14 +87,15 @@ contains !> level chosen. !> Initializes FFTW. !-------------------------------------------------------------------------------------------------- -subroutine Utilities_init() +subroutine utilities_init() use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use IO, only: & + IO_error use numerics, only: & DAMASK_NumThreadsInt, & fftw_planner_flag, & fftw_timelimit, & memory_efficient - use debug, only: & debug_level, & debug_spectral, & @@ -91,24 +103,25 @@ subroutine Utilities_init() debug_spectralDivergence, & debug_spectralRestart, & debug_spectralFFTW - - use mesh, only : & + use mesh, only: & + res, & + res1_red, & virt_dim + use math ! must use the whole module for use of FFTW implicit none - integer(pInt) :: i, j, k + integer(pInt) :: i, j, k integer(pInt), dimension(3) :: k_s - !$ integer(pInt) :: ierr - type(C_PTR) :: tensorField ! field in real and fourier space - type(C_PTR) :: scalarField_realC, scalarField_fourierC - type(C_PTR) :: divergence + type(C_PTR) :: & + tensorField, & !< field cotaining data for FFTW in real and fourier space (in place) + scalarField_realC, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place) + scalarField_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place) + divergence !< field cotaining data for FFTW in real and fourier space when debugging divergence (in place) - - write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectral_utilities init -+>>>' - write(6,'(a)') ' $Id$' + write(6,'(/,a)') ' <<<+- DAMASK_spectral_utilities init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" - write(6,'(a)') '' + write(6,'(a)') '' !-------------------------------------------------------------------------------------------------- ! set debugging parameters @@ -119,34 +132,34 @@ subroutine Utilities_init() !-------------------------------------------------------------------------------------------------- ! allocation - allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) ! start out isothermally - tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) - call c_f_pointer(tensorField, field_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField - call c_f_pointer(tensorField, field_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField + allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) ! start out isothermally + tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) + call c_f_pointer(tensorField, field_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField + call c_f_pointer(tensorField, field_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField !-------------------------------------------------------------------------------------------------- -! general initialization of fftw (see manual on fftw.org for more details) +! general initialization of FFTW (see manual on fftw.org for more details) if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=808_pInt) ! check for correct precision in C !$ if(DAMASK_NumThreadsInt > 0_pInt) then -!$ ierr = fftw_init_threads() -!$ if (ierr == 0_pInt) call IO_error(error_ID = 809_pInt) +!$ i = fftw_init_threads() ! returns 0 in case of problem +!$ if (i == 0_pInt) call IO_error(error_ID = 809_pInt) !$ call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) !$ endif call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation !-------------------------------------------------------------------------------------------------- ! creating plans - plan_forward = fftw_plan_many_dft_r2c(3,[ res(3),res(2) ,res(1)],9,& ! dimensions , length in each dimension in reversed order - field_real,[ res(3),res(2) ,res(1)+2_pInt],& ! input data , physical length in each dimension in reversed order - 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions - field_fourier,[ res(3),res(2) ,res1_red],& - 1, res(3)*res(2)* res1_red,fftw_planner_flag) + plan_forward = fftw_plan_many_dft_r2c(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms + field_real, [res(3),res(2) ,res(1)+2_pInt],& ! input data, physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding, product of physical length in the 3 dimensions + field_fourier, [res(3),res(2) ,res1_red],& ! output data, physical length in each dimension in reversed order + 1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision - plan_backward =fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],9,& - field_fourier,[ res(3),res(2) ,res1_red],& - 1, res(3)*res(2)* res1_red,& - field_real,[ res(3),res(2) ,res(1)+2_pInt],& - 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + plan_backward = fftw_plan_many_dft_c2r(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms + field_fourier, [res(3),res(2) ,res1_red],& ! input data, physical length in each dimension in reversed order + 1, res(3)*res(2)* res1_red,& ! striding, product of physical length in the 3 dimensions + field_real, [res(3),res(2) ,res(1)+2_pInt],& ! output data, physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt), fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision !-------------------------------------------------------------------------------------------------- ! depending on (debug) options, allocate more memory and create additional plans @@ -154,7 +167,7 @@ subroutine Utilities_init() divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T)) call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3]) call c_f_pointer(divergence, divergence_fourier, [ res1_red, res(2),res(3),3]) - allocate (divergence_post(res(1),res(2),res(3),3)); divergence_post = 0.0_pReal + allocate (divergence_post(res(1),res(2),res(3),3),source = 0.0_pReal) plan_divergence = fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],3,& divergence_fourier,[ res(3),res(2) ,res1_red],& 1, res(3)*res(2)* res1_red,& @@ -163,14 +176,14 @@ subroutine Utilities_init() endif if (debugFFTW) then - scalarField_realC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! do not do an inplace transform - scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) - call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) - call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) - plan_scalarField_forth = fftw_plan_dft_3d(res(3),res(2),res(1),& !reversed order - scalarField_real,scalarField_fourier,-1,fftw_planner_flag) - plan_scalarField_back = fftw_plan_dft_3d(res(3),res(2),res(1),& !reversed order - scalarField_fourier,scalarField_real,+1,fftw_planner_flag) + scalarField_realC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for real representation (no in place transform) + scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for fourier representation (no in place transform) + call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) ! place a pointer for a real representation + call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) ! place a pointer for a fourier representation + plan_scalarField_forth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) + scalarField_real,scalarField_fourier,-1,fftw_planner_flag) ! input, output, forward FFT(-1), planner precision + plan_scalarField_back = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) + scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision endif if (debugGeneral) write(6,'(a)') 'FFTW initialized' @@ -179,58 +192,78 @@ subroutine Utilities_init() ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) do k = 1_pInt, res(3) k_s(3) = k - 1_pInt - if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) + if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do j = 1_pInt, res(2) k_s(2) = j - 1_pInt - if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) + if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1 do i = 1_pInt, res1_red - k_s(1) = i - 1_pInt - xi(1:3,i,j,k) = real(k_s, pReal)/virt_dim + k_s(1) = i - 1_pInt ! symmetry, junst running from 0,1,...,N/2,N/2+1 + xi(1:3,i,j,k) = real(k_s, pReal)/virt_dim ! if divergence_correction is set, frequencies are calculated on unit length enddo; enddo; enddo if(memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(3,3,3,3,1,1,1), source = 0.0_pReal) else ! precalculation of gamma_hat field - allocate (gamma_hat(3,3,3,3,res1_red ,res(2),res(3)), source =0.0_pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + allocate (gamma_hat(3,3,3,3,res1_red ,res(2),res(3)), source =0.0_pReal) endif -end subroutine Utilities_init +end subroutine utilities_init + !-------------------------------------------------------------------------------------------------- !> @brief updates references stiffness and potentially precalculated gamma operator !> @details Sets the current reference stiffness to the stiffness given as an argument. !> If the gamma operator is precalculated, it is calculated with this stiffness. !> In case of a on-the-fly calculation, only the reference stiffness is updated. -!> The gamma operator is filtered depening on the filter selected in numerics +!> The gamma operator is filtered depening on the filter selected in numerics. +!> Also writes out the current reference stiffness for restart. !-------------------------------------------------------------------------------------------------- -subroutine Utilities_updateGamma(C) - +subroutine utilities_updateGamma(C,saveReference) + use IO, only: & + IO_write_jobBinaryFile use numerics, only: & memory_efficient - + use math, only: & + math_inv33 + use mesh, only: & + res, & + res1_red + implicit none - real(pReal), dimension(3,3,3,3), intent(in) :: C - real(pReal), dimension(3,3) :: temp33_Real, xiDyad - real(pReal) :: filter - integer(pInt) :: i, j, k, l, m, n, o + real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness + logical , intent(in) :: saveReference !< save reference stiffness to file for restart + real(pReal), dimension(3,3) :: temp33_Real, xiDyad + real(pReal) :: filter !< weighting of current component + integer(pInt) :: & + i, j, k, & + l, m, n, o C_ref = C + if (saveReference) then + write(6,'(a)') 'writing reference stiffness to file' + call IO_write_jobBinaryFile(777,'C_ref',size(C_ref)) + write (777,rec=1) C_ref + close(777) + endif + if(.not. memory_efficient) then do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red - if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & temp33_Real(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) temp33_Real = math_inv33(temp33_Real) - filter = Utilities_getFilter(xi(1:3,i,j,k)) + filter = utilities_getFilter(xi(1:3,i,j,k)) ! weighting factor computed by getFilter function forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& gamma_hat(l,m,n,o, i,j,k) = filter*temp33_Real(l,n)*xiDyad(m,o) endif enddo; enddo; enddo - gamma_hat(1:3,1:3,1:3,1:3, 1,1,1) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + gamma_hat(1:3,1:3,1:3,1:3, 1,1,1) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 endif -end subroutine Utilities_updateGamma + +end subroutine utilities_updateGamma + !-------------------------------------------------------------------------------------------------- !> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed @@ -238,30 +271,31 @@ end subroutine Utilities_updateGamma !> In case of debugging the FFT, also one component of the tensor (specified by row and column) !> is independetly transformed complex to complex and compared to the whole tensor transform !-------------------------------------------------------------------------------------------------- -subroutine Utilities_FFTforward(row,column) +subroutine utilities_FFTforward(row,column) + use math use mesh, only : & - virt_dim - use math, only: & - math_divergenceFFT - - implicit none - integer(pInt), intent(in), optional :: row, column + virt_dim, & + res, & + res1_red + + implicit none + integer(pInt), intent(in), optional :: row, column !< if debug FFTW, compare 3D array field of row and column !-------------------------------------------------------------------------------------------------- ! copy one component of the stress field to to a single FT and check for mismatch if (debugFFTW) then if (.not. present(row) .or. .not. present(column)) stop - scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component + scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) endif !-------------------------------------------------------------------------------------------------- ! call function to calculate divergence from math (for post processing) to check results if (debugDivergence) & - divergence_post = math_divergenceFFT(virt_dim,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3)) ! padding + divergence_post = math_divergenceFFT(virt_dim,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3)) ! some elements are padded !-------------------------------------------------------------------------------------------------- -! doing the FT +! doing the FFT call fftw_execute_dft_r2c(plan_forward,field_real,field_fourier) !-------------------------------------------------------------------------------------------------- @@ -269,7 +303,7 @@ subroutine Utilities_FFTforward(row,column) if (debugFFTW) then call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) write(6,'(a,i1,1x,i1)') 'checking FT results of compontent ', row, column - write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',& + write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',& ! print real and imaginary part seperately maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & @@ -277,16 +311,17 @@ subroutine Utilities_FFTforward(row,column) field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) endif + !-------------------------------------------------------------------------------------------------- ! removing highest frequencies field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) field_fourier (1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(3)>1_pInt) & + if(res(3)>1_pInt) & ! do not delete the whole slice in case of 2D calculation field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) -end subroutine Utilities_FFTforward +end subroutine utilities_FFTforward !-------------------------------------------------------------------------------------------------- @@ -296,32 +331,38 @@ end subroutine Utilities_FFTforward !> is independetly transformed complex to complex and compared to the whole tensor transform !> results is weighted by number of points stored in wgt !-------------------------------------------------------------------------------------------------- -subroutine Utilities_FFTbackward(row,column) +subroutine utilities_FFTbackward(row,column) + use math !< must use the whole module for use of FFTW + use mesh, only: & + wgt, & + res, & + res1_red - implicit none - integer(pInt), intent(in), optional :: row, column - integer(pInt) :: i, j, k, m, n + implicit none + integer(pInt), intent(in), optional :: row, column !< if debug FFTW, compare 3D array field of row and column + integer(pInt) :: i, j, k, m, n !-------------------------------------------------------------------------------------------------- -! comparing 1 and 3x3 inverse FT results - if (debugFFTW) then - scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) - do i = 0_pInt, res(1)/2_pInt-2_pInt ! unpack fft data for conj complex symmetric part - m = 1_pInt - do k = 1_pInt, res(3) - n = 1_pInt - do j = 1_pInt, res(2) - scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) - if(n == 1_pInt) n = res(2) + 1_pInt - n = n-1_pInt - enddo - if(m == 1_pInt) m = res(3) + 1_pInt - m = m -1_pInt - enddo; enddo - endif - - - call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient +! unpack FFT data for conj complex symmetric part + if (debugFFTW) then + scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) + do i = 0_pInt, res(1)/2_pInt-2_pInt + m = 1_pInt + do k = 1_pInt, res(3) + n = 1_pInt + do j = 1_pInt, res(2) + scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) + if(n == 1_pInt) n = res(2) + 1_pInt + n = n-1_pInt + enddo + if(m == 1_pInt) m = res(3) + 1_pInt + m = m -1_pInt + enddo; enddo + endif + +!-------------------------------------------------------------------------------------------------- +! doing the iFFT + call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results @@ -333,40 +374,47 @@ subroutine Utilities_FFTbackward(row,column) field_real(1:res(1),1:res(2),1:res(3),row,column))/& real(scalarField_real(1:res(1),1:res(2),1:res(3)))) endif - field_real = field_real * wgt + + field_real = field_real * wgt ! normalize the result by number of elements -end subroutine Utilities_FFTbackward +end subroutine utilities_FFTbackward !-------------------------------------------------------------------------------------------------- -!> @brief doing convolution gamma_hat * field_real with average value given by fieldAim +!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !-------------------------------------------------------------------------------------------------- -subroutine Utilities_fourierConvolution(fieldAim) - +subroutine utilities_fourierConvolution(fieldAim) use numerics, only: & memory_efficient - - implicit none - real(pReal), dimension(3,3), intent(in) :: fieldAim - real(pReal), dimension(3,3) :: xiDyad, temp33_Real - real(pReal) :: filter - integer(pInt) :: i, j, k, l, m, n, o - complex(pReal), dimension(3,3) :: temp33_complex + use math, only: & + math_inv33 + use mesh, only: & + mesh_NcpElems, & + res, & + res1_red + implicit none + real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution + real(pReal), dimension(3,3) :: xiDyad, temp33_Real + real(pReal) :: filter !< weighting of current component + complex(pReal), dimension(3,3) :: temp33_complex + integer(pInt) :: & + i, j, k, & + l, m, n, o write(6,'(/,a)') '... doing convolution .....................................................' !-------------------------------------------------------------------------------------------------- ! to the actual spectral method calculation (mechanical equilibrium) - if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat + if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red - if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & temp33_Real(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) temp33_Real = math_inv33(temp33_Real) - filter = Utilities_getFilter(xi(1:3,i,j,k)) + filter = utilities_getFilter(xi(1:3,i,j,k)) ! weighting factor computed by getFilter function forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& gamma_hat(l,m,n,o, 1,1,1) = filter*temp33_Real(l,n)*xiDyad(m,o) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & @@ -374,99 +422,118 @@ subroutine Utilities_fourierConvolution(fieldAim) field_fourier(i,j,k,1:3,1:3) = temp33_Complex endif enddo; enddo; enddo - else ! use precalculated gamma-operator + else ! use precalculated gamma-operator do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) & temp33_Complex(m,n) = sum(gamma_hat(m,n,1:3,1:3, i,j,k) * field_fourier(i,j,k,1:3,1:3)) field_fourier(i,j,k, 1:3,1:3) = temp33_Complex enddo; enddo; enddo endif + field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(mesh_NcpElems,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(mesh_NcpElems,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - -end subroutine Utilities_fourierConvolution +end subroutine utilities_fourierConvolution !-------------------------------------------------------------------------------------------------- !> @brief calculate root mean square of divergence of field_fourier !-------------------------------------------------------------------------------------------------- -real(pReal) function Utilities_divergenceRMS() - implicit none - integer(pInt) :: i, j, k - real(pReal) :: err_div_RMS, err_real_div_RMS, err_post_div_RMS,& - err_div_max, err_real_div_max - complex(pReal), dimension(3) :: temp3_complex +real(pReal) function utilities_divergenceRMS() + use math !< must use the whole module for use of FFTW + use mesh, only: & + wgt, & + res, & + res1_red - write(6,'(/,a)') '... calculating divergence ................................................' + implicit none + integer(pInt) :: i, j, k + real(pReal) :: & + err_div_RMS, & !< RMS of divergence in Fourier space + err_real_div_RMS, & !< RMS of divergence in real space + err_post_div_RMS, & !< RMS of divergence in Fourier space, calculated using function for post processing + err_div_max, & !< maximum value of divergence in Fourier space + err_real_div_max !< maximum value of divergence in real space + complex(pReal), dimension(3) :: temp3_complex + + write(6,'(/,a)') '... calculating divergence ................................................' !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space - Utilities_divergenceRMS = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2) - do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. - Utilities_divergenceRMS = Utilities_divergenceRMS & - + 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again - xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector - +sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& - xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) - enddo - Utilities_divergenceRMS = Utilities_divergenceRMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart - + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) - enddo; enddo + utilities_divergenceRMS = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2) + do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. + utilities_divergenceRMS = utilities_divergenceRMS & + + 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector + +sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& + xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) + enddo + utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart + + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + + sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& + xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& + xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) + enddo; enddo - Utilities_divergenceRMS = sqrt(Utilities_divergenceRMS) *wgt ! RMS in real space calculated with Parsevals theorem from Fourier space - + utilities_divergenceRMS = sqrt(utilities_divergenceRMS) *wgt ! RMS in real space calculated with Parsevals theorem from Fourier space + !-------------------------------------------------------------------------------------------------- ! calculate additional divergence criteria and report - if (debugDivergence) then ! calculate divergence again - err_div_max = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red - temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier - xi(1:3,i,j,k))*TWOPIIMG - err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) - divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared - enddo; enddo; enddo - - call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted + if (debugDivergence) then ! calculate divergence again + err_div_max = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red + temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier + xi(1:3,i,j,k))*TWOPIIMG + err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) + divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared + enddo; enddo; enddo + + call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted - err_real_div_RMS = sqrt(wgt*sum(divergence_real**2.0_pReal)) ! RMS in real space - err_post_div_RMS = sqrt(wgt*sum(divergence_post**2.0_pReal)) ! RMS in real space - err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space - err_div_max = sqrt( err_div_max) ! max in Fourier space - - write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS - write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS - write(6,'(1x,a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS - write(6,'(1x,a,es11.4)') 'error divergence FT max = ',err_div_max - write(6,'(1x,a,es11.4)') 'error divergence Real max = ',err_real_div_max - endif + err_real_div_RMS = sqrt(wgt*sum(divergence_real**2.0_pReal)) ! RMS in real space + err_post_div_RMS = sqrt(wgt*sum(divergence_post**2.0_pReal)) ! RMS in real space from funtion in math.f90 + err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space + err_div_max = sqrt( err_div_max) ! max in Fourier space + + write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS + write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS + write(6,'(1x,a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS + write(6,'(1x,a,es11.4)') 'error divergence FT max = ',err_div_max + write(6,'(1x,a,es11.4)') 'error divergence Real max = ',err_real_div_max + endif -end function Utilities_divergenceRMS +end function utilities_divergenceRMS !-------------------------------------------------------------------------------------------------- -!> @brief calculates mask compliance +!> @brief calculates mask compliance tensor !-------------------------------------------------------------------------------------------------- -function Utilities_maskedCompliance(rot_BC,mask_stress,C) - +function utilities_maskedCompliance(rot_BC,mask_stress,C) + use IO, only: & + IO_error + use math, only: & + math_Plain3333to99, & + math_plain99to3333, & + math_rotate_forward3333, & + math_rotate_forward33, & + math_invert + implicit none - real(pReal), dimension(3,3,3,3) :: Utilities_maskedCompliance - real(pReal), dimension(3,3,3,3), intent(in) :: C + real(pReal), dimension(3,3,3,3) :: utilities_maskedCompliance !< masked compliance + real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness + real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame + logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC integer(pInt) :: j, k, m, n - real(pReal), dimension(3,3), intent(in) :: rot_BC - logical, dimension(3,3), intent(in) :: mask_stress logical, dimension(9) :: mask_stressVector - real(pReal), dimension(3,3,3,3) :: C_lastInc real(pReal), dimension(9,9) :: temp99_Real integer(pInt) :: size_reduced = 0_pInt - real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced, sTimesC ! reduced compliance and stiffness (only for stress BC) + real(pReal), dimension(:,:), allocatable :: & + s_reduced, & !< reduced compliance matrix (depending on number of stress BC) + c_reduced, & !< reduced stiffness (depending on number of stress BC) + sTimesC !< temp variable to check inversion logical :: errmatinv character(len=1024):: formatString @@ -477,11 +544,11 @@ function Utilities_maskedCompliance(rot_BC,mask_stress,C) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) - C_lastInc = math_rotate_forward3333(C,rot_BC) ! calculate stiffness from former inc - print*,'C' - print*,C_lastInc - temp99_Real = math_Plain3333to99(C_lastInc) - k = 0_pInt ! build reduced stiffness + temp99_Real = math_Plain3333to99(math_rotate_forward3333(C,rot_BC)) + if(debugGeneral) & + write(6,'(a,/,9(9(2x,f12.7,1x)/))',advance='no') 'Stiffness C rotated / GPa =',& + transpose(temp99_Real)/1.e9_pReal + k = 0_pInt ! calculate reduced stiffness do n = 1_pInt,9_pInt if(mask_stressVector(n)) then k = k + 1_pInt @@ -493,7 +560,7 @@ function Utilities_maskedCompliance(rot_BC,mask_stress,C) endif; enddo; endif; enddo call math_invert(size_reduced, c_reduced, s_reduced, errmatinv) ! invert reduced stiffness if(errmatinv) call IO_error(error_ID=400_pInt) - temp99_Real = 0.0_pReal ! build full compliance + temp99_Real = 0.0_pReal ! fill up compliance with zeros k = 0_pInt do n = 1_pInt,9_pInt if(mask_stressVector(n)) then @@ -504,14 +571,16 @@ function Utilities_maskedCompliance(rot_BC,mask_stress,C) j = j + 1_pInt temp99_Real(n,m) = s_reduced(k,j) endif; enddo; endif; enddo +!-------------------------------------------------------------------------------------------------- +! check if inversion was successfull sTimesC = matmul(c_reduced,s_reduced) do m=1_pInt, size_reduced do n=1_pInt, size_reduced - if(m==n .and. abs(sTimesC(m,n)) > (1.0_pReal + 10.0e-12_pReal)) errmatinv = .true. - if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. + if(m==n .and. abs(sTimesC(m,n)) > (1.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! diagonal elements of S*C should be 1 + if(m/=n .and. abs(sTimesC(m,n)) > (0.0_pReal + 10.0e-12_pReal)) errmatinv = .true. ! off diagonal elements of S*C should be 0 enddo enddo - if(debugGeneral .or. errmatinv) then + if(debugGeneral .or. errmatinv) then ! report write(formatString, '(I16.16)') size_reduced formatString = '(a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' write(6,trim(formatString),advance='no') 'C * S', transpose(matmul(c_reduced,s_reduced)) @@ -524,161 +593,201 @@ function Utilities_maskedCompliance(rot_BC,mask_stress,C) else temp99_real = 0.0_pReal endif - Utilities_maskedCompliance = math_Plain99to3333(temp99_Real) - print*,'masked S' - print*,Utilities_maskedCompliance -end function Utilities_maskedCompliance + if(debugGeneral) & ! report + write(6,'(a,/,9(9(2x,f12.7,1x)/))',advance='no') 'Masked Compliance * GPa =', & + transpose(temp99_Real*1.e9_pReal) + utilities_maskedCompliance = math_Plain99to3333(temp99_Real) -subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& - P,C,P_av,ForwardData,rotation_BC) - use debug, only: & - debug_reset, & - debug_info - use CPFEM, only: & - CPFEM_general - use FEsolving, only: restartWrite - - implicit none - real(pReal), dimension(res(1),res(2),res(3)) :: temperature - real(pReal), dimension(res(1),res(2),res(3),3), intent(in) :: coordinates - - real(pReal), dimension(3,3,res(1),res(2),res(3)), intent(in) :: F,F_lastInc - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P - real(pReal),intent(in) :: timeinc - logical, intent(in) :: forwardData - integer(pInt) :: i, j, k, ielem - integer(pInt) :: calcMode, collectMode - real(pReal), dimension(3,3,3,3) :: dPdF - real(pReal), dimension(3,3,3,3),intent(out) :: C - real(pReal), dimension(6) :: sigma ! cauchy stress - real(pReal), dimension(6,6) :: dsde - real(pReal), dimension(3,3), intent(in) :: rotation_BC - real(pReal), dimension(3,3),intent(out) :: P_av - - write(6,'(/,a,/)') '... evaluating constitutive response ......................................' - if (forwardData) then - calcMode = 1_pInt - collectMode = 4_pInt - else - calcMode = 2_pInt - collectMode = 3_pInt - endif - if (cutBack) then - calcMode = 2_pInt - collectMode = 5_pInt - endif - - if (DebugGeneral) then - write(6,*) 'collect mode', collectMode - write(6,*) 'calc mode', calcMode - endif - - ielem = 0_pInt - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - ielem = ielem + 1_pInt - call CPFEM_general(collectMode,& ! collect cycle - coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & - temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) - collectMode = 3_pInt - enddo; enddo; enddo - - P = 0.0_pReal ! needed because of the padding for FFTW - C = 0.0_pReal - ielem = 0_pInt - call debug_reset() - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - ielem = ielem + 1_pInt - call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1, - coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) - temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) - calcMode = 2_pInt - C = C + dPdF - enddo; enddo; enddo - call debug_info() - - P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) !average of P rotated - restartWrite = .false. - cutBack = .false. - - write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& - math_transpose33(P_av)/1.e6_pReal - - C = C * wgt -end subroutine Utilities_constitutiveResponse +end function utilities_maskedCompliance -function Utilities_calculateRate(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field) +!-------------------------------------------------------------------------------------------------- +!> @brief calculates constitutive response +!-------------------------------------------------------------------------------------------------- +subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& + P,C,P_av,forwardData,rotation_BC) + use debug, only: & + debug_reset, & + debug_info + use math, only: & + math_transpose33, & + math_rotate_forward33 + use FEsolving, only: & + restartWrite + use mesh, only: & + res, & + wgt + use CPFEM, only: & + CPFEM_general + implicit none - real(pReal), intent(in), dimension(3,3) :: delta_aim - real(pReal), intent(in) :: timeinc, timeinc_old, guessmode - real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: Utilities_calculateRate + real(pReal), intent(inout), dimension(res(1),res(2),res(3)) :: temperature !< temperature field + real(pReal), intent(in), dimension(res(1),res(2),res(3),3) :: coordinates !< coordinates field + real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + F_lastInc, & !< target deformation gradient + F !< previous deformation gradient + real(pReal), intent(in) :: timeinc !< loading time + logical, intent(in) :: forwardData !< age results + real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame + + real(pReal),intent(out), dimension(3,3,3,3) :: C !< average stiffness + real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress + real(pReal),intent(out), dimension(3,3,res(1),res(2),res(3)) :: P !< PK stress + + integer(pInt) :: & + i, j, k, & + ielem, & + calcMode, & !< CPFEM mode for calculation + collectMode !< CPFEM mode for collection + real(pReal), dimension(3,3,3,3) :: dPdF !< d P / d F + real(pReal), dimension(6) :: sigma !< cauchy stress in mandel notation + real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon + + write(6,'(/,a,/)') '... evaluating constitutive response ......................................' + if (forwardData) then ! aging results + calcMode = 1_pInt + collectMode = 4_pInt + else ! normal calculation + calcMode = 2_pInt + collectMode = 3_pInt + endif + if (cutBack) then ! restore saved variables + calcMode = 2_pInt + collectMode = 5_pInt + endif + + if (DebugGeneral) write(6,*) 'collect mode: ', collectMode,' calc mode: ', calcMode + + ielem = 0_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ielem = ielem + 1_pInt + call CPFEM_general(collectMode,& ! collect cycle + coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & + temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) + collectMode = 3_pInt + enddo; enddo; enddo + + P = 0.0_pReal ! needed because of the padding for FFTW + C = 0.0_pReal + ielem = 0_pInt + call debug_reset() + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ielem = ielem + 1_pInt + call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1, + coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) + temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) + calcMode = 2_pInt + C = C + dPdF + enddo; enddo; enddo + C = C * wgt + call debug_info() + + P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) ! average of P rotated + restartWrite = .false. ! reset restartWrite status + cutBack = .false. ! reset cutBack status + + write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& + math_transpose33(P_av)/1.e6_pReal +end subroutine utilities_constitutiveResponse + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates forward rate, either guessing or just add delta/timeinc +!-------------------------------------------------------------------------------------------------- +pure function utilities_calculateRate(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field) + use mesh, only: & + res + + implicit none + real(pReal), intent(in), dimension(3,3) :: delta_aim !< homogeneous addon + real(pReal), intent(in) :: & + timeinc, & !< timeinc of current step + timeinc_old, & !< timeinc of last step + guessmode !< timeinc of current step + real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + field_lastInc, & !< data of previous step + field !< data of current step + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_calculateRate if (guessmode == 1.0_pReal) then - Utilities_calculateRate = (field-field_lastInc) / timeinc_old + utilities_calculateRate = (field-field_lastInc) / timeinc_old else - Utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc + utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc endif -end function Utilities_calculateRate +end function utilities_calculateRate -function Utilities_forwardField(timeinc,aim,field_lastInc,rate) +!-------------------------------------------------------------------------------------------------- +!> @brief forwards a field with a pointwise given rate, ensures that the average matches the aim +!-------------------------------------------------------------------------------------------------- +pure function utilities_forwardField(timeinc,aim,field_lastInc,rate) + use mesh, only: & + res, & + wgt + implicit none - real(pReal), intent(in) :: timeinc - real(pReal), intent(in), dimension(3,3) :: aim - real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,rate - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: Utilities_forwardField - real(pReal), dimension(3,3) :: fieldDiff + real(pReal), intent(in) :: timeinc !< timeinc of current step + real(pReal), intent(in), dimension(3,3) :: aim !< average field value aim + real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + field_lastInc,& !< initial field + rate !< rate by which to forward + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_forwardField + real(pReal), dimension(3,3) :: fieldDiff !< - aim - Utilities_forwardField = field_lastInc + rate*timeinc - fieldDiff = sum(sum(sum(Utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim - Utilities_forwardField = Utilities_forwardField - & + utilities_forwardField = field_lastInc + rate*timeinc + fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim + utilities_forwardField = utilities_forwardField - & spread(spread(spread(fieldDiff,3,res(1)),4,res(2)),5,res(3)) -end function Utilities_forwardField +end function utilities_forwardField -real(pReal) function Utilities_getFilter(k) - use numerics, only: & +!-------------------------------------------------------------------------------------------------- +!> @brief calculates filter for fourier convolution depending on type given in numerics.config +!-------------------------------------------------------------------------------------------------- +real(pReal) function utilities_getFilter(k) + use IO, only: & + IO_error + use numerics, only: & myfilter + use mesh, only: & + res + use math, only: & + PI - implicit none + implicit none + real(pReal),intent(in), dimension(3) :: k !< indices of frequency - real(pReal), dimension(3),intent(in) :: k - - select case (myfilter) - + select case (myfilter) case ('none') - Utilities_getFilter = 1.0_pReal - - case ('cosine') - Utilities_getFilter = (1.0_pReal + cos(pi*k(3)/res(3))) & - *(1.0_pReal + cos(pi*k(2)/res(2))) & - *(1.0_pReal + cos(pi*k(1)/res(1)))/8.0_pReal - + utilities_getFilter = 1.0_pReal + case ('cosine') !< cosine curve with 1 for avg and zero for highest freq + utilities_getFilter = (1.0_pReal + cos(PI*k(3)/res(3))) & + *(1.0_pReal + cos(PI*k(2)/res(2))) & + *(1.0_pReal + cos(PI*k(1)/res(1)))/8.0_pReal case default call IO_error(error_ID = 892_pInt, ext_msg = trim(myfilter)) - end select -end function Utilities_getFilter +end function utilities_getFilter -subroutine Utilities_destroy() +!-------------------------------------------------------------------------------------------------- +!> @brief cleans up +!-------------------------------------------------------------------------------------------------- +subroutine utilities_destroy() + use math + implicit none + if (debugDivergence) call fftw_destroy_plan(plan_divergence) - implicit none - if (debugDivergence) call fftw_destroy_plan(plan_divergence) + if (debugFFTW) call fftw_destroy_plan(plan_scalarField_forth) + if (debugFFTW) call fftw_destroy_plan(plan_scalarField_back) - if (debugFFTW) then - call fftw_destroy_plan(plan_scalarField_forth) - call fftw_destroy_plan(plan_scalarField_back) - endif - - call fftw_destroy_plan(plan_forward) - call fftw_destroy_plan(plan_backward) + call fftw_destroy_plan(plan_forward) + call fftw_destroy_plan(plan_backward) -end subroutine Utilities_destroy +end subroutine utilities_destroy -end module DAMASK_spectral_Utilities +end module DAMASK_spectral_utilities