diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 deleted file mode 100644 index cf8856392..000000000 --- a/code/DAMASK_spectral_solverBasic.f90 +++ /dev/null @@ -1,391 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -! $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 Basic scheme solver -!> @details this solver follows closely the original large strain formulation presented by -!> Suquet. The iterative procedure is solved using a fix-point iteration. -!-------------------------------------------------------------------------------------------------- -module DAMASK_spectral_SolverBasic - use prec, only: & - pInt, & - pReal - use math, only: & - math_I3 - use DAMASK_spectral_Utilities, only: & - tSolutionState - - implicit none - private - character (len=*), parameter, public :: & - DAMASK_spectral_SolverBasic_label = 'basic' - -!-------------------------------------------------------------------------------------------------- -! pointwise global data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: & - F, & !< deformation gradient field - F_lastInc, & !< deformation gradient field last increment - Fdot !< assumed rate for F n to F n+1 - -!-------------------------------------------------------------------------------------------------- -! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aim = math_I3, & !< deformation gradient aim - F_aim_lastInc = math_I3, & !< deformation gradient aim last increment - F_aimDot = 0.0_pReal !< assumed rate - real(pReal), private,dimension(3,3,3,3) :: & - C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, & !< average stiffness - C_lastInc = 0.0_pReal !< average stiffness last increment - public :: & - basic_init, & - basic_solution, & - basic_destroy -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields and fills them with data, potentially from restart info -!-------------------------------------------------------------------------------------------------- -subroutine basic_init(temperature) - 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_read_realFile, & - IO_intOut, & - IO_timeStamp - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRestart - use FEsolving, only: & - restartInc - use DAMASK_interface, only: & - getSolverJobName - use DAMASK_spectral_Utilities, only: & - Utilities_init, & - Utilities_constitutiveResponse, & - Utilities_updateGamma, & - grid, & - wgt, & - geomSize - use mesh, only: & - mesh_ipCoordinates, & - mesh_deformedCoordsFFT - - implicit none - real(pReal), intent(inout) :: & - temperature - real(pReal), dimension(:,:,:,:,:), allocatable :: P - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - - call Utilities_Init() - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' - write(6,'(a)') ' $Id$' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - allocate (P (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) -!-------------------------------------------------------------------------------------------------- -! allocate global fields - allocate (F (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) - allocate (F_lastInc (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) - allocate (Fdot (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal) - -!-------------------------------------------------------------------------------------------------- -! init fields and average quantities - if (restartInc == 1_pInt) then ! no deformation (no restart) - F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity - F_lastInc = F - elseif (restartInc > 1_pInt) then ! using old values from file - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment ', restartInc - 1_pInt, ' from file' - flush(6) - call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) - call IO_read_realFile(777,'F_lastInc',trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) - call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot)) - read (777,rec=1) f_aimDot - close (777) - F_aim = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt ! average of F - F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc - endif - - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)]) - call Utilities_constitutiveResponse(F_lastInc,F,temperature,0.0_pReal,P,C,C_minMaxAvg,& - temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness - - if (restartInc > 1_pInt) then ! using old values from files - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) & - write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading more values of increment', restartInc - 1_pInt, 'from file' - flush(6) - call IO_read_realFile(777,'C',trim(getSolverJobName()),size(C)) - read (777,rec=1) C - close (777) - call IO_read_realFile(777,'C_lastInc',trim(getSolverJobName()),size(C_lastInc)) - read (777,rec=1) C_lastInc - close (777) - call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) - read (777,rec=1) C_minMaxAvg - close (777) - endif - - call Utilities_updateGamma(C_minmaxAvg,.True.) - -end subroutine basic_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief solution for the basic scheme with internal iterations -!-------------------------------------------------------------------------------------------------- -type(tSolutionState) function basic_solution(& - incInfo,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC) - use numerics, only: & - itmax, & - itmin, & - update_gamma - use math, only: & - math_mul33x33 ,& - math_rotate_backward33, & - math_transpose33, & - math_mul3333xx33 - use mesh, only: & - mesh_ipCoordinates,& - mesh_deformedCoordsFFT - use IO, only: & - IO_write_jobRealFile, & - IO_intOut - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use DAMASK_spectral_Utilities, only: & - tBoundaryCondition, & - field_real, & - Utilities_forwardField, & - Utilities_maskedCompliance, & - Utilities_FFTforward, & - Utilities_divergenceRMS, & - Utilities_fourierConvolution, & - Utilities_FFTbackward, & - Utilities_updateGamma, & - Utilities_constitutiveResponse, & - Utilities_calculateRate, & - grid,& - geomSize, & - wgt - use FEsolving, only: & - restartWrite, & - restartRead, & - terminallyIll - use DAMASK_spectral_Utilities, only: & - cutBack - - implicit none -!-------------------------------------------------------------------------------------------------- -! input data for solution - real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old, & !< increment in time of last increment - loadCaseTime !< remaining time of current load case - logical, intent(in) :: & - guess !< if .false., assume homogeneous addon - type(tBoundaryCondition), intent(in) :: & - P_BC,& !< stress boundary conditions - F_BC !< deformation boundary conditions - character(len=*), intent(in) :: & - incInfo !< string with information of current increment for output to screen - real(pReal), dimension(3,3), intent(in) :: & - rotation_BC !< rotation load to lab - real(pReal), intent(inout) :: & - temperature_bc - real(pReal), dimension(3,3,3,3) :: & - S !< current average compliance - real(pReal), dimension(3,3) :: & - F_aim_lastIter, & !< aim of last iteration - P_av - real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: P -!-------------------------------------------------------------------------------------------------- -! loop variables, convergence etc. - real(pReal) :: err_div, err_stress - integer(pInt) :: iter - logical :: ForwardData - -!-------------------------------------------------------------------------------------------------- -! write restart information for spectral solver - if (restartWrite) then - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - call IO_write_jobRealFile(777,'F',size(F)) ! writing deformation gradient field to file - write (777,rec=1) F - close (777) - call IO_write_jobRealFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file - write (777,rec=1) F_lastInc - close (777) - call IO_write_jobRealFile(777,'F_aimDot',size(f_aimDot)) - write (777,rec=1) f_aimDot - close(777) - call IO_write_jobRealFile(777,'C',size(C)) - write (777,rec=1) C - close(777) - call IO_write_jobRealFile(777,'C_lastInc',size(C_lastInc)) - write (777,rec=1) C_lastInc - close(777) - endif - -!-------------------------------------------------------------------------------------------------- -! forward data - if (cutBack) then - F_aim = F_aim_lastInc - F = F_lastInc - C = C_lastInc - else - C_lastInc = C - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)]) - -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values - F_aim)/loadCaseTime - endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old - F_aim_lastInc = F_aim - -!-------------------------------------------------------------------------------------------------- -! update rate and forward last inc - Fdot = utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lastInc,F) - F_lastInc = F - endif - F_aim = F_aim + f_aimDot * timeinc - F = Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim,rotation_BC)) - -!-------------------------------------------------------------------------------------------------- -! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) - if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) - -!-------------------------------------------------------------------------------------------------- -! iteration till converged - if (.not. restartRead) ForwardData = .True. - iter = 0_pInt - basic_solution%iterationsNeeded = itmax - convergenceLoop: do while(iter < itmax) - iter = iter + 1_pInt -!-------------------------------------------------------------------------------------------------- -! report begin of new iteration - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & - ' @ Iteration ', itmin, '≤',iter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab)=', & - math_transpose33(math_rotate_backward33(F_aim,rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) - flush(6) -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - F_aim_lastIter = F_aim - basic_solution%termIll = .false. - call Utilities_constitutiveResponse(F_lastInc,F,temperature_bc,timeinc,& - P,C,C_minmaxAvg,P_av,ForwardData,rotation_BC) - basic_solution%termIll = terminallyIll - terminallyIll = .false. - ForwardData = .false. - -!-------------------------------------------------------------------------------------------------- -! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) ! S = 0.0 for no bc - err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc - -!-------------------------------------------------------------------------------------------------- -! updated deformation gradient using fix point algorithm of basic scheme - field_real = 0.0_pReal - field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(P,[grid(1),grid(2),grid(3),3,3],& - order=[4,5,1,2,3]) ! field real has a different order - call Utilities_FFTforward() - err_div = Utilities_divergenceRMS() - call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,rotation_BC)) - call Utilities_FFTbackward() - F = F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization - basic_solution%converged = basic_Converged(err_div,P_av,err_stress) - write(6,'(/,a)') ' ===========================================================================' - flush(6) - if ((basic_solution%converged .and. iter >= itmin) .or. basic_solution%termIll) then - basic_solution%iterationsNeeded = iter - exit - endif - enddo convergenceLoop - -end function basic_solution - - -!-------------------------------------------------------------------------------------------------- -!> @brief convergence check for basic scheme based on div of P and deviation from stress aim -!-------------------------------------------------------------------------------------------------- -logical function basic_Converged(err_div,P_av,err_stress) - use numerics, only: & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_stress_tolrel, & - err_stress_tolabs - use math, only: & - math_mul33x33, & - math_eigenvalues33, & - math_transpose33 - - implicit none - real(pReal), dimension(3,3), intent(in) :: & - P_av - - real(pReal), intent(in) :: & - err_div, & - err_stress - - real(pReal) :: & - stressTol, & - divTol - - divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) - stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) - - basic_Converged = all([ err_div/divTol,& - err_stress/stressTol ] < 1.0_pReal) - -!-------------------------------------------------------------------------------------------------- -! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div,' / m, tol =',divTol,')' - write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', & - err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')' - flush(6) - -end function basic_Converged - - -!-------------------------------------------------------------------------------------------------- -!> @brief does the cleaning up after the simulation has finished -!-------------------------------------------------------------------------------------------------- -subroutine basic_destroy() - -use DAMASK_spectral_Utilities, only: & - Utilities_destroy - - implicit none - call Utilities_destroy() - -end subroutine basic_destroy - - -end module DAMASK_spectral_SolverBasic diff --git a/code/numerics.f90 b/code/numerics.f90 index 2f77d20ea..db2d90ec4 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -87,7 +87,7 @@ module numerics character(len=64), private :: & fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag character(len=64), protected, public :: & - spectral_solver = 'basic' , & !< spectral solution method + spectral_solver = 'basicpetsc' , & !< spectral solution method spectral_filter = 'none' !< spectral filtering method character(len=1024), protected, public :: & petsc_options = '-snes_type ngmres &