diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index dcbbe444a..a0cd44faa 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,8 @@ if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240") + SET_SOURCE_FILES_PROPERTIES( "DAMASK_interface.f90" PROPERTIES + COMPILE_FLAGS "-ffree-line-length-164") # long lines for interaction matrix endif() @@ -172,17 +174,17 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral") list(APPEND OBJECTFILES $) add_library(SPECTRAL_SOLVER OBJECT - "spectral_thermal.f90" - "spectral_damage.f90" - "spectral_mech_Polarisation.f90" - "spectral_mech_Basic.f90") + "grid_thermal_spectral.f90" + "grid_damage_spectral.f90" + "grid_mech_spectral_basic.f90" + "grid_mech_spectral_polarisation.f90") add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES) list(APPEND OBJECTFILES $) if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY") - add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES}) + add_executable(DAMASK_spectral "DAMASK_grid.f90" ${OBJECTFILES}) else() - add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90") + add_library(DAMASK_spectral OBJECT "DAMASK_grid.f90") endif() add_dependencies(DAMASK_spectral SPECTRAL_SOLVER) diff --git a/src/DAMASK_abaqus.f b/src/DAMASK_abaqus.f index 74f1a292f..514307fe8 100644 --- a/src/DAMASK_abaqus.f +++ b/src/DAMASK_abaqus.f @@ -49,15 +49,15 @@ subroutine DAMASK_interface_init write(6,'(/,a)') ' <<<+- DAMASK_abaqus init -+>>>' - write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2018' + write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(/,a)') ' Version: '//DAMASKVERSION ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md #if __INTEL_COMPILER >= 1800 - write(6,'(/,a)') 'Compiled with: '//compiler_version() - write(6,'(a)') 'Compiler options: '//compiler_options() + write(6,'(/,a)') ' Compiled with: '//compiler_version() + write(6,'(a)') ' Compiler options: '//compiler_options() #else write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& ', build date :', __INTEL_COMPILER_BUILD_DATE diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_grid.f90 similarity index 96% rename from src/DAMASK_spectral.f90 rename to src/DAMASK_grid.f90 index c53d339df..496bfd0de 100644 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_grid.f90 @@ -7,11 +7,6 @@ !> results !-------------------------------------------------------------------------------------------------- program DAMASK_spectral -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif #include use PETScsys use prec, only: & @@ -35,8 +30,7 @@ program DAMASK_spectral IO_error, & IO_lc, & IO_intOut, & - IO_warning, & - IO_timeStamp + IO_warning use debug, only: & debug_level, & debug_spectral, & @@ -77,10 +71,10 @@ program DAMASK_spectral FIELD_MECH_ID, & FIELD_THERMAL_ID, & FIELD_DAMAGE_ID - use spectral_mech_Basic - use spectral_mech_Polarisation - use spectral_damage - use spectral_thermal + use grid_mech_spectral_basic + use grid_mech_spectral_polarisation + use grid_damage_spectral + use grid_thermal_spectral use results implicit none @@ -141,11 +135,11 @@ program DAMASK_spectral integer(pInt), parameter :: maxRealOut = maxByteOut/pReal integer(pLongInt), dimension(2) :: outputIndex PetscErrorCode :: ierr - procedure(basic_init), pointer :: & + procedure(grid_mech_spectral_basic_init), pointer :: & mech_init - procedure(basic_forward), pointer :: & + procedure(grid_mech_spectral_basic_forward), pointer :: & mech_forward - procedure(basic_solution), pointer :: & + procedure(grid_mech_spectral_basic_solution), pointer :: & mech_solution external :: & @@ -155,10 +149,9 @@ program DAMASK_spectral ! init DAMASK (all modules) call CPFEM_initAll write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>' - write(6,'(/,a,/)') ' Roters et al., Computational Materials Science, 2018' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" + write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' + write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' call results_openJobFile() call results_closeJobFile() @@ -173,17 +166,17 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type select case (spectral_solver) - case (DAMASK_spectral_SolverBasic_label) - mech_init => basic_init - mech_forward => basic_forward - mech_solution => basic_solution + case (GRID_MECH_SPECTRAL_BASIC_LABEL) + mech_init => grid_mech_spectral_basic_init + mech_forward => grid_mech_spectral_basic_forward + mech_solution => grid_mech_spectral_basic_solution - case (DAMASK_spectral_SolverPolarisation_label) + case (GRID_MECH_SPECTRAL_POLARISATION_LABEL) if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) & call IO_warning(42_pInt, ext_msg='debug Divergence') - mech_init => polarisation_init - mech_forward => polarisation_forward - mech_solution => polarisation_solution + mech_init => grid_mech_spectral_polarisation_init + mech_forward => grid_mech_spectral_polarisation_forward + mech_solution => grid_mech_spectral_polarisation_solution case default call IO_error(error_ID = 891_pInt, ext_msg = trim(spectral_solver)) @@ -365,10 +358,10 @@ program DAMASK_spectral call mech_init case(FIELD_THERMAL_ID) - call spectral_thermal_init + call grid_thermal_spectral_init case(FIELD_DAMAGE_ID) - call spectral_damage_init + call grid_damage_spectral_init end select enddo @@ -510,8 +503,8 @@ program DAMASK_spectral stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rotation) - case(FIELD_THERMAL_ID); call spectral_thermal_forward() - case(FIELD_DAMAGE_ID); call spectral_damage_forward() + case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward + case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward end select enddo @@ -529,10 +522,10 @@ program DAMASK_spectral rotation_BC = loadCases(currentLoadCase)%rotation) case(FIELD_THERMAL_ID) - solres(field) = spectral_thermal_solution(timeinc,timeIncOld,remainingLoadCaseTime) + solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld,remainingLoadCaseTime) case(FIELD_DAMAGE_ID) - solres(field) = spectral_damage_solution(timeinc,timeIncOld,remainingLoadCaseTime) + solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld,remainingLoadCaseTime) end select diff --git a/src/DAMASK_interface.f90 b/src/DAMASK_interface.f90 index d0c5c426c..143561f5d 100644 --- a/src/DAMASK_interface.f90 +++ b/src/DAMASK_interface.f90 @@ -152,15 +152,24 @@ subroutine DAMASK_interface_init() write(6,'(/,a)') ' <<<+- DAMASK_interface init -+>>>' - write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2018' + ! http://patorjk.com/software/taag/#p=display&f=Lean&t=DAMASK + write(6,*) achar(27)//'[94m' + write(6,*) ' _/_/_/ _/_/ _/ _/ _/_/ _/_/_/ _/ _/' + write(6,*) ' _/ _/ _/ _/ _/_/ _/_/ _/ _/ _/ _/ _/' + write(6,*) ' _/ _/ _/_/_/_/ _/ _/ _/ _/_/_/_/ _/_/ _/_/' + write(6,*) ' _/ _/ _/ _/ _/ _/ _/ _/ _/ _/ _/' + write(6,*) ' _/_/_/ _/ _/ _/ _/ _/ _/ _/_/_/ _/ _/' + write(6,*) achar(27)//'[0m' + + write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(/,a)') ' Version: '//DAMASKVERSION ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - write(6,'(/,a)') 'Compiled with: '//compiler_version() - write(6,'(a)') 'Compiler options: '//compiler_options() + write(6,'(/,a)') ' Compiled with: '//compiler_version() + write(6,'(a)') ' Compiler options: '//compiler_options() #elif defined(__INTEL_COMPILER) write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& ', build date :', __INTEL_COMPILER_BUILD_DATE diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index b447690bf..edf00a203 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -60,15 +60,15 @@ subroutine DAMASK_interface_init write(6,'(/,a)') ' <<<+- DAMASK_marc init -+>>>' - write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2018' + write(6,'(/,a)') ' Roters et al., Computational Materials Science 158:420–478, 2019' write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2018.04.030' write(6,'(/,a)') ' Version: '//DAMASKVERSION ! https://github.com/jeffhammond/HPCInfo/blob/master/docs/Preprocessor-Macros.md #if __INTEL_COMPILER >= 1800 - write(6,'(/,a)') 'Compiled with: '//compiler_version() - write(6,'(a)') 'Compiler options: '//compiler_options() + write(6,'(/,a)') ' Compiled with: '//compiler_version() + write(6,'(a)') ' Compiler options: '//compiler_options() #else write(6,'(/,a,i4.4,a,i8.8)') ' Compiled with Intel fortran version :', __INTEL_COMPILER,& ', build date :', __INTEL_COMPILER_BUILD_DATE diff --git a/src/FEM_mech.f90 b/src/FEM_mech.f90 index b6d9ac17f..68a24d548 100644 --- a/src/FEM_mech.f90 +++ b/src/FEM_mech.f90 @@ -14,10 +14,7 @@ module FEM_mech use PETScDMplex use PETScDT use prec, only: & - pInt, & pReal - use math, only: & - math_I3 use FEM_utilities, only: & tSolutionState, & tFieldBC, & @@ -88,8 +85,8 @@ subroutine FEM_mech_init(fieldBC) PetscDS :: mechDS PetscDualSpace :: mechDualSpace DMLabel :: BCLabel - PetscInt, allocatable, target :: numComp(:), numDoF(:), bcField(:) - PetscInt, pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:) + PetscInt, dimension(:), allocatable, target :: numComp, numDoF, bcField + PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt :: numBC, bcSize, nc IS :: bcPoint IS, allocatable, target :: bcComps(:), bcPoints(:) @@ -467,20 +464,20 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) PetscReal, target :: v0(dimPlex), cellJ(dimPlex*dimPlex), & invcellJ(dimPlex*dimPlex) - PetscReal, pointer :: pV0(:), pCellJ(:), pInvcellJ(:) - PetscReal, dimension(:), pointer :: basisField, basisFieldDer + PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & + pV0, pCellJ, pInvcellJ PetscInt :: cellStart, cellEnd, cell, field, face, & qPt, basis, comp, cidx - PetscScalar, target :: K_e (cellDof,cellDof), & - K_eA (cellDof,cellDof), & - K_eB (cellDof,cellDof), & - K_eVec(cellDof*cellDof) + PetscScalar,dimension(cellDOF,cellDOF), target :: K_e, & + K_eA , & + K_eB + PetscScalar, target :: K_eVec(cellDof*cellDof) PetscReal :: BMat (dimPlex*dimPlex,cellDof), & BMatAvg(dimPlex*dimPlex,cellDof), & MatA (dimPlex*dimPlex,cellDof), & MatB (1 ,cellDof) PetscScalar, dimension(:), pointer :: pK_e, x_scal - PetscReal, dimension(3,3) :: F = math_I3, FAvg, FInv + PetscReal, dimension(3,3) :: F, FAvg, FInv PetscObject :: dummy PetscInt :: bcSize IS :: bcPoints @@ -619,11 +616,11 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) if (guess .and. .not. cutBack) then ForwardData = .True. materialpoint_F0 = materialpoint_F - call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local + call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call VecSet(x_local,0.0,ierr); CHKERRQ(ierr) - call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector + call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector CHKERRQ(ierr) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr) CHKERRQ(ierr) diff --git a/src/FEM_zoo.f90 b/src/FEM_zoo.f90 index 6abdfe883..53ef4655d 100644 --- a/src/FEM_zoo.f90 +++ b/src/FEM_zoo.f90 @@ -34,20 +34,9 @@ contains !> @brief initializes FEM interpolation data !-------------------------------------------------------------------------------------------------- subroutine FEM_Zoo_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use IO, only: & - IO_timeStamp - implicit none write(6,'(/,a)') ' <<<+- FEM_Zoo init -+>>>' - write(6,'(a)') ' $Id: FEM_Zoo.f90 4354 2015-08-04 15:04:53Z MPIE\p.shanthraj $' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" !-------------------------------------------------------------------------------------------------- ! 2D linear diff --git a/src/FEsolving.f90 b/src/FEsolving.f90 index f31500c26..d63617135 100644 --- a/src/FEsolving.f90 +++ b/src/FEsolving.f90 @@ -43,11 +43,6 @@ contains !> solver the information is provided by the interface module !-------------------------------------------------------------------------------------------------- subroutine FE_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif use debug, only: & debug_level, & debug_FEsolving, & @@ -61,8 +56,7 @@ subroutine FE_init IO_open_inputFile, & IO_open_logFile, & #endif - IO_warning, & - IO_timeStamp + IO_warning use DAMASK_interface implicit none @@ -75,8 +69,6 @@ subroutine FE_init #endif write(6,'(/,a)') ' <<<+- FEsolving init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" modelName = getSolverJobName() diff --git a/src/IO.f90 b/src/IO.f90 index 983465caa..cc90cbbb2 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -22,7 +22,6 @@ module IO public :: & IO_init, & IO_read_ASCII, & - IO_recursiveRead, & IO_open_file, & IO_open_jobFile_binary, & IO_write_jobFile, & @@ -35,8 +34,7 @@ module IO IO_lc, & IO_error, & IO_warning, & - IO_intOut, & - IO_timeStamp + IO_intOut #if defined(Marc4DAMASK) || defined(Abaqus) public :: & IO_open_inputFile, & @@ -160,89 +158,6 @@ function IO_read_ASCII(fileName) result(fileContent) end function IO_read_ASCII -!-------------------------------------------------------------------------------------------------- -!> @brief recursively reads a text file. -!! Recursion is triggered by "{path/to/inputfile}" in a line -!-------------------------------------------------------------------------------------------------- -recursive function IO_recursiveRead(fileName,cnt) result(fileContent) - - implicit none - character(len=*), intent(in) :: fileName - integer(pInt), intent(in), optional :: cnt !< recursion counter - character(len=256), dimension(:), allocatable :: fileContent !< file content, separated per lines - character(len=256), dimension(:), allocatable :: includedContent - character(len=256) :: line - character(len=256), parameter :: dummy = 'https://damask.mpie.de' !< to fill up remaining array - character(len=:), allocatable :: rawData - integer(pInt) :: & - fileLength, & - fileUnit, & - startPos, endPos, & - myTotalLines, & !< # lines read from file without include statements - l,i, & - myStat - logical :: warned - - if (present(cnt)) then - if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) - endif - -!-------------------------------------------------------------------------------------------------- -! read data as stream - inquire(file = fileName, size=fileLength) - if (fileLength == 0) then - allocate(fileContent(0)) - return - endif - open(newunit=fileUnit, file=fileName, access='stream',& - status='old', position='rewind', action='read',iostat=myStat) - if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(fileName)) - allocate(character(len=fileLength)::rawData) - read(fileUnit) rawData - close(fileUnit) - -!-------------------------------------------------------------------------------------------------- -! count lines to allocate string array - myTotalLines = 1_pInt - do l=1_pInt, len(rawData) - if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1 - enddo - allocate(fileContent(myTotalLines)) - -!-------------------------------------------------------------------------------------------------- -! split raw data at end of line and handle includes - warned = .false. - startPos = 1_pInt - l = 1_pInt - do while (l <= myTotalLines) - endPos = merge(startPos + scan(rawData(startPos:),new_line('')) - 2_pInt,len(rawData),l /= myTotalLines) - if (endPos - startPos > 255_pInt) then - line = rawData(startPos:startPos+255_pInt) - if (.not. warned) then - call IO_warning(207_pInt,ext_msg=trim(fileName),el=l) - warned = .true. - endif - else - line = rawData(startPos:endpos) - endif - startPos = endPos + 2_pInt ! jump to next line start - - recursion: if (scan(trim(adjustl(line)),'{') == 1 .and. scan(trim(line),'}') > 2) then - includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), & - merge(cnt,1_pInt,present(cnt))) ! to track recursion depth - fileContent = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,myTotalLines-l)] ] ! add content and grow array - myTotalLines = myTotalLines - 1_pInt + size(includedContent) - l = l - 1_pInt + size(includedContent) - else recursion - fileContent(l) = line - l = l + 1_pInt - endif recursion - - enddo - -end function IO_recursiveRead - - !-------------------------------------------------------------------------------------------------- !> @brief opens existing file for reading to given unit. Path to file is relative to working !! directory @@ -717,21 +632,6 @@ pure function IO_intOut(intToPrint) end function IO_intOut -!-------------------------------------------------------------------------------------------------- -!> @brief returns time stamp -!-------------------------------------------------------------------------------------------------- -function IO_timeStamp() - - implicit none - character(len=10) :: IO_timeStamp - integer(pInt), dimension(8) :: values - - call DATE_AND_TIME(VALUES=values) - write(IO_timeStamp,'(i2.2,a1,i2.2,a1,i2.2)') values(5),':',values(6),':',values(7) - -end function IO_timeStamp - - !-------------------------------------------------------------------------------------------------- !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx !> in ABAQUS either time step is reduced or execution terminated diff --git a/src/Lambert.f90 b/src/Lambert.f90 index dc2626296..c570a600e 100644 --- a/src/Lambert.f90 +++ b/src/Lambert.f90 @@ -38,29 +38,30 @@ !> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014). !-------------------------------------------------------------------------- module Lambert - use math - use prec, only: & - pReal + use prec, only: & + pReal + use math, only: & + PI - implicit none - private - real(pReal), parameter, private :: & - SPI = sqrt(PI), & - PREF = sqrt(6.0_pReal/PI), & - A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & - AP = PI**(2.0_pReal/3.0_pReal), & - SC = A/AP, & - BETA = A/2.0_pReal, & - R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), & - R2 = sqrt(2.0_pReal), & - PI12 = PI/12.0_pReal, & - PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA - - public :: & - LambertCubeToBall, & - LambertBallToCube - private :: & - GetPyramidOrder + implicit none + private + real(pReal), parameter, private :: & + SPI = sqrt(PI), & + PREF = sqrt(6.0_pReal/PI), & + A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), & + AP = PI**(2.0_pReal/3.0_pReal), & + SC = A/AP, & + BETA = A/2.0_pReal, & + R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), & + R2 = sqrt(2.0_pReal), & + PI12 = PI/12.0_pReal, & + PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA + + public :: & + LambertCubeToBall, & + LambertBallToCube + private :: & + GetPyramidOrder contains @@ -71,56 +72,56 @@ contains !> @brief map from 3D cubic grid to 3D ball !-------------------------------------------------------------------------- function LambertCubeToBall(cube) result(ball) - use, intrinsic :: IEEE_ARITHMETIC - use prec, only: & - pInt, & - dEq0 - - implicit none - real(pReal), intent(in), dimension(3) :: cube - real(pReal), dimension(3) :: ball, LamXYZ, XYZ - real(pReal) :: T(2), c, s, q - real(pReal), parameter :: eps = 1.0e-8_pReal - integer(pInt), dimension(3) :: p - integer(pInt), dimension(2) :: order + use, intrinsic :: IEEE_ARITHMETIC + use prec, only: & + dEq0 - if (maxval(abs(cube)) > AP/2.0+eps) then - ball = IEEE_value(cube,IEEE_positive_inf) - return - end if + implicit none + real(pReal), intent(in), dimension(3) :: cube + real(pReal), dimension(3) :: ball, LamXYZ, XYZ + real(pReal), dimension(2) :: T + real(pReal) :: c, s, q + real(pReal), parameter :: eps = 1.0e-8_pReal + integer, dimension(3) :: p + integer, dimension(2) :: order + + if (maxval(abs(cube)) > AP/2.0+eps) then + ball = IEEE_value(cube,IEEE_positive_inf) + return + end if + + ! transform to the sphere grid via the curved square, and intercept the zero point + center: if (all(dEq0(cube))) then + ball = 0.0_pReal + else center + ! get pyramide and scale by grid parameter ratio + p = GetPyramidOrder(cube) + XYZ = cube(p) * sc + + ! intercept all the points along the z-axis + special: if (all(dEq0(XYZ(1:2)))) then + LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ] + else special + order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ + q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger + c = cos(q) + s = sin(q) + q = prek * XYZ(order(2))/ sqrt(R2-c) + T = [ (R2*c - 1.0), R2 * s] * q + + ! transform to sphere grid (inverse Lambert) + ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] + c = sum(T**2) + s = Pi * c/(24.0*XYZ(3)**2) + c = sPi * c / sqrt(24.0_pReal) / XYZ(3) + q = sqrt( 1.0 - s ) + LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ] + endif special + + ! reverse the coordinates back to the regular order according to the original pyramid number + ball = LamXYZ(p) - ! transform to the sphere grid via the curved square, and intercept the zero point - center: if (all(dEq0(cube))) then - ball = 0.0_pReal - else center - ! get pyramide and scale by grid parameter ratio - p = GetPyramidOrder(cube) - XYZ = cube(p) * sc - - ! intercept all the points along the z-axis - special: if (all(dEq0(XYZ(1:2)))) then - LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ] - else special - order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ - q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger - c = cos(q) - s = sin(q) - q = prek * XYZ(order(2))/ sqrt(R2-c) - T = [ (R2*c - 1.0), R2 * s] * q - - ! transform to sphere grid (inverse Lambert) - ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] - c = sum(T**2) - s = Pi * c/(24.0*XYZ(3)**2) - c = sPi * c / sqrt(24.0_pReal) / XYZ(3) - q = sqrt( 1.0 - s ) - LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ] - endif special - - ! reverse the coordinates back to the regular order according to the original pyramid number - ball = LamXYZ(p) - - endif center + endif center end function LambertCubeToBall @@ -131,57 +132,58 @@ end function LambertCubeToBall !> @brief map from 3D ball to 3D cubic grid !-------------------------------------------------------------------------- pure function LambertBallToCube(xyz) result(cube) - use, intrinsic :: IEEE_ARITHMETIC, only:& - IEEE_positive_inf, & - IEEE_value - use prec, only: & - pInt, & - dEq0 - - implicit none - real(pReal), intent(in), dimension(3) :: xyz - real(pReal), dimension(3) :: cube, xyz1, xyz3 - real(pReal), dimension(2) :: Tinv, xyz2 - real(pReal) :: rs, qxy, q2, sq2, q, tt - integer(pInt), dimension(3) :: p + use, intrinsic :: IEEE_ARITHMETIC, only:& + IEEE_positive_inf, & + IEEE_value + use prec, only: & + dEq0 + use math, only: & + math_clip - rs = norm2(xyz) - if (rs > R1) then - cube = IEEE_value(cube,IEEE_positive_inf) - return - endif + implicit none + real(pReal), intent(in), dimension(3) :: xyz + real(pReal), dimension(3) :: cube, xyz1, xyz3 + real(pReal), dimension(2) :: Tinv, xyz2 + real(pReal) :: rs, qxy, q2, sq2, q, tt + integer, dimension(3) :: p + + rs = norm2(xyz) + if (rs > R1) then + cube = IEEE_value(cube,IEEE_positive_inf) + return + endif + + center: if (all(dEq0(xyz))) then + cube = 0.0_pReal + else center + p = GetPyramidOrder(xyz) + xyz3 = xyz(p) + + ! inverse M_3 + xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) + + ! inverse M_2 + qxy = sum(xyz2**2) + + special: if (dEq0(qxy)) then + Tinv = 0.0_pReal + else special + q2 = qxy + maxval(abs(xyz2))**2 + sq2 = sqrt(q2) + q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) + tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy + Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & + [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & + abs(xyz2(2)) <= abs(xyz2(1))) + endif special + + ! inverse M_1 + xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc + + ! reverst the coordinates back to the regular order according to the original pyramid number + cube = xyz1(p) - center: if (all(dEq0(xyz))) then - cube = 0.0_pReal - else center - p = GetPyramidOrder(xyz) - xyz3 = xyz(p) - - ! inverse M_3 - xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) - - ! inverse M_2 - qxy = sum(xyz2**2) - - special: if (dEq0(qxy)) then - Tinv = 0.0_pReal - else special - q2 = qxy + maxval(abs(xyz2))**2 - sq2 = sqrt(q2) - q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) - tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy - Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & - [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & - abs(xyz2(2)) <= abs(xyz2(1))) - endif special - - ! inverse M_1 - xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc - - ! reverst the coordinates back to the regular order according to the original pyramid number - cube = xyz1(p) - - endif center + endif center end function LambertBallToCube @@ -192,25 +194,23 @@ end function LambertBallToCube !> @brief determine to which pyramid a point in a cubic grid belongs !-------------------------------------------------------------------------- pure function GetPyramidOrder(xyz) - use prec, only: & - pInt - implicit none - real(pReal),intent(in),dimension(3) :: xyz - integer(pInt), dimension(3) :: GetPyramidOrder - - if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & - ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then - GetPyramidOrder = [1,2,3] - else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & - ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then - GetPyramidOrder = [2,3,1] - else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & - ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then - GetPyramidOrder = [3,1,2] - else - GetPyramidOrder = -1 ! should be impossible, but might simplify debugging - end if + implicit none + real(pReal),intent(in),dimension(3) :: xyz + integer, dimension(3) :: GetPyramidOrder + + if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & + ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then + GetPyramidOrder = [1,2,3] + else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & + ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then + GetPyramidOrder = [2,3,1] + else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & + ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then + GetPyramidOrder = [3,1,2] + else + GetPyramidOrder = -1 ! should be impossible, but might simplify debugging + end if end function GetPyramidOrder diff --git a/src/compilation_info.f90 b/src/compilation_info.f90 deleted file mode 100644 index e69de29bb..000000000 diff --git a/src/config.f90 b/src/config.f90 index 94514529d..65f9f8fad 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -7,14 +7,13 @@ !-------------------------------------------------------------------------------------------------- module config use prec, only: & - pReal, & - pInt + pReal implicit none private type, private :: tPartitionedString character(len=:), allocatable :: val - integer(pInt), dimension(:), allocatable :: pos + integer, dimension(:), allocatable :: pos end type tPartitionedString type, private :: tPartitionedStringList @@ -51,6 +50,10 @@ module config config_homogenization, & config_texture, & config_crystallite + + type(tPartitionedStringList), public, protected :: & + config_numerics, & + config_debug character(len=64), dimension(:), allocatable, public, protected :: & phase_name, & !< name of each phase @@ -61,7 +64,7 @@ module config ! ToDo: Remove, use size(config_phase) etc - integer(pInt), public, protected :: & + integer, public, protected :: & material_Nphase, & !< number of phases material_Nhomogenization !< number of homogenizations @@ -74,15 +77,15 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief reads material.config and stores its content per part !-------------------------------------------------------------------------------------------------- -subroutine config_init() +subroutine config_init use prec, only: & pStringLen use DAMASK_interface, only: & getSolverJobName use IO, only: & + IO_read_ASCII, & IO_error, & IO_lc, & - IO_recursiveRead, & IO_getTag use debug, only: & debug_level, & @@ -90,7 +93,7 @@ subroutine config_init() debug_levelBasic implicit none - integer(pInt) :: myDebug,i + integer :: myDebug,i character(len=pStringLen) :: & line, & @@ -104,36 +107,38 @@ subroutine config_init() inquire(file=trim(getSolverJobName())//'.materialConfig',exist=fileExists) if(fileExists) then - fileContent = IO_recursiveRead(trim(getSolverJobName())//'.materialConfig') + write(6,'(/,a)') ' reading '//trim(getSolverJobName())//'.materialConfig'; flush(6) + fileContent = read_materialConfig(trim(getSolverJobName())//'.materialConfig') else inquire(file='material.config',exist=fileExists) - if(.not. fileExists) call IO_error(100_pInt,ext_msg='material.config') - fileContent = IO_recursiveRead('material.config') + if(.not. fileExists) call IO_error(100,ext_msg='material.config') + write(6,'(/,a)') ' reading material.config'; flush(6) + fileContent = read_materialConfig('material.config') endif - do i = 1_pInt, size(fileContent) + do i = 1, size(fileContent) line = trim(fileContent(i)) part = IO_lc(IO_getTag(line,'<','>')) select case (trim(part)) case (trim('phase')) - call parseFile(phase_name,config_phase,line,fileContent(i+1:)) + call parse_materialConfig(phase_name,config_phase,line,fileContent(i+1:)) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Phase parsed'; flush(6) case (trim('microstructure')) - call parseFile(microstructure_name,config_microstructure,line,fileContent(i+1:)) + call parse_materialConfig(microstructure_name,config_microstructure,line,fileContent(i+1:)) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Microstructure parsed'; flush(6) case (trim('crystallite')) - call parseFile(crystallite_name,config_crystallite,line,fileContent(i+1:)) + call parse_materialConfig(crystallite_name,config_crystallite,line,fileContent(i+1:)) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Crystallite parsed'; flush(6) case (trim('homogenization')) - call parseFile(homogenization_name,config_homogenization,line,fileContent(i+1:)) + call parse_materialConfig(homogenization_name,config_homogenization,line,fileContent(i+1:)) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Homogenization parsed'; flush(6) case (trim('texture')) - call parseFile(texture_name,config_texture,line,fileContent(i+1:)) + call parse_materialConfig(texture_name,config_texture,line,fileContent(i+1:)) if (iand(myDebug,debug_levelBasic) /= 0) write(6,'(a)') ' Texture parsed'; flush(6) end select @@ -143,49 +148,143 @@ subroutine config_init() material_Nhomogenization = size(config_homogenization) material_Nphase = size(config_phase) - if (material_Nhomogenization < 1) call IO_error(160_pInt,ext_msg='') - if (size(config_microstructure) < 1) call IO_error(160_pInt,ext_msg='') - if (size(config_crystallite) < 1) call IO_error(160_pInt,ext_msg='') - if (material_Nphase < 1) call IO_error(160_pInt,ext_msg='') - if (size(config_texture) < 1) call IO_error(160_pInt,ext_msg='') + if (material_Nhomogenization < 1) call IO_error(160,ext_msg='') + if (size(config_microstructure) < 1) call IO_error(160,ext_msg='') + if (size(config_crystallite) < 1) call IO_error(160,ext_msg='') + if (material_Nphase < 1) call IO_error(160,ext_msg='') + if (size(config_texture) < 1) call IO_error(160,ext_msg='') + + + inquire(file='numerics.config', exist=fileExists) + if (fileExists) then + write(6,'(/,a)') ' reading numerics.config'; flush(6) + fileContent = IO_read_ASCII('numerics.config') + call parse_debugAndNumericsConfig(config_numerics,fileContent) + endif + + inquire(file='debug.config', exist=fileExists) + if (fileExists) then + write(6,'(/,a)') ' reading debug.config'; flush(6) + fileContent = IO_read_ASCII('debug.config') + call parse_debugAndNumericsConfig(config_debug,fileContent) + endif -end subroutine config_init +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief reads material.config +!! Recursion is triggered by "{path/to/inputfile}" in a line +!-------------------------------------------------------------------------------------------------- +recursive function read_materialConfig(fileName,cnt) result(fileContent) + use IO, only: & + IO_warning + + implicit none + character(len=*), intent(in) :: fileName + integer, intent(in), optional :: cnt !< recursion counter + character(len=pStringLen), dimension(:), allocatable :: fileContent !< file content, separated per lines + character(len=pStringLen), dimension(:), allocatable :: includedContent + character(len=pStringLen) :: line + character(len=pStringLen), parameter :: dummy = 'https://damask.mpie.de' !< to fill up remaining array + character(len=:), allocatable :: rawData + integer :: & + fileLength, & + fileUnit, & + startPos, endPos, & + myTotalLines, & !< # lines read from file without include statements + l,i, & + myStat + logical :: warned + + if (present(cnt)) then + if (cnt>10) call IO_error(106,ext_msg=trim(fileName)) + endif + +!-------------------------------------------------------------------------------------------------- +! read data as stream + inquire(file = fileName, size=fileLength) + if (fileLength == 0) then + allocate(fileContent(0)) + return + endif + open(newunit=fileUnit, file=fileName, access='stream',& + status='old', position='rewind', action='read',iostat=myStat) + if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName)) + allocate(character(len=fileLength)::rawData) + read(fileUnit) rawData + close(fileUnit) + +!-------------------------------------------------------------------------------------------------- +! count lines to allocate string array + myTotalLines = 1 + do l=1, len(rawData) + if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1 + enddo + allocate(fileContent(myTotalLines)) + +!-------------------------------------------------------------------------------------------------- +! split raw data at end of line and handle includes + warned = .false. + startPos = 1 + l = 1 + do while (l <= myTotalLines) + endPos = merge(startPos + scan(rawData(startPos:),new_line('')) - 2,len(rawData),l /= myTotalLines) + if (endPos - startPos > pStringLen -1) then + line = rawData(startPos:startPos+pStringLen-1) + if (.not. warned) then + call IO_warning(207,ext_msg=trim(fileName),el=l) + warned = .true. + endif + else + line = rawData(startPos:endpos) + endif + startPos = endPos + 2 ! jump to next line start + + recursion: if (scan(trim(adjustl(line)),'{') == 1 .and. scan(trim(line),'}') > 2) then + includedContent = read_materialConfig(trim(line(scan(line,'{')+1:scan(line,'}')-1)), & + merge(cnt,1,present(cnt))) ! to track recursion depth + fileContent = [ fileContent(1:l-1), includedContent, [(dummy,i=1,myTotalLines-l)] ] ! add content and grow array + myTotalLines = myTotalLines - 1 + size(includedContent) + l = l - 1 + size(includedContent) + else recursion + fileContent(l) = line + l = l + 1 + endif recursion + + enddo + +end function read_materialConfig !-------------------------------------------------------------------------------------------------- !> @brief parses the material.config file !-------------------------------------------------------------------------------------------------- -subroutine parseFile(sectionNames,part,line, & - fileContent) - use prec, only: & - pStringLen - use IO, only: & - IO_error, & - IO_getTag - +subroutine parse_materialConfig(sectionNames,part,line, & + fileContent) implicit none character(len=64), allocatable, dimension(:), intent(out) :: sectionNames type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part character(len=pStringLen), intent(inout) :: line character(len=pStringLen), dimension(:), intent(in) :: fileContent - integer(pInt), allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section - integer(pInt) :: i, j + integer, allocatable, dimension(:) :: partPosition ! position of [] tags + last line in section + integer :: i, j logical :: echo echo = .false. - if (allocated(part)) call IO_error(161_pInt,ext_msg=trim(line)) + if (allocated(part)) call IO_error(161,ext_msg=trim(line)) allocate(partPosition(0)) - do i = 1_pInt, size(fileContent) + do i = 1, size(fileContent) line = trim(fileContent(i)) if (IO_getTag(line,'<','>') /= '') exit nextSection: if (IO_getTag(line,'[',']') /= '') then partPosition = [partPosition, i] cycle endif nextSection - if (size(partPosition) < 1_pInt) & + if (size(partPosition) < 1) & echo = (trim(IO_getTag(line,'/','/')) == 'echo') .or. echo enddo @@ -194,9 +293,9 @@ subroutine parseFile(sectionNames,part,line, & partPosition = [partPosition, i] ! needed when actually storing content - do i = 1_pInt, size(partPosition) -1_pInt + do i = 1, size(partPosition) -1 sectionNames(i) = trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']'))) - do j = partPosition(i) + 1_pInt, partPosition(i+1) -1_pInt + do j = partPosition(i) + 1, partPosition(i+1) -1 call part(i)%add(trim(adjustl(fileContent(j)))) enddo if (echo) then @@ -205,7 +304,27 @@ subroutine parseFile(sectionNames,part,line, & endif enddo -end subroutine parseFile +end subroutine parse_materialConfig + + +!-------------------------------------------------------------------------------------------------- +!> @brief parses the material.config file +!-------------------------------------------------------------------------------------------------- +subroutine parse_debugAndNumericsConfig(config_list, & + fileContent) + implicit none + type(tPartitionedStringList), intent(out) :: config_list + character(len=pStringLen), dimension(:), intent(in) :: fileContent + integer :: i + + do i = 1, size(fileContent) + call config_list%add(trim(adjustl(fileContent(i)))) + enddo + +end subroutine parse_debugAndNumericsConfig + +end subroutine config_init + !-------------------------------------------------------------------------------------------------- !> @brief deallocates the linked lists that store the content of the configuration files @@ -233,9 +352,15 @@ subroutine config_deallocate(what) case('material.config/texture') deallocate(config_texture) - + + case('debug.config') + call config_debug%free + + case('numerics.config') + call config_numerics%free + case default - call IO_error(0_pInt,ext_msg='config_deallocate') + call IO_error(0,ext_msg='config_deallocate') end select @@ -375,7 +500,7 @@ end function keyExists !> @brief count number of key appearances !> @details traverses list and counts each occurrence of specified key !-------------------------------------------------------------------------------------------------- -integer(pInt) function countKeys(this,key) +integer function countKeys(this,key) use IO, only: & IO_stringValue @@ -385,12 +510,12 @@ integer(pInt) function countKeys(this,key) character(len=*), intent(in) :: key type(tPartitionedStringList), pointer :: item - countKeys = 0_pInt + countKeys = 0 item => this do while (associated(item%next)) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) & - countKeys = countKeys + 1_pInt + countKeys = countKeys + 1 item => item%next enddo @@ -422,13 +547,13 @@ real(pReal) function getFloat(this,key,defaultVal) do while (associated(item%next)) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then found = .true. - if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) + if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) getFloat = IO_FloatValue(item%string%val,item%string%pos,2) endif item => item%next enddo - if (.not. found) call IO_error(140_pInt,ext_msg=key) + if (.not. found) call IO_error(140,ext_msg=key) end function getFloat @@ -438,7 +563,7 @@ end function getFloat !> @details gets the last value if the key occurs more than once. If key is not found exits with !! error unless default is given !-------------------------------------------------------------------------------------------------- -integer(pInt) function getInt(this,key,defaultVal) +integer function getInt(this,key,defaultVal) use IO, only: & IO_error, & IO_stringValue, & @@ -447,7 +572,7 @@ integer(pInt) function getInt(this,key,defaultVal) implicit none class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key - integer(pInt), intent(in), optional :: defaultVal + integer, intent(in), optional :: defaultVal type(tPartitionedStringList), pointer :: item logical :: found @@ -458,13 +583,13 @@ integer(pInt) function getInt(this,key,defaultVal) do while (associated(item%next)) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then found = .true. - if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) + if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) getInt = IO_IntValue(item%string%val,item%string%pos,2) endif item => item%next enddo - if (.not. found) call IO_error(140_pInt,ext_msg=key) + if (.not. found) call IO_error(140,ext_msg=key) end function getInt @@ -497,14 +622,14 @@ character(len=65536) function getString(this,key,defaultVal,raw) found = present(defaultVal) if (found) then getString = trim(defaultVal) - if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0_pInt,ext_msg='getString') + if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0,ext_msg='getString') endif item => this do while (associated(item%next)) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then found = .true. - if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) + if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) if (whole) then getString = trim(item%string%val(item%string%pos(4):)) ! raw string starting a second chunk @@ -515,7 +640,7 @@ character(len=65536) function getString(this,key,defaultVal,raw) item => item%next enddo - if (.not. found) call IO_error(140_pInt,ext_msg=key) + if (.not. found) call IO_error(140,ext_msg=key) end function getString @@ -536,9 +661,9 @@ function getFloats(this,key,defaultVal,requiredSize) class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key real(pReal), dimension(:), intent(in), optional :: defaultVal - integer(pInt), intent(in), optional :: requiredSize + integer, intent(in), optional :: requiredSize type(tPartitionedStringList), pointer :: item - integer(pInt) :: i + integer :: i logical :: found, & cumulative @@ -552,8 +677,8 @@ function getFloats(this,key,defaultVal,requiredSize) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then found = .true. if (.not. cumulative) getFloats = [real(pReal)::] - if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) - do i = 2_pInt, item%string%pos(1) + if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) + do i = 2, item%string%pos(1) getFloats = [getFloats,IO_FloatValue(item%string%val,item%string%pos,i)] enddo endif @@ -561,7 +686,7 @@ function getFloats(this,key,defaultVal,requiredSize) enddo if (.not. found) then - if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif + if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140,ext_msg=key); endif endif if (present(requiredSize)) then if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key) @@ -582,13 +707,13 @@ function getInts(this,key,defaultVal,requiredSize) IO_IntValue implicit none - integer(pInt), dimension(:), allocatable :: getInts + integer, dimension(:), allocatable :: getInts class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key - integer(pInt), dimension(:), intent(in), optional :: defaultVal - integer(pInt), intent(in), optional :: requiredSize + integer, dimension(:), intent(in), optional :: defaultVal + integer, intent(in), optional :: requiredSize type(tPartitionedStringList), pointer :: item - integer(pInt) :: i + integer :: i logical :: found, & cumulative @@ -601,9 +726,9 @@ function getInts(this,key,defaultVal,requiredSize) do while (associated(item%next)) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then found = .true. - if (.not. cumulative) getInts = [integer(pInt)::] - if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) - do i = 2_pInt, item%string%pos(1) + if (.not. cumulative) getInts = [integer::] + if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) + do i = 2, item%string%pos(1) getInts = [getInts,IO_IntValue(item%string%val,item%string%pos,i)] enddo endif @@ -611,7 +736,7 @@ function getInts(this,key,defaultVal,requiredSize) enddo if (.not. found) then - if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif + if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140,ext_msg=key); endif endif if (present(requiredSize)) then if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key) @@ -639,7 +764,7 @@ function getStrings(this,key,defaultVal,raw) logical, intent(in), optional :: raw type(tPartitionedStringList), pointer :: item character(len=65536) :: str - integer(pInt) :: i + integer :: i logical :: found, & whole, & cumulative @@ -657,16 +782,16 @@ function getStrings(this,key,defaultVal,raw) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then found = .true. if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings) - if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key) + if (item%string%pos(1) < 2) call IO_error(143,ext_msg=key) notAllocated: if (.not. allocated(getStrings)) then if (whole) then str = item%string%val(item%string%pos(4):) getStrings = [str] else - str = IO_StringValue(item%string%val,item%string%pos,2_pInt) + str = IO_StringValue(item%string%val,item%string%pos,2) allocate(getStrings(1),source=str) - do i=3_pInt,item%string%pos(1) + do i=3,item%string%pos(1) str = IO_StringValue(item%string%val,item%string%pos,i) getStrings = [getStrings,str] enddo @@ -676,7 +801,7 @@ function getStrings(this,key,defaultVal,raw) str = item%string%val(item%string%pos(4):) getStrings = [getStrings,str] else - do i=2_pInt,item%string%pos(1) + do i=2,item%string%pos(1) str = IO_StringValue(item%string%val,item%string%pos,i) getStrings = [getStrings,str] enddo @@ -687,7 +812,7 @@ function getStrings(this,key,defaultVal,raw) enddo if (.not. found) then - if (present(defaultVal)) then; getStrings = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif + if (present(defaultVal)) then; getStrings = defaultVal; else; call IO_error(140,ext_msg=key); endif endif end function getStrings diff --git a/src/grid_damage_spectral.f90 b/src/grid_damage_spectral.f90 new file mode 100644 index 000000000..0663545d3 --- /dev/null +++ b/src/grid_damage_spectral.f90 @@ -0,0 +1,364 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH +!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Spectral solver for nonlocal damage +!-------------------------------------------------------------------------------------------------- +module grid_damage_spectral +#include +#include + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + +!-------------------------------------------------------------------------------------------------- +! PETSc data + SNES, private :: damage_snes + Vec, private :: solution_vec + PetscInt, private :: xstart, xend, ystart, yend, zstart, zend + real(pReal), private, dimension(:,:,:), allocatable :: & + damage_current, & !< field of current damage + damage_lastInc, & !< field of previous damage + damage_stagInc !< field of staggered damage + +!-------------------------------------------------------------------------------------------------- +! reference diffusion tensor, mobility etc. + integer, private :: totalIter = 0 !< total iteration in current increment + real(pReal), dimension(3,3), private :: D_ref + real(pReal), private :: mobility_ref + + public :: & + grid_damage_spectral_init, & + grid_damage_spectral_solution, & + grid_damage_spectral_forward + private :: & + formResidual + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data +! ToDo: Restart not implemented +!-------------------------------------------------------------------------------------------------- +subroutine grid_damage_spectral_init + use spectral_utilities, only: & + wgt + use mesh, only: & + grid, & + grid3 + use damage_nonlocal, only: & + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + + implicit none + PetscInt, dimension(worldsize) :: localK + integer :: i, j, k, cell + DM :: damage_grid + Vec :: uBound, lBound + PetscErrorCode :: ierr + character(len=100) :: snes_type + + write(6,'(/,a)') ' <<<+- grid_spectral_damage init -+>>>' + + write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' + write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' + +!-------------------------------------------------------------------------------------------------- +! set default and user defined options for PETSc + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-damage_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc + call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) + localK = 0 + localK(worldrank+1) = grid3 + call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) + call DMDACreate3D(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1, 1, worldsize, & + 1, 0, & ! #dof (damage phase field), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + damage_grid,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) ! connect snes to da + call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr) + call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) + call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments + call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) + if (trim(snes_type) == 'vinewtonrsls' .or. & + trim(snes_type) == 'vinewtonssls') then + call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) + call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) + call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) + call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) + call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc. + call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) + call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) + endif + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) + CHKERRQ(ierr) + xend = xstart + xend - 1 + yend = ystart + yend - 1 + zend = zstart + zend - 1 + allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) + allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) + allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) + call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! damage reference diffusion update + cell = 0 + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) + mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + +end subroutine grid_damage_spectral_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the spectral damage scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(solution) + use numerics, only: & + itmax, & + err_damage_tolAbs, & + err_damage_tolRel + use mesh, only: & + grid, & + grid3 + use damage_nonlocal, only: & + damage_nonlocal_putNonLocalDamage + + implicit none + 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 + integer :: i, j, k, cell + type(tSolutionState) :: solution + PetscInt ::position + PetscReal :: minDamage, maxDamage, stagNorm, solnNorm + + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + solution%converged =.false. + +!-------------------------------------------------------------------------------------------------- +! set module wide availabe data + params%timeinc = timeinc + params%timeincOld = timeinc_old + + call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) + call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) + + if (reason < 1) then + solution%converged = .false. + solution%iterationsNeeded = itmax + else + solution%converged = .true. + solution%iterationsNeeded = totalIter + endif + stagNorm = maxval(abs(damage_current - damage_stagInc)) + solnNorm = maxval(abs(damage_current)) + call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) + damage_stagInc = damage_current + solution%stagConverged = stagNorm < min(err_damage_tolAbs, err_damage_tolRel*solnNorm) + +!-------------------------------------------------------------------------------------------------- +! updating damage state + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) + enddo; enddo; enddo + + call VecMin(solution_vec,position,minDamage,ierr); CHKERRQ(ierr) + call VecMax(solution_vec,position,maxDamage,ierr); CHKERRQ(ierr) + if (solution%converged) & + write(6,'(/,a)') ' ... nonlocal damage converged .....................................' + write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& + minDamage, maxDamage, stagNorm + write(6,'(/,a)') ' ===========================================================================' + flush(6) + +end function grid_damage_spectral_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief spectral damage forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine grid_damage_spectral_forward + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + cutBack, & + wgt + use damage_nonlocal, only: & + damage_nonlocal_putNonLocalDamage, & + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + integer :: i, j, k, cell + DM :: dm_local + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr + + if (cutBack) then + damage_current = damage_lastInc + damage_stagInc = damage_lastInc +!-------------------------------------------------------------------------------------------------- +! reverting damage field state + cell = 0 + call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current + call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) + enddo; enddo; enddo + else +!-------------------------------------------------------------------------------------------------- +! update rate and forward last inc + damage_lastInc = damage_current + cell = 0 + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) + mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + endif + +end subroutine grid_damage_spectral_forward + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the spectral damage residual vector +!-------------------------------------------------------------------------------------------------- +subroutine formResidual(in,x_scal,f_scal,dummy,ierr) + use numerics, only: & + residualStiffness + use mesh, only: & + grid, & + grid3 + use math, only: & + math_mul33x3 + use spectral_utilities, only: & + scalarField_real, & + vectorField_real, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGreenConvolution, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence + use damage_nonlocal, only: & + damage_nonlocal_getSourceAndItsTangent,& + damage_nonlocal_getDiffusion33, & + damage_nonlocal_getMobility + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension( & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & + x_scal + PetscScalar, dimension( & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & + f_scal + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: i, j, k, cell + real(pReal) :: phiDot, dPhiDot_dPhi, mobility + + damage_current = x_scal +!-------------------------------------------------------------------------------------------------- +! evaluate polarization field + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current + call utilities_FFTscalarForward + call utilities_fourierScalarGradient !< calculate gradient of damage field + call utilities_FFTvectorBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & + vectorField_real(1:3,i,j,k)) + enddo; enddo; enddo + call utilities_FFTvectorForward + call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field + call utilities_FFTscalarBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell) + mobility = damage_nonlocal_getMobility(1,cell) + scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & + params%timeinc*phiDot + & + mobility*damage_lastInc(i,j,k) - & + mobility*damage_current(i,j,k) + & + mobility_ref*damage_current(i,j,k) + enddo; enddo; enddo + +!-------------------------------------------------------------------------------------------------- +! convolution of damage field with green operator + call utilities_FFTscalarForward + call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) + call utilities_FFTscalarBackward + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc + where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & + scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness + +!-------------------------------------------------------------------------------------------------- +! constructing residual + f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current + +end subroutine formResidual + +end module grid_damage_spectral diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 new file mode 100644 index 000000000..ebcc28b5e --- /dev/null +++ b/src/grid_mech_spectral_basic.f90 @@ -0,0 +1,545 @@ +!-------------------------------------------------------------------------------------------------- +!> @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 Grid solver for mechanics: Spectral basic +!-------------------------------------------------------------------------------------------------- +module grid_mech_spectral_basic +#include +#include + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use math, only: & + math_I3 + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private + + character (len=*), parameter, public :: & + GRID_MECH_SPECTRAL_BASIC_LABEL = 'basic' + +!-------------------------------------------------------------------------------------------------- +! derived types + type(tSolutionParams), private :: params + +!-------------------------------------------------------------------------------------------------- +! PETSc data + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec + +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), private, dimension(3,3) :: & + F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient + F_aim = math_I3, & !< current prescribed deformation gradient + F_aim_lastInc = math_I3, & !< previous average deformation gradient + P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress + + character(len=1024), private :: incInfo !< time and increment information + + real(pReal), private, dimension(3,3,3,3) :: & + C_volAvg = 0.0_pReal, & !< current volume average stiffness + C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness + C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness + C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness + S = 0.0_pReal !< current compliance (filled up with zeros) + + real(pReal), private :: & + err_BC, & !< deviation from stress BC + err_div !< RMS of div of P + + integer, private :: & + totalIter = 0 !< total iteration in current increment + + public :: & + grid_mech_spectral_basic_init, & + grid_mech_spectral_basic_solution, & + grid_mech_spectral_basic_forward + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all necessary fields and fills them with data, potentially from restart info +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_spectral_basic_init + use IO, only: & + IO_intOut, & + IO_error, & + IO_open_jobFile_binary + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRestart + use FEsolving, only: & + restartInc + use numerics, only: & + worldrank, & + worldsize, & + petsc_options + use homogenization, only: & + materialpoint_F0 + use DAMASK_interface, only: & + getSolverJobName + use spectral_utilities, only: & + utilities_constitutiveResponse, & + utilities_updateGamma, & + utilities_updateIPcoords, & + wgt + use mesh, only: & + grid, & + grid3 + use math, only: & + math_invSym3333 + + implicit none + real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P + real(pReal), dimension(3,3) :: & + temp33_Real = 0.0_pReal + + PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: F + PetscInt, dimension(worldsize) :: localK + integer :: fileUnit + character(len=1024) :: rankStr + + write(6,'(/,a)') ' <<<+- grid_mech_spectral_basic init -+>>>' + + write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' + write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' + + write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' + write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' + +!-------------------------------------------------------------------------------------------------- +! set default and user defined options for PETSc + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! allocate global fields + allocate (F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) + +!-------------------------------------------------------------------------------------------------- +! initialize solver specific parts of PETSc + call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) + localK = 0 + localK(worldrank+1) = grid3 + call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) + call DMDACreate3d(PETSC_COMM_WORLD, & + DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary + DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point + grid(1),grid(2),grid(3), & ! global grid + 1 , 1, worldsize, & + 9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + da,ierr) ! handle, error + CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da + call DMsetFromOptions(da,ierr); CHKERRQ(ierr) + call DMsetUp(da,ierr); CHKERRQ(ierr) + call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) + call DMDASNESsetFunctionLocal(da,INSERT_VALUES,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr)! residual vector of same shape as solution vector + CHKERRQ(ierr) + call SNESsetConvergenceTest(snes,grid_mech_spectral_basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "_converged" + CHKERRQ(ierr) + call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments + +!-------------------------------------------------------------------------------------------------- +! init fields + call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with + + restart: if (restartInc > 0) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading values of increment ', restartInc, ' from file' + flush(6) + endif + + fileUnit = IO_open_jobFile_binary('F_aimDot') + read(fileUnit) F_aimDot; close(fileUnit) + + write(rankStr,'(a1,i0)')'_',worldrank + + fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) + read(fileUnit) F; close (fileUnit) + fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) + read(fileUnit) F_lastInc; close (fileUnit) + + F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F + call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0) call IO_error(894, ext_msg='F_aim') + F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc + call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + if(ierr /=0) call IO_error(894, ext_msg='F_aim_lastInc') + elseif (restartInc == 0) then restart + F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity + F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) + endif restart + + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) + call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + reshape(F,shape(F_lastInc)), & ! target F + 0.0_pReal, & ! time increment + math_I3) ! no rotation of boundary condition + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc + ! QUESTION: why not writing back right after reading (l.189)? + + restartRead: if (restartInc > 0) then + if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) & + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & + 'reading more values of increment ', restartInc, ' from file' + flush(6) + fileUnit = IO_open_jobFile_binary('C_volAvg') + read(fileUnit) C_volAvg; close(fileUnit) + fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') + read(fileUnit) C_volAvgLastInc; close(fileUnit) + fileUnit = IO_open_jobFile_binary('C_ref') + read(fileUnit) C_minMaxAvg; close(fileUnit) + endif restartRead + + call Utilities_updateGamma(C_minMaxAvg,.true.) + +end subroutine grid_mech_spectral_basic_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief solution for the Basic scheme with internal iterations +!-------------------------------------------------------------------------------------------------- +function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) + use numerics, only: & + update_gamma + use spectral_utilities, only: & + tBoundaryCondition, & + utilities_maskedCompliance, & + utilities_updateGamma + use FEsolving, only: & + restartWrite, & + terminallyIll + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + character(len=*), intent(in) :: & + incInfoIn + real(pReal), intent(in) :: & + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment + type(tBoundaryCondition), intent(in) :: & + stress_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + type(tSolutionState) :: & + solution + +!-------------------------------------------------------------------------------------------------- +! PETSc Data + PetscErrorCode :: ierr + SNESConvergedReason :: reason + + incInfo = incInfoIn + +!-------------------------------------------------------------------------------------------------- +! update stiffness (and gamma operator) + S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) + +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + +!-------------------------------------------------------------------------------------------------- +! solve BVP + call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) + +!-------------------------------------------------------------------------------------------------- +! check convergence + call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) + + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll + terminallyIll = .false. + + +end function grid_mech_spectral_basic_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forms the basic residual vector +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_spectral_basic_formResidual(in, F, & + residuum, dummy, ierr) + use numerics, only: & + itmax, & + itmin + use mesh, only: & + grid, & + grid3 + use math, only: & + math_rotate_backward33, & + math_mul3333xx33 + use debug, only: & + debug_level, & + debug_spectral, & + debug_spectralRotation + use spectral_utilities, only: & + tensorField_real, & + utilities_FFTtensorForward, & + utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & + utilities_constitutiveResponse, & + utilities_divergenceRMS + use IO, only: & + IO_intOut + use FEsolving, only: & + terminallyIll + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) + PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), & + intent(in) :: F !< deformation gradient field + PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), & + intent(out) :: residuum !< residuum field + real(pReal), dimension(3,3) :: & + deltaF_aim + PetscInt :: & + PETScIter, & + nfuncs + PetscObject :: dummy + PetscErrorCode :: ierr + + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + + if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment +!-------------------------------------------------------------------------------------------------- +! begin of new iteration + newIteration: if (totalIter <= PETScIter) then + totalIter = totalIter + 1 + write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & + trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax + if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) + write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & + ' deformation gradient aim =', transpose(F_aim) + flush(6) + endif newIteration + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) + P_av,C_volAvg,C_minMaxAvg, & + F,params%timeinc,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) + F_aim = F_aim - deltaF_aim + err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + +!-------------------------------------------------------------------------------------------------- +! updated deformation gradient using fix point algorithm of basic scheme + tensorField_real = 0.0_pReal + tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform + call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" + err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use + call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg + call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier + +!-------------------------------------------------------------------------------------------------- +! constructing residual + residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + +end subroutine grid_mech_spectral_basic_formResidual + + +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_spectral_basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt :: PETScIter + PetscReal :: & + xnorm, & ! not used + snorm, & ! not used + fnorm ! not used + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + divTol, & + BCTol + + divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) + + converged: if ((totalIter >= itmin .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then converged + reason = -1 + else converged + reason = 0 + endif converged + +!-------------------------------------------------------------------------------------------------- +! 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_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) + +end subroutine grid_mech_spectral_basic_converged + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!> @details find new boundary conditions and best F estimate for end of current timestep +!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates +!-------------------------------------------------------------------------------------------------- +subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use numerics, only: & + worldrank + use homogenization, only: & + materialpoint_F0 + use mesh, only: & + grid, & + grid3 + use CPFEM2, only: & + CPFEM_age + use spectral_utilities, only: & + utilities_calculateRate, & + utilities_forwardField, & + utilities_updateIPcoords, & + tBoundaryCondition, & + cutBack + use IO, only: & + IO_open_jobFile_binary + use FEsolving, only: & + restartWrite + + implicit none + logical, intent(in) :: & + guess + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + stress_BC, & + deformation_BC + real(pReal), dimension(3,3), intent(in) :: & + rotation_BC + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: F + + integer :: fileUnit + character(len=32) :: rankStr + + call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + + if (cutBack) then + C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? + C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? + else + !-------------------------------------------------------------------------------------------------- + ! restart information for spectral solver + if (restartWrite) then ! QUESTION: where is this logical properly set? + write(6,'(/,a)') ' writing converged results for restart' + flush(6) + + if (worldrank == 0) then + fileUnit = IO_open_jobFile_binary('C_volAvg','w') + write(fileUnit) C_volAvg; close(fileUnit) + fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') + write(fileUnit) C_volAvgLastInc; close(fileUnit) + fileUnit = IO_open_jobFile_binary('F_aimDot','w') + write(fileUnit) F_aimDot; close(fileUnit) + endif + + write(rankStr,'(a1,i0)')'_',worldrank + fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') + write(fileUnit) F; close (fileUnit) + fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') + write(fileUnit) F_lastInc; close (fileUnit) + endif + + call CPFEM_age ! age state and kinematics + call utilities_updateIPcoords(F) + + C_volAvgLastInc = C_volAvg + C_minMaxAvgLastInc = C_minMaxAvg + + F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aim_lastInc = F_aim + + !-------------------------------------------------------------------------------------------------- + ! calculate rate for aim + if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) + elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * deformation_BC%values + elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed + F_aimDot = & + F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + endif + + + Fdot = Utilities_calculateRate(guess, & + F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & + math_rotate_backward33(F_aimDot,rotation_BC)) + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + endif + +!-------------------------------------------------------------------------------------------------- +! update average and local deformation gradients + F_aim = F_aim_lastInc + F_aimDot * timeinc + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + +end subroutine grid_mech_spectral_basic_forward + +end module grid_mech_spectral_basic diff --git a/src/spectral_mech_Polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 similarity index 93% rename from src/spectral_mech_Polarisation.f90 rename to src/grid_mech_spectral_polarisation.f90 index d6e297c72..4746670d5 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Polarisation scheme solver !-------------------------------------------------------------------------------------------------- -module spectral_mech_Polarisation +module grid_mech_spectral_polarisation #include #include use PETScdmda @@ -22,7 +22,7 @@ module spectral_mech_Polarisation private character (len=*), parameter, public :: & - DAMASK_spectral_solverPolarisation_label = 'polarisation' + GRID_MECH_SPECTRAL_POLARISATION_LABEL = 'polarisation' !-------------------------------------------------------------------------------------------------- ! derived types @@ -70,16 +70,16 @@ module spectral_mech_Polarisation totalIter = 0_pInt !< total iteration in current increment public :: & - Polarisation_init, & - Polarisation_solution, & - Polarisation_forward + grid_mech_spectral_polarisation_init, & + grid_mech_spectral_polarisation_solution, & + grid_mech_spectral_polarisation_forward contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- -subroutine Polarisation_init +subroutine grid_mech_spectral_polarisation_init use IO, only: & IO_intOut, & IO_error, & @@ -92,7 +92,8 @@ subroutine Polarisation_init restartInc use numerics, only: & worldrank, & - worldsize + worldsize, & + petsc_options use homogenization, only: & materialpoint_F0 use DAMASK_interface, only: & @@ -127,6 +128,13 @@ subroutine Polarisation_init write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' +!-------------------------------------------------------------------------------------------------- +! set default and user defined options for PETSc + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) + !-------------------------------------------------------------------------------------------------- ! allocate global fields allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) @@ -154,9 +162,9 @@ subroutine Polarisation_init call DMsetFromOptions(da,ierr); CHKERRQ(ierr) call DMsetUp(da,ierr); CHKERRQ(ierr) call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) - call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector + call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,Polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + call SNESsetConvergenceTest(snes,grid_mech_spectral_polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -227,13 +235,13 @@ subroutine Polarisation_init C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) -end subroutine Polarisation_init +end subroutine grid_mech_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) +function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) use IO, only: & IO_error use numerics, only: & @@ -255,12 +263,13 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol character(len=*), intent(in) :: & incInfoIn real(pReal), intent(in) :: & - timeinc, & !< increment time for current solution - timeinc_old !< increment time of last successful increment + timeinc, & !< time increment of current solution + timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC - + type(tSolutionState) :: & + solution !-------------------------------------------------------------------------------------------------- ! PETSc Data PetscErrorCode :: ierr @@ -293,19 +302,19 @@ type(tSolutionState) function Polarisation_solution(incInfoIn,timeinc,timeinc_ol ! check convergence call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - Polarisation_solution%converged = reason > 0 - Polarisation_solution%iterationsNeeded = totalIter - Polarisation_solution%termIll = terminallyIll + solution%converged = reason > 0 + solution%iterationsNeeded = totalIter + solution%termIll = terminallyIll terminallyIll = .false. if (reason == -4) call IO_error(893_pInt) ! MPI error -end function Polarisation_solution +end function grid_mech_spectral_polarisation_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the Polarisation residual vector !-------------------------------------------------------------------------------------------------- -subroutine Polarisation_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) +subroutine formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) FandF_tau, & ! defgrad fields on grid residuum, & ! residuum fields on grid dummy, & @@ -449,13 +458,13 @@ subroutine Polarisation_formResidual(in, & nullify(F_tau) nullify(residual_F) nullify(residual_F_tau) -end subroutine Polarisation_formResidual +end subroutine formResidual !-------------------------------------------------------------------------------------------------- !> @brief convergence check !-------------------------------------------------------------------------------------------------- -subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) +subroutine grid_mech_spectral_polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) use numerics, only: & itmax, & itmin, & @@ -521,14 +530,14 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, write(6,'(/,a)') ' ===========================================================================' flush(6) -end subroutine Polarisation_converged +end subroutine grid_mech_spectral_polarisation_converged !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !> @details find new boundary conditions and best F estimate for end of current timestep !> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates !-------------------------------------------------------------------------------------------------- -subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) +subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) use math, only: & math_mul33x33, & math_mul3333xx33, & @@ -670,6 +679,6 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformati nullify(F_tau) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) -end subroutine Polarisation_forward +end subroutine grid_mech_spectral_polarisation_forward -end module spectral_mech_Polarisation +end module grid_mech_spectral_polarisation diff --git a/src/spectral_thermal.f90 b/src/grid_thermal_spectral.f90 similarity index 53% rename from src/spectral_thermal.f90 rename to src/grid_thermal_spectral.f90 index 3e2a4b1f9..69e23c86b 100644 --- a/src/spectral_thermal.f90 +++ b/src/grid_thermal_spectral.f90 @@ -1,56 +1,56 @@ !-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH !> @brief Spectral solver for thermal conduction !-------------------------------------------------------------------------------------------------- -module spectral_thermal +module grid_thermal_spectral #include #include - use PETScdmda - use PETScsnes - use prec, only: & - pReal - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private - - character (len=*), parameter, public :: & - spectral_thermal_label = 'spectralthermal' - + use PETScdmda + use PETScsnes + use prec, only: & + pReal + use spectral_utilities, only: & + tSolutionState, & + tSolutionParams + + implicit none + private !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - SNES, private :: thermal_snes - Vec, private :: solution - PetscInt, private :: xstart, xend, ystart, yend, zstart, zend - real(pReal), private, dimension(:,:,:), allocatable :: & - temperature_current, & !< field of current temperature - temperature_lastInc, & !< field of previous temperature - temperature_stagInc !< field of staggered temperature + SNES, private :: thermal_snes + Vec, private :: solution_vec + PetscInt, private :: xstart, xend, ystart, yend, zstart, zend + real(pReal), private, dimension(:,:,:), allocatable :: & + temperature_current, & !< field of current temperature + temperature_lastInc, & !< field of previous temperature + temperature_stagInc !< field of staggered temperature !-------------------------------------------------------------------------------------------------- ! reference diffusion tensor, mobility etc. - integer, private :: totalIter = 0 !< total iteration in current increment - real(pReal), dimension(3,3), private :: D_ref - real(pReal), private :: mobility_ref - - public :: & - spectral_thermal_init, & - spectral_thermal_solution, & - spectral_thermal_forward + integer, private :: totalIter = 0 !< total iteration in current increment + real(pReal), dimension(3,3), private :: D_ref + real(pReal), private :: mobility_ref + + public :: & + grid_thermal_spectral_init, & + grid_thermal_spectral_solution, & + grid_thermal_spectral_forward + private :: & + formResidual contains !-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!> @brief allocates all neccessary fields and fills them with data +! ToDo: Restart not implemented !-------------------------------------------------------------------------------------------------- -subroutine spectral_thermal_init +subroutine grid_thermal_spectral_init use spectral_utilities, only: & wgt use mesh, only: & @@ -66,19 +66,27 @@ subroutine spectral_thermal_init thermalMapping use numerics, only: & worldrank, & - worldsize + worldsize, & + petsc_options implicit none - integer, dimension(worldsize) :: localK + PetscInt, dimension(worldsize) :: localK integer :: i, j, k, cell DM :: thermal_grid - PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr - write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' + write(6,'(/,a)') ' <<<+- grid_thermal_spectral init -+>>>' write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' + +!-------------------------------------------------------------------------------------------------- +! set default and user defined options for PETSc + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-thermal_snes_type ngmres',ierr) + CHKERRQ(ierr) + call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc @@ -92,16 +100,15 @@ subroutine spectral_thermal_init DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid 1, 1, worldsize, & - 1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & !< local grid - thermal_grid,ierr) !< handle, error + 1, 0, & ! #dof (thermal phase field), ghost boundary width (domain overlap) + [grid(1)],[grid(2)],localK, & ! local grid + thermal_grid,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr) call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) - call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,& - PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector + call DMCreateGlobalVector(thermal_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) + call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -116,16 +123,16 @@ subroutine spectral_thermal_init allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 temperature_current(i,j,k) = temperature(material_homogenizationAt(cell))% & p(thermalMapping(material_homogenizationAt(cell))%p(1,cell)) temperature_lastInc(i,j,k) = temperature_current(i,j,k) temperature_stagInc(i,j,k) = temperature_current(i,j,k) enddo; enddo; enddo - call DMDAVecGetArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current - call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! thermal reference diffusion update @@ -143,12 +150,13 @@ subroutine spectral_thermal_init mobility_ref = mobility_ref*wgt call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) -end subroutine spectral_thermal_init +end subroutine grid_thermal_spectral_init + !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,loadCaseTime) +function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(solution) use numerics, only: & itmax, & err_thermal_tolAbs, & @@ -160,42 +168,41 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load thermal_conduction_putTemperatureAndItsRate implicit none - 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 + timeinc, & !< increment in time for current solution + timeinc_old, & !< increment in time of last increment + loadCaseTime !< remaining time of current load case integer :: i, j, k, cell + type(tSolutionState) :: solution PetscInt :: position PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm PetscErrorCode :: ierr SNESConvergedReason :: reason - spectral_thermal_solution%converged =.false. + solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc params%timeincOld = timeinc_old - call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) + call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) if (reason < 1) then - spectral_thermal_solution%converged = .false. - spectral_thermal_solution%iterationsNeeded = itmax + solution%converged = .false. + solution%iterationsNeeded = itmax else - spectral_thermal_solution%converged = .true. - spectral_thermal_solution%iterationsNeeded = totalIter + solution%converged = .true. + solution%iterationsNeeded = totalIter endif stagNorm = maxval(abs(temperature_current - temperature_stagInc)) solnNorm = maxval(abs(temperature_current)) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) temperature_stagInc = temperature_current - spectral_thermal_solution%stagConverged = stagNorm < err_thermal_tolAbs & - .or. stagNorm < err_thermal_tolRel*solnNorm + solution%stagConverged = stagNorm < min(err_thermal_tolAbs, err_thermal_tolRel*solnNorm) !-------------------------------------------------------------------------------------------------- ! updating thermal state @@ -207,157 +214,158 @@ type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,load 1,cell) enddo; enddo; enddo - call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr) - call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr) - if (spectral_thermal_solution%converged) & + call VecMin(solution_vec,position,minTemperature,ierr); CHKERRQ(ierr) + call VecMax(solution_vec,position,maxTemperature,ierr); CHKERRQ(ierr) + if (solution%converged) & write(6,'(/,a)') ' ... thermal conduction converged ..................................' write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& minTemperature, maxTemperature, stagNorm write(6,'(/,a)') ' ===========================================================================' flush(6) -end function spectral_thermal_solution +end function grid_thermal_spectral_solution + + +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine grid_thermal_spectral_forward + use mesh, only: & + grid, & + grid3 + use spectral_utilities, only: & + cutBack, & + wgt + use thermal_conduction, only: & + thermal_conduction_putTemperatureAndItsRate, & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + + implicit none + integer :: i, j, k, cell + DM :: dm_local + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr + + if (cutBack) then + temperature_current = temperature_lastInc + temperature_stagInc = temperature_lastInc + +!-------------------------------------------------------------------------------------------------- +! reverting thermal field state + cell = 0 + call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with + x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current + call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & + (temperature_current(i,j,k) - & + temperature_lastInc(i,j,k))/params%timeinc, & + 1,cell) + enddo; enddo; enddo + else +!-------------------------------------------------------------------------------------------------- +! update rate and forward last inc + temperature_lastInc = temperature_current + cell = 0 + D_ref = 0.0_pReal + mobility_ref = 0.0_pReal + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) + mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & + thermal_conduction_getSpecificHeat(1,cell) + enddo; enddo; enddo + D_ref = D_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + mobility_ref = mobility_ref*wgt + call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + endif + +end subroutine grid_thermal_spectral_forward !-------------------------------------------------------------------------------------------------- !> @brief forms the spectral thermal residual vector !-------------------------------------------------------------------------------------------------- -subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) - use mesh, only: & - grid, & - grid3 - use math, only: & - math_mul33x3 - use spectral_utilities, only: & - scalarField_real, & - vectorField_real, & - utilities_FFTvectorForward, & - utilities_FFTvectorBackward, & - utilities_FFTscalarForward, & - utilities_FFTscalarBackward, & - utilities_fourierGreenConvolution, & - utilities_fourierScalarGradient, & - utilities_fourierVectorDivergence - use thermal_conduction, only: & - thermal_conduction_getSourceAndItsTangent, & - thermal_conduction_getConductivity33, & - thermal_conduction_getMassDensity, & - thermal_conduction_getSpecificHeat +subroutine formResidual(in,x_scal,f_scal,dummy,ierr) + use mesh, only: & + grid, & + grid3 + use math, only: & + math_mul33x3 + use spectral_utilities, only: & + scalarField_real, & + vectorField_real, & + utilities_FFTvectorForward, & + utilities_FFTvectorBackward, & + utilities_FFTscalarForward, & + utilities_FFTscalarBackward, & + utilities_fourierGreenConvolution, & + utilities_fourierScalarGradient, & + utilities_fourierVectorDivergence + use thermal_conduction, only: & + thermal_conduction_getSourceAndItsTangent, & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + + implicit none + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & + in + PetscScalar, dimension( & + XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & + x_scal + PetscScalar, dimension( & + X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & + f_scal + PetscObject :: dummy + PetscErrorCode :: ierr + integer :: i, j, k, cell + real(pReal) :: Tdot, dTdot_dT - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & - in - PetscScalar, dimension( & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, dimension( & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal - PetscObject :: dummy - PetscErrorCode :: ierr - integer :: i, j, k, cell - real(pReal) :: Tdot, dTdot_dT - - temperature_current = x_scal + temperature_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current - call utilities_FFTscalarForward() - call utilities_fourierScalarGradient() !< calculate gradient of damage field - call utilities_FFTvectorBackward() - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & - vectorField_real(1:3,i,j,k)) - enddo; enddo; enddo - call utilities_FFTvectorForward() - call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field - call utilities_FFTscalarBackward() - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) - scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & - params%timeinc*Tdot + & - thermal_conduction_getMassDensity (1,cell)* & - thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - & - temperature_current(i,j,k)) + & - mobility_ref*temperature_current(i,j,k) - enddo; enddo; enddo + scalarField_real = 0.0_pReal + scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current + call utilities_FFTscalarForward + call utilities_fourierScalarGradient !< calculate gradient of damage field + call utilities_FFTvectorBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & + vectorField_real(1:3,i,j,k)) + enddo; enddo; enddo + call utilities_FFTvectorForward + call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field + call utilities_FFTscalarBackward + cell = 0 + do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) + cell = cell + 1 + call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) + scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & + params%timeinc*Tdot + & + thermal_conduction_getMassDensity (1,cell)* & + thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - & + temperature_current(i,j,k)) + & + mobility_ref*temperature_current(i,j,k) + enddo; enddo; enddo !-------------------------------------------------------------------------------------------------- ! convolution of damage field with green operator - call utilities_FFTscalarForward() - call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) - call utilities_FFTscalarBackward() + call utilities_FFTscalarForward + call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) + call utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) + f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3) -end subroutine spectral_thermal_formResidual +end subroutine formResidual -!-------------------------------------------------------------------------------------------------- -!> @brief forwarding routine -!-------------------------------------------------------------------------------------------------- -subroutine spectral_thermal_forward() - use mesh, only: & - grid, & - grid3 - use spectral_utilities, only: & - cutBack, & - wgt - use thermal_conduction, only: & - thermal_conduction_putTemperatureAndItsRate, & - thermal_conduction_getConductivity33, & - thermal_conduction_getMassDensity, & - thermal_conduction_getSpecificHeat - - implicit none - integer :: i, j, k, cell - DM :: dm_local - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr - - if (cutBack) then - temperature_current = temperature_lastInc - temperature_stagInc = temperature_lastInc - -!-------------------------------------------------------------------------------------------------- -! reverting thermal field state - cell = 0 - call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with - x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current - call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & - (temperature_current(i,j,k) - & - temperature_lastInc(i,j,k))/params%timeinc, & - 1,cell) - enddo; enddo; enddo - else -!-------------------------------------------------------------------------------------------------- -! update rate and forward last inc - temperature_lastInc = temperature_current - cell = 0 - D_ref = 0.0_pReal - mobility_ref = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) - mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & - thermal_conduction_getSpecificHeat(1,cell) - enddo; enddo; enddo - D_ref = D_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - mobility_ref = mobility_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - endif - -end subroutine spectral_thermal_forward - -end module spectral_thermal +end module grid_thermal_spectral diff --git a/src/math.f90 b/src/math.f90 index 660f76190..f77f81f8b 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -7,8 +7,7 @@ !-------------------------------------------------------------------------------------------------- module math use prec, only: & - pReal, & - pInt + pReal implicit none private @@ -34,37 +33,37 @@ module math 1.0_pReal, 1.0_pReal, 1.0_pReal, & 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward) - integer(pInt), dimension (2,6), parameter, private :: & + integer, dimension (2,6), parameter, private :: & mapNye = reshape([& - 1_pInt,1_pInt, & - 2_pInt,2_pInt, & - 3_pInt,3_pInt, & - 1_pInt,2_pInt, & - 2_pInt,3_pInt, & - 1_pInt,3_pInt & + 1,1, & + 2,2, & + 3,3, & + 1,2, & + 2,3, & + 1,3 & ],[2,6]) !< arrangement in Nye notation. - integer(pInt), dimension (2,6), parameter, private :: & + integer, dimension (2,6), parameter, private :: & mapVoigt = reshape([& - 1_pInt,1_pInt, & - 2_pInt,2_pInt, & - 3_pInt,3_pInt, & - 2_pInt,3_pInt, & - 1_pInt,3_pInt, & - 1_pInt,2_pInt & + 1,1, & + 2,2, & + 3,3, & + 2,3, & + 1,3, & + 1,2 & ],[2,6]) !< arrangement in Voigt notation - integer(pInt), dimension (2,9), parameter, private :: & + integer, dimension (2,9), parameter, private :: & mapPlain = reshape([& - 1_pInt,1_pInt, & - 1_pInt,2_pInt, & - 1_pInt,3_pInt, & - 2_pInt,1_pInt, & - 2_pInt,2_pInt, & - 2_pInt,3_pInt, & - 3_pInt,1_pInt, & - 3_pInt,2_pInt, & - 3_pInt,3_pInt & + 1,1, & + 1,2, & + 1,3, & + 2,1, & + 2,2, & + 2,3, & + 3,1, & + 3,2, & + 3,3 & ],[2,9]) !< arrangement in Plain notation !-------------------------------------------------------------------------------------------------- @@ -184,18 +183,17 @@ subroutine math_init randomSeed implicit none - integer(pInt) :: i - real(pReal), dimension(4) :: randTest + integer :: i + real(pReal), dimension(4) :: randTest integer :: randSize integer, dimension(:), allocatable :: randInit write(6,'(/,a)') ' <<<+- math init -+>>>' call random_seed(size=randSize) - if (allocated(randInit)) deallocate(randInit) allocate(randInit(randSize)) - if (randomSeed > 0_pInt) then - randInit(1:randSize) = int(randomSeed) ! randomSeed is of type pInt, randInit not + if (randomSeed > 0) then + randInit = randomSeed call random_seed(put=randInit) else call random_seed() @@ -204,12 +202,12 @@ subroutine math_init call random_seed(put = randInit) endif - do i = 1_pInt, 4_pInt + do i = 1, 4 call random_number(randTest(i)) enddo write(6,'(a,I2)') ' size of random seed: ', randSize - do i = 1_pInt,randSize + do i = 1,randSize write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i) enddo write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest @@ -244,7 +242,7 @@ subroutine math_check any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) & 'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif ! +++ q -> R -> q +++ @@ -254,7 +252,7 @@ subroutine math_check any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) & 'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif ! +++ q -> euler -> q +++ @@ -264,7 +262,7 @@ subroutine math_check any(abs(-q-q2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) & 'quat -> euler -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif ! +++ R -> euler -> R +++ @@ -273,32 +271,32 @@ subroutine math_check if ( any(abs( R-R2) > tol_math_check) ) then write (error_msg, '(a,e14.6)' ) & 'R -> euler -> R maximum deviation ',maxval(abs( R-R2)) - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif ! +++ check rotation sense of q and R +++ - v = halton([2_pInt,8_pInt,5_pInt]) ! random vector + v = halton([2,8,5]) ! random vector R = math_qToR(q) if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v' - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif ! +++ check vector expansion +++ if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & - math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt,3_pInt,0_pInt])) > tol_math_check)) then + math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) then write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - & - math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt])) > tol_math_check)) then + math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2])) > tol_math_check)) then write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]' - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - & - math_expand([1.0_pReal,2.0_pReal],[1_pInt,2_pInt,3_pInt])) > tol_math_check)) then + math_expand([1.0_pReal,2.0_pReal],[1,2,3])) > tol_math_check)) then write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' - call IO_error(401_pInt,ext_msg=error_msg) + call IO_error(401,ext_msg=error_msg) endif end subroutine math_check @@ -312,9 +310,9 @@ end subroutine math_check recursive subroutine math_qsort(a, istart, iend, sortDim) implicit none - integer(pInt), dimension(:,:), intent(inout) :: a - integer(pInt), intent(in),optional :: istart,iend, sortDim - integer(pInt) :: ipivot,s,e,d + integer, dimension(:,:), intent(inout) :: a + integer, intent(in),optional :: istart,iend, sortDim + integer :: ipivot,s,e,d if(present(istart)) then s = istart @@ -336,8 +334,8 @@ recursive subroutine math_qsort(a, istart, iend, sortDim) if (s < e) then ipivot = qsort_partition(a,s, e, d) - call math_qsort(a, s, ipivot-1_pInt, d) - call math_qsort(a, ipivot+1_pInt, e, d) + call math_qsort(a, s, ipivot-1, d) + call math_qsort(a, ipivot+1, e, d) endif !-------------------------------------------------------------------------------------------------- @@ -346,17 +344,17 @@ recursive subroutine math_qsort(a, istart, iend, sortDim) !------------------------------------------------------------------------------------------------- !> @brief Partitioning required for quicksort !------------------------------------------------------------------------------------------------- - integer(pInt) function qsort_partition(a, istart, iend, sort) + integer function qsort_partition(a, istart, iend, sort) implicit none - integer(pInt), dimension(:,:), intent(inout) :: a - integer(pInt), intent(in) :: istart,iend,sort - integer(pInt), dimension(size(a,1)) :: tmp - integer(pInt) :: i,j + integer, dimension(:,:), intent(inout) :: a + integer, intent(in) :: istart,iend,sort + integer, dimension(size(a,1)) :: tmp + integer :: i,j do ! find the first element on the right side less than or equal to the pivot point - do j = iend, istart, -1_pInt + do j = iend, istart, -1 if (a(sort,j) <= a(sort,istart)) exit enddo ! find the first element on the left side greater than the pivot point @@ -390,15 +388,15 @@ pure function math_expand(what,how) implicit none real(pReal), dimension(:), intent(in) :: what - integer(pInt), dimension(:), intent(in) :: how + integer, dimension(:), intent(in) :: how real(pReal), dimension(sum(how)) :: math_expand - integer(pInt) :: i + integer :: i - if (sum(how) == 0_pInt) & + if (sum(how) == 0) & return - do i = 1_pInt, size(how) - math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1_pInt,size(what))+1_pInt) + do i = 1, size(how) + math_expand(sum(how(1:i-1))+1:sum(how(1:i))) = what(mod(i-1,size(what))+1) enddo end function math_expand @@ -410,11 +408,11 @@ end function math_expand pure function math_range(N) implicit none - integer(pInt), intent(in) :: N !< length of range - integer(pInt) :: i - integer(pInt), dimension(N) :: math_range + integer, intent(in) :: N !< length of range + integer :: i + integer, dimension(N) :: math_range - math_range = [(i,i=1_pInt,N)] + math_range = [(i,i=1,N)] end function math_range @@ -425,12 +423,12 @@ end function math_range pure function math_identity2nd(dimen) implicit none - integer(pInt), intent(in) :: dimen !< tensor dimension - integer(pInt) :: i + integer, intent(in) :: dimen !< tensor dimension + integer :: i real(pReal), dimension(dimen,dimen) :: math_identity2nd math_identity2nd = 0.0_pReal - forall(i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal + forall(i=1:dimen) math_identity2nd(i,i) = 1.0_pReal end function math_identity2nd @@ -441,13 +439,13 @@ end function math_identity2nd pure function math_identity4th(dimen) implicit none - integer(pInt), intent(in) :: dimen !< tensor dimension - integer(pInt) :: i,j,k,l + integer, intent(in) :: dimen !< tensor dimension + integer :: i,j,k,l real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th real(pReal), dimension(dimen,dimen) :: identity2nd identity2nd = math_identity2nd(dimen) - forall(i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt:dimen) & + forall(i=1:dimen,j=1:dimen,k=1:dimen,l=1:dimen) & math_identity4th(i,j,k,l) = 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) end function math_identity4th @@ -462,15 +460,15 @@ end function math_identity4th real(pReal) pure function math_civita(i,j,k) implicit none - integer(pInt), intent(in) :: i,j,k + integer, intent(in) :: i,j,k math_civita = 0.0_pReal - if (((i == 1_pInt).and.(j == 2_pInt).and.(k == 3_pInt)) .or. & - ((i == 2_pInt).and.(j == 3_pInt).and.(k == 1_pInt)) .or. & - ((i == 3_pInt).and.(j == 1_pInt).and.(k == 2_pInt))) math_civita = 1.0_pReal - if (((i == 1_pInt).and.(j == 3_pInt).and.(k == 2_pInt)) .or. & - ((i == 2_pInt).and.(j == 1_pInt).and.(k == 3_pInt)) .or. & - ((i == 3_pInt).and.(j == 2_pInt).and.(k == 1_pInt))) math_civita = -1.0_pReal + if (((i == 1).and.(j == 2).and.(k == 3)) .or. & + ((i == 2).and.(j == 3).and.(k == 1)) .or. & + ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal + if (((i == 1).and.(j == 3).and.(k == 2)) .or. & + ((i == 2).and.(j == 1).and.(k == 3)) .or. & + ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal end function math_civita @@ -484,7 +482,7 @@ end function math_civita real(pReal) pure function math_delta(i,j) implicit none - integer(pInt), intent (in) :: i,j + integer, intent (in) :: i,j math_delta = merge(0.0_pReal, 1.0_pReal, i /= j) @@ -515,9 +513,9 @@ pure function math_outer(A,B) implicit none real(pReal), dimension(:), intent(in) :: A,B real(pReal), dimension(size(A,1),size(B,1)) :: math_outer - integer(pInt) :: i,j + integer :: i,j - forall(i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_outer(i,j) = A(i)*B(j) + forall(i=1:size(A,1),j=1:size(B,1)) math_outer(i,j) = A(i)*B(j) end function math_outer @@ -543,10 +541,10 @@ real(pReal) pure function math_mul33xx33(A,B) implicit none real(pReal), dimension(3,3), intent(in) :: A,B - integer(pInt) :: i,j + integer :: i,j real(pReal), dimension(3,3) :: C - forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) C(i,j) = A(i,j) * B(i,j) + forall(i=1:3,j=1:3) C(i,j) = A(i,j) * B(i,j) math_mul33xx33 = sum(C) end function math_mul33xx33 @@ -561,9 +559,9 @@ pure function math_mul3333xx33(A,B) real(pReal), dimension(3,3) :: math_mul3333xx33 real(pReal), dimension(3,3,3,3), intent(in) :: A real(pReal), dimension(3,3), intent(in) :: B - integer(pInt) :: i,j + integer :: i,j - forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) + forall(i = 1:3,j = 1:3) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) end function math_mul3333xx33 @@ -574,12 +572,12 @@ end function math_mul3333xx33 pure function math_mul3333xx3333(A,B) implicit none - integer(pInt) :: i,j,k,l + integer :: i,j,k,l real(pReal), dimension(3,3,3,3), intent(in) :: A real(pReal), dimension(3,3,3,3), intent(in) :: B real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 - forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt, k = 1_pInt:3_pInt, l= 1_pInt:3_pInt) & + forall(i = 1:3,j = 1:3, k = 1:3, l= 1:3) & math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) end function math_mul3333xx3333 @@ -593,9 +591,9 @@ pure function math_mul33x33(A,B) implicit none real(pReal), dimension(3,3) :: math_mul33x33 real(pReal), dimension(3,3), intent(in) :: A,B - integer(pInt) :: i,j + integer :: i,j - forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + forall(i=1:3,j=1:3) math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) end function math_mul33x33 @@ -608,9 +606,9 @@ pure function math_mul66x66(A,B) implicit none real(pReal), dimension(6,6) :: math_mul66x66 real(pReal), dimension(6,6), intent(in) :: A,B - integer(pInt) :: i,j + integer :: i,j - forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) & + forall(i=1:6,j=1:6) & math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) & + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) @@ -625,9 +623,9 @@ pure function math_mul99x99(A,B) implicit none real(pReal), dimension(9,9) :: math_mul99x99 real(pReal), dimension(9,9), intent(in) :: A,B - integer(pInt) i,j + integer i,j - forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) & + forall(i=1:9,j=1:9) & math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) & + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) & + A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j) @@ -644,9 +642,9 @@ pure function math_mul33x3(A,B) real(pReal), dimension(3) :: math_mul33x3 real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3), intent(in) :: B - integer(pInt) :: i + integer :: i - forall (i=1_pInt:3_pInt) math_mul33x3(i) = sum(A(i,1:3)*B) + forall (i=1:3) math_mul33x3(i) = sum(A(i,1:3)*B) end function math_mul33x3 @@ -660,9 +658,9 @@ pure function math_mul33x3_complex(A,B) complex(pReal), dimension(3) :: math_mul33x3_complex complex(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3), intent(in) :: B - integer(pInt) :: i + integer :: i - forall (i=1_pInt:3_pInt) math_mul33x3_complex(i) = sum(A(i,1:3)*cmplx(B,0.0_pReal,pReal)) + forall (i=1:3) math_mul33x3_complex(i) = sum(A(i,1:3)*cmplx(B,0.0_pReal,pReal)) end function math_mul33x3_complex @@ -676,9 +674,9 @@ pure function math_mul66x6(A,B) real(pReal), dimension(6) :: math_mul66x6 real(pReal), dimension(6,6), intent(in) :: A real(pReal), dimension(6), intent(in) :: B - integer(pInt) :: i + integer :: i - forall (i=1_pInt:6_pInt) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) & + forall (i=1:6) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) & + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6) end function math_mul66x6 @@ -690,8 +688,8 @@ end function math_mul66x6 pure function math_exp33(A,n) implicit none - integer(pInt) :: i - integer(pInt), intent(in), optional :: n + integer :: i + integer, intent(in), optional :: n real(pReal), dimension(3,3), intent(in) :: A real(pReal), dimension(3,3) :: B, math_exp33 real(pReal) :: invFac @@ -700,7 +698,7 @@ pure function math_exp33(A,n) invFac = 1.0_pReal ! 0! math_exp33 = B ! A^0 = eye2 - do i = 1_pInt, merge(n,5_pInt,present(n)) + do i = 1, merge(n,5,present(n)) invFac = invFac/real(i,pReal) ! invfac = 1/i! B = math_mul33x33(B,A) math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i! @@ -717,9 +715,9 @@ pure function math_transpose33(A) implicit none real(pReal),dimension(3,3) :: math_transpose33 real(pReal),dimension(3,3),intent(in) :: A - integer(pInt) :: i,j + integer :: i,j - forall(i=1_pInt:3_pInt, j=1_pInt:3_pInt) math_transpose33(i,j) = A(j,i) + forall(i=1:3, j=1:3) math_transpose33(i,j) = A(j,i) end function math_transpose33 @@ -814,10 +812,10 @@ function math_invSym3333(A) real(pReal),dimension(3,3,3,3),intent(in) :: A - integer(pInt) :: ierr - integer(pInt), dimension(6) :: ipiv6 - real(pReal), dimension(6,6) :: temp66_Real - real(pReal), dimension(6) :: work6 + integer :: ierr + integer, dimension(6) :: ipiv6 + real(pReal), dimension(6,6) :: temp66_Real + real(pReal), dimension(6) :: work6 external :: & dgetrf, & dgetri @@ -825,10 +823,10 @@ function math_invSym3333(A) temp66_real = math_sym3333to66(A) call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) - if (ierr == 0_pInt) then + if (ierr == 0) then math_invSym3333 = math_66toSym3333(temp66_real) else - call IO_error(400_pInt, ext_msg = 'math_invSym3333') + call IO_error(400, ext_msg = 'math_invSym3333') endif end function math_invSym3333 @@ -859,13 +857,13 @@ end subroutine math_invert2 subroutine math_invert(myDim,A, InvA, error) implicit none - integer(pInt), intent(in) :: myDim + integer, intent(in) :: myDim real(pReal), dimension(myDim,myDim), intent(in) :: A - integer(pInt) :: ierr - integer(pInt), dimension(myDim) :: ipiv - real(pReal), dimension(myDim) :: work + integer :: ierr + integer, dimension(myDim) :: ipiv + real(pReal), dimension(myDim) :: work real(pReal), dimension(myDim,myDim), intent(out) :: invA logical, intent(out) :: error @@ -876,7 +874,7 @@ subroutine math_invert(myDim,A, InvA, error) invA = A call dgetrf(myDim,myDim,invA,myDim,ipiv,ierr) call dgetri(myDim,InvA,myDim,ipiv,work,myDim,ierr) - error = merge(.true.,.false., ierr /= 0_pInt) + error = merge(.true.,.false., ierr /= 0) end subroutine math_invert @@ -1029,8 +1027,8 @@ real(pReal) pure function math_detSym33(m) implicit none real(pReal), dimension(3,3), intent(in) :: m - math_detSym33 = -(m(1,1)*m(2,3)**2_pInt + m(2,2)*m(1,3)**2_pInt + m(3,3)*m(1,2)**2_pInt) & - + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3) + math_detSym33 = -(m(1,1)*m(2,3)**2 + m(2,2)*m(1,3)**2 + m(3,3)*m(1,2)**2) & + + m(1,1)*m(2,2)*m(3,3) + 2.0_pReal * m(1,2)*m(1,3)*m(2,3) end function math_detSym33 @@ -1044,9 +1042,9 @@ pure function math_33to9(m33) real(pReal), dimension(9) :: math_33to9 real(pReal), dimension(3,3), intent(in) :: m33 - integer(pInt) :: i + integer :: i - forall(i=1_pInt:9_pInt) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) + forall(i=1:9) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) end function math_33to9 @@ -1060,9 +1058,9 @@ pure function math_9to33(v9) real(pReal), dimension(3,3) :: math_9to33 real(pReal), dimension(9), intent(in) :: v9 - integer(pInt) :: i + integer :: i - forall(i=1_pInt:9_pInt) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) + forall(i=1:9) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) end function math_9to33 @@ -1081,7 +1079,7 @@ pure function math_sym33to6(m33,weighted) logical, optional, intent(in) :: weighted real(pReal), dimension(6) :: w - integer(pInt) :: i + integer :: i if(present(weighted)) then w = merge(nrmMandel,1.0_pReal,weighted) @@ -1089,7 +1087,7 @@ pure function math_sym33to6(m33,weighted) w = nrmMandel endif - forall(i=1_pInt:6_pInt) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i)) + forall(i=1:6) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i)) end function math_sym33to6 @@ -1108,7 +1106,7 @@ pure function math_6toSym33(v6,weighted) logical, optional, intent(in) :: weighted real(pReal), dimension(6) :: w - integer(pInt) :: i + integer :: i if(present(weighted)) then w = merge(invnrmMandel,1.0_pReal,weighted) @@ -1116,7 +1114,7 @@ pure function math_6toSym33(v6,weighted) w = invnrmMandel endif - do i=1_pInt,6_pInt + do i=1,6 math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i) math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i) enddo @@ -1133,9 +1131,9 @@ pure function math_3333to99(m3333) real(pReal), dimension(9,9) :: math_3333to99 real(pReal), dimension(3,3,3,3), intent(in) :: m3333 - integer(pInt) :: i,j + integer :: i,j - forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) & + forall(i=1:9,j=1:9) & math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) end function math_3333to99 @@ -1150,9 +1148,9 @@ pure function math_99to3333(m99) real(pReal), dimension(3,3,3,3) :: math_99to3333 real(pReal), dimension(9,9), intent(in) :: m99 - integer(pInt) :: i,j + integer :: i,j - forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) & + forall(i=1:9,j=1:9) & math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j) end function math_99to3333 @@ -1172,7 +1170,7 @@ pure function math_sym3333to66(m3333,weighted) logical, optional, intent(in) :: weighted real(pReal), dimension(6) :: w - integer(pInt) :: i,j + integer :: i,j if(present(weighted)) then w = merge(nrmMandel,1.0_pReal,weighted) @@ -1180,7 +1178,7 @@ pure function math_sym3333to66(m3333,weighted) w = nrmMandel endif - forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) & + forall(i=1:6,j=1:6) & math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) end function math_sym3333to66 @@ -1200,7 +1198,7 @@ pure function math_66toSym3333(m66,weighted) logical, optional, intent(in) :: weighted real(pReal), dimension(6) :: w - integer(pInt) :: i,j + integer :: i,j if(present(weighted)) then w = merge(invnrmMandel,1.0_pReal,weighted) @@ -1208,7 +1206,7 @@ pure function math_66toSym3333(m66,weighted) w = invnrmMandel endif - do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt + do i=1,6; do j=1, 6 math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j) math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j) @@ -1226,9 +1224,9 @@ pure function math_Voigt66to3333(m66) implicit none real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 real(pReal), dimension(6,6), intent(in) :: m66 - integer(pInt) :: i,j + integer :: i,j - do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt + do i=1,6; do j=1, 6 math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j) math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j) @@ -1250,7 +1248,7 @@ function math_qRand() real(pReal), dimension(4) :: math_qRand real(pReal), dimension(3) :: rnd - rnd = halton([8_pInt,4_pInt,9_pInt]) + rnd = halton([8,4,9]) math_qRand = [cos(2.0_pReal*PI*rnd(1))*sqrt(rnd(3)), & sin(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), & cos(2.0_pReal*PI*rnd(2))*sqrt(1.0_pReal-rnd(3)), & @@ -1346,10 +1344,10 @@ pure function math_qRot(Q,v) real(pReal), dimension(3), intent(in) :: v real(pReal), dimension(3) :: math_qRot real(pReal), dimension(4,4) :: T - integer(pInt) :: i, j + integer :: i, j - do i = 1_pInt,4_pInt - do j = 1_pInt,i + do i = 1,4 + do j = 1,i T(i,j) = Q(i) * Q(j) enddo enddo @@ -1408,7 +1406,7 @@ pure function math_RtoQ(R) real(pReal), dimension(3,3), intent(in) :: R real(pReal), dimension(4) :: absQ, math_RtoQ real(pReal) :: max_absQ - integer, dimension(1) :: largest !no pInt, maxloc returns integer default + integer, dimension(1) :: largest math_RtoQ = 0.0_pReal @@ -1639,9 +1637,9 @@ pure function math_qToR(q) implicit none real(pReal), dimension(4), intent(in) :: q real(pReal), dimension(3,3) :: math_qToR, T,S - integer(pInt) :: i, j + integer :: i, j - forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) T(i,j) = q(i+1_pInt) * q(j+1_pInt) + forall(i = 1:3, j = 1:3) T(i,j) = q(i+1) * q(j+1) S = reshape( [0.0_pReal, -q(4), q(3), & q(4), 0.0_pReal, -q(2), & @@ -1772,7 +1770,7 @@ function math_sampleRandomOri() implicit none real(pReal), dimension(3) :: math_sampleRandomOri, rnd - rnd = halton([1_pInt,7_pInt,3_pInt]) + rnd = halton([1,7,3]) math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, & acos(2.0_pReal*rnd(2)-1.0_pReal), & rnd(3)*2.0_pReal*PI] @@ -1800,7 +1798,7 @@ function math_sampleGaussOri(center,FWHM) math_sampleGaussOri = center else GaussConvolution: do - rnd = halton([8_pInt,3_pInt,6_pInt,11_pInt]) + rnd = halton([8,3,6,11]) axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1] axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),& sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis @@ -1830,10 +1828,10 @@ function math_sampleFiberOri(alpha,beta,FWHM) u real(pReal), dimension(3) :: rnd real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt - integer(pInt), dimension(:),allocatable :: idx !< components of 2D vector + integer, dimension(:),allocatable :: idx !< components of 2D vector real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components) real(pReal):: angle,c - integer(pInt):: j,& !< index of smallest component + integer:: j,& !< index of smallest component i allocate(a(0)) @@ -1843,11 +1841,11 @@ function math_sampleFiberOri(alpha,beta,FWHM) R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system - rnd = halton([7_pInt,10_pInt,3_pInt]) + rnd = halton([7,10,3]) R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis if (FWHM > 0.1_pReal*INRAD) then - reducedTo2D: do i=1_pInt,3_pInt + reducedTo2D: do i=1,3 if (i /= minloc(abs(fInS),1)) then a=[a,fInS(i)] idx=[idx,i] @@ -1868,7 +1866,7 @@ function math_sampleFiberOri(alpha,beta,FWHM) R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component exit endif rejectionSampling - rnd = halton([7_pInt,10_pInt,3_pInt]) + rnd = halton([7,10,3]) enddo GaussConvolution endif math_sampleFiberOri = math_RtoEuler(R) @@ -1897,7 +1895,7 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width) myWidth = merge(width,3.0_pReal,present(width)) ! use +-3*sigma as default value for scatter if not given do - rnd = halton([6_pInt,2_pInt]) + rnd = halton([6,2]) scatter = myWidth * (2.0_pReal * rnd(1) - 1.0_pReal) if (rnd(2) <= exp(-0.5_pReal * scatter ** 2.0_pReal)) exit ! test if scattered value is drawn enddo @@ -1915,7 +1913,7 @@ end function math_sampleGaussVar pure function math_symmetricEulers(sym,Euler) implicit none - integer(pInt), intent(in) :: sym !< symmetry Class + integer, intent(in) :: sym !< symmetry Class real(pReal), dimension(3), intent(in) :: Euler real(pReal), dimension(3,3) :: math_symmetricEulers @@ -1926,9 +1924,9 @@ pure function math_symmetricEulers(sym,Euler) math_symmetricEulers = modulo(math_symmetricEulers,2.0_pReal*pi) select case (sym) - case (4_pInt) ! orthotropic: all done + case (4) ! orthotropic: all done - case (2_pInt) ! monoclinic: return only first + case (2) ! monoclinic: return only first math_symmetricEulers(1:3,2:3) = 0.0_pReal case default ! triclinic: return blank @@ -1949,14 +1947,14 @@ subroutine math_eigenValuesVectorsSym(m,values,vectors,error) real(pReal), dimension(size(m,1)), intent(out) :: values real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: vectors logical, intent(out) :: error - integer(pInt) :: info + integer :: info real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f external :: & dsyev vectors = m ! copy matrix to input (doubles as output) array call dsyev('V','U',size(m,1),vectors,size(m,1),values,work,(64+2)*size(m,1),info) - error = (info == 0_pInt) + error = (info == 0) end subroutine math_eigenValuesVectorsSym @@ -1982,11 +1980,11 @@ subroutine math_eigenValuesVectorsSym33(m,values,vectors) vectors(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), & m(1, 3) * m(1, 2) - m(2, 3) * m(1, 1), & - m(1, 2)**2_pInt] + m(1, 2)**2] T = maxval(abs(values)) - U = max(T, T**2_pInt) - threshold = sqrt(5.68e-14_pReal * U**2_pInt) + U = max(T, T**2) + threshold = sqrt(5.68e-14_pReal * U**2) ! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2 vectors(1:3,1) = [ vectors(1,2) + m(1, 3) * values(1), & @@ -2030,13 +2028,13 @@ function math_eigenvectorBasisSym(m) real(pReal), dimension(size(m,1),size(m,1)) :: vectors real(pReal), dimension(size(m,1),size(m,1)) :: math_eigenvectorBasisSym logical :: error - integer(pInt) :: i + integer :: i math_eigenvectorBasisSym = 0.0_pReal call math_eigenValuesVectorsSym(m,values,vectors,error) if(error) return - do i=1_pInt, size(m,1) + do i=1, size(m,1) math_eigenvectorBasisSym = math_eigenvectorBasisSym & + sqrt(values(i)) * math_outer(vectors(:,i),vectors(:,i)) enddo @@ -2193,7 +2191,7 @@ function math_rotationalPart33(m) inversionFailed: if (all(dEq0(Uinv))) then math_rotationalPart33 = math_I3 - call IO_warning(650_pInt) + call IO_warning(650) else inversionFailed math_rotationalPart33 = math_mul33x33(m,Uinv) endif inversionFailed @@ -2213,14 +2211,14 @@ function math_eigenvaluesSym(m) real(pReal), dimension(:,:), intent(in) :: m real(pReal), dimension(size(m,1)) :: math_eigenvaluesSym real(pReal), dimension(size(m,1),size(m,1)) :: vectors - integer(pInt) :: info + integer :: info real(pReal), dimension((64+2)*size(m,1)) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f external :: & dsyev vectors = m ! copy matrix to input (doubles as output) array call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info) - if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) + if (info /= 0) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN) end function math_eigenvaluesSym @@ -2294,19 +2292,19 @@ end function math_invariantsSym33 function halton(bases) implicit none - integer(pInt), intent(in), dimension(:):: & + integer, intent(in), dimension(:):: & bases !< bases (prime number ID) real(pReal), dimension(size(bases)) :: & halton - integer(pInt), save :: & - current = 1_pInt + integer, save :: & + current = 1 real(pReal), dimension(size(bases)) :: & base_inv - integer(pInt), dimension(size(bases)) :: & + integer, dimension(size(bases)) :: & base, & t - integer(pInt), dimension(0:1600), parameter :: & - prime = int([& + integer, dimension(0:1600), parameter :: & + prime = [& 1, & 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, & 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, & @@ -2482,9 +2480,9 @@ function halton(bases) 13121, 13127, 13147, 13151, 13159, 13163, 13171, 13177, 13183, 13187, & 13217, 13219, 13229, 13241, 13249, 13259, 13267, 13291, 13297, 13309, & 13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399, 13411, & - 13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499],pInt) + 13417, 13421, 13441, 13451, 13457, 13463, 13469, 13477, 13487, 13499] - current = current + 1_pInt + current = current + 1 base = prime(bases) base_inv = 1.0_pReal/real(base,pReal) @@ -2492,7 +2490,7 @@ function halton(bases) halton = 0.0_pReal t = current - do while (any( t /= 0_pInt) ) + do while (any( t /= 0) ) halton = halton + real(mod(t,base), pReal) * base_inv base_inv = base_inv / real(base, pReal) t = t / base @@ -2504,11 +2502,11 @@ end function halton !-------------------------------------------------------------------------------------------------- !> @brief factorial !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function math_factorial(n) +integer pure function math_factorial(n) implicit none - integer(pInt), intent(in) :: n - integer(pInt) :: i + integer, intent(in) :: n + integer :: i math_factorial = product([(i, i=1,n)]) @@ -2518,11 +2516,11 @@ end function math_factorial !-------------------------------------------------------------------------------------------------- !> @brief binomial coefficient !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function math_binomial(n,k) +integer pure function math_binomial(n,k) implicit none - integer(pInt), intent(in) :: n, k - integer(pInt) :: i, j + integer, intent(in) :: n, k + integer :: i, j j = min(k,n-k) math_binomial = product([(i, i=n, n-j+1, -1)])/math_factorial(j) @@ -2533,13 +2531,13 @@ end function math_binomial !-------------------------------------------------------------------------------------------------- !> @brief multinomial coefficient !-------------------------------------------------------------------------------------------------- -integer(pInt) pure function math_multinomial(alpha) +integer pure function math_multinomial(alpha) implicit none - integer(pInt), intent(in), dimension(:) :: alpha - integer(pInt) :: i + integer, intent(in), dimension(:) :: alpha + integer :: i - math_multinomial = 1_pInt + math_multinomial = 1 do i = 1, size(alpha) math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i)) enddo @@ -2616,11 +2614,11 @@ pure function math_rotate_forward3333(tensor,rot_tensor) real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333 real(pReal), dimension(3,3), intent(in) :: rot_tensor real(pReal), dimension(3,3,3,3), intent(in) :: tensor - integer(pInt) :: i,j,k,l,m,n,o,p + integer :: i,j,k,l,m,n,o,p math_rotate_forward3333 = 0.0_pReal - do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt;do k = 1_pInt,3_pInt;do l = 1_pInt,3_pInt - do m = 1_pInt,3_pInt;do n = 1_pInt,3_pInt;do o = 1_pInt,3_pInt;do p = 1_pInt,3_pInt + do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 + do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 math_rotate_forward3333(i,j,k,l) & = math_rotate_forward3333(i,j,k,l) & + rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p) diff --git a/src/mesh_FEM.f90 b/src/mesh_FEM.f90 index 29f0bc682..af17b2866 100644 --- a/src/mesh_FEM.f90 +++ b/src/mesh_FEM.f90 @@ -120,9 +120,7 @@ subroutine tMesh_FEM_init(self,dimen,order,nodes) !-------------------------------------------------------------------------------------------------- subroutine mesh_init() use DAMASK_interface - 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_timeStamp, & IO_error, & IO_open_file, & IO_stringPos, & @@ -161,8 +159,6 @@ subroutine mesh_init() write(6,'(/,a)') ' <<<+- mesh init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" ! read in file call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr) diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 0b3df43ef..f8a64b55b 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -19,19 +19,10 @@ contains !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- subroutine plastic_none_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use prec, only: & - pInt use debug, only: & debug_level, & debug_constitutive, & debug_levelBasic - use IO, only: & - IO_timeStamp use material, only: & phase_plasticity, & material_allocatePlasticState, & @@ -41,28 +32,26 @@ subroutine plastic_none_init plasticState implicit none - integer(pInt) :: & + integer :: & Ninstance, & p, & NipcMyPhase write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - Ninstance = int(count(phase_plasticity == PLASTICITY_NONE_ID),pInt) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - do p = 1_pInt, size(phase_plasticity) + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phase == p) - call material_allocatePlasticState(p,NipcMyPhase,0_pInt,0_pInt,0_pInt, & - 0_pInt,0_pInt,0_pInt) - plasticState(p)%sizePostResults = 0_pInt + call material_allocatePlasticState(p,NipcMyPhase,0,0,0, & + 0,0,0) + plasticState(p)%sizePostResults = 0 enddo diff --git a/src/quit.f90 b/src/quit.f90 index ad61943e4..0215737e5 100644 --- a/src/quit.f90 +++ b/src/quit.f90 @@ -11,22 +11,19 @@ subroutine quit(stop_id) use MPI, only: & MPI_finalize #endif - use prec, only: & - pInt use PetscSys use hdf5 implicit none - integer(pInt), intent(in) :: stop_id + integer, intent(in) :: stop_id integer, dimension(8) :: dateAndTime ! type default integer - integer :: hdferr - integer(pInt) :: error = 0_pInt + integer :: error PetscErrorCode :: ierr = 0 - call h5open_f(hdferr) - if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5open_f',hdferr ! prevents error if not opened yet - call h5close_f(hdferr) - if (hdferr /= 0) write(6,'(a,i5)') ' Error in h5close_f',hdferr + call h5open_f(error) + if (error /= 0) write(6,'(a,i5)') ' Error in h5open_f ',error ! prevents error if not opened yet + call h5close_f(error) + if (error /= 0) write(6,'(a,i5)') ' Error in h5close_f ',error call PETScFinalize(ierr) CHKERRQ(ierr) @@ -45,8 +42,8 @@ subroutine quit(stop_id) dateAndTime(6),':',& dateAndTime(7) - if (stop_id == 0_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 0 ! normal termination - if (stop_id == 2_pInt .and. ierr == 0_pInt .and. error == 0_pInt) stop 2 ! not all incs converged + if (stop_id == 0 .and. ierr == 0 .and. error == 0) stop 0 ! normal termination + if (stop_id == 2 .and. ierr == 0 .and. error == 0) stop 2 ! not all incs converged stop 1 ! error (message from IO_error) end subroutine quit diff --git a/src/rotations.f90 b/src/rotations.f90 index 55602b557..470d82efa 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -383,15 +383,13 @@ end function om2eu !> @brief convert axis angle pair to orientation matrix !--------------------------------------------------------------------------------------------------- pure function ax2om(ax) result(om) - use prec, only: & - pInt implicit none real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3,3) :: om real(pReal) :: q, c, s, omc - integer(pInt) :: i + integer :: i c = cos(ax(4)) s = sin(ax(4)) @@ -476,13 +474,12 @@ end function ax2ho !--------------------------------------------------------------------------------------------------- pure function ho2ax(ho) result(ax) use prec, only: & - pInt, & dEq0 implicit none real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(4) :: ax - integer(pInt) :: i + integer :: i real(pReal) :: hmag_squared, s, hm real(pReal), parameter, dimension(16) :: & tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & @@ -519,7 +516,6 @@ end function ho2ax !--------------------------------------------------------------------------------------------------- function om2ax(om) result(ax) use prec, only: & - pInt, & dEq0, & cEq, & dNeq0 @@ -537,7 +533,7 @@ function om2ax(om) result(ax) real(pReal), dimension(3) :: Wr, Wi real(pReal), dimension(10) :: WORK real(pReal), dimension(3,3) :: VR, devNull, o - integer(pInt) :: INFO, LWORK, i + integer :: INFO, LWORK, i external :: dgeev,sgeev @@ -557,7 +553,7 @@ function om2ax(om) result(ax) ! call the eigenvalue solver call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO) - if (INFO /= 0) call IO_error(0_pInt,ext_msg='Error in om2ax DGEEV return not zero') + if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero') i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc ax(1:3) = VR(1:3,i) where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & diff --git a/src/spectral_damage.f90 b/src/spectral_damage.f90 deleted file mode 100644 index 212d36938..000000000 --- a/src/spectral_damage.f90 +++ /dev/null @@ -1,358 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Spectral solver for nonlocal damage -!-------------------------------------------------------------------------------------------------- -module spectral_damage -#include -#include - use PETScdmda - use PETScsnes - use prec, only: & - pInt, & - pReal - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private - - character (len=*), parameter, public :: & - spectral_damage_label = 'spectraldamage' - -!-------------------------------------------------------------------------------------------------- -! derived types - type(tSolutionParams), private :: params - -!-------------------------------------------------------------------------------------------------- -! PETSc data - SNES, private :: damage_snes - Vec, private :: solution - PetscInt, private :: xstart, xend, ystart, yend, zstart, zend - real(pReal), private, dimension(:,:,:), allocatable :: & - damage_current, & !< field of current damage - damage_lastInc, & !< field of previous damage - damage_stagInc !< field of staggered damage - -!-------------------------------------------------------------------------------------------------- -! reference diffusion tensor, mobility etc. - integer(pInt), private :: totalIter = 0 !< total iteration in current increment - real(pReal), dimension(3,3), private :: D_ref - real(pReal), private :: mobility_ref - - public :: & - spectral_damage_init, & - spectral_damage_solution, & - spectral_damage_forward - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields and fills them with data, potentially from restart info -!-------------------------------------------------------------------------------------------------- -subroutine spectral_damage_init() - use IO, only: & - IO_intOut - use spectral_utilities, only: & - wgt - use mesh, only: & - grid, & - grid3 - use damage_nonlocal, only: & - damage_nonlocal_getDiffusion33, & - damage_nonlocal_getMobility - use numerics, only: & - worldrank, & - worldsize - - implicit none - PetscInt, dimension(worldsize) :: localK - integer :: i, j, k, cell - DM :: damage_grid - Vec :: uBound, lBound - PetscErrorCode :: ierr - character(len=100) :: snes_type - - write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' - - write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' - write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' - -!-------------------------------------------------------------------------------------------------- -! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 - call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) - call DMDACreate3D(PETSC_COMM_WORLD, & - DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary - DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point - grid(1),grid(2),grid(3), & !< global grid - 1, 1, worldsize, & - 1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & !< local grid - damage_grid,ierr) !< handle, error - CHKERRQ(ierr) - call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da - call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr) - call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) - call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor) - call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,& - PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector - CHKERRQ(ierr) - call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments - call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) - if (trim(snes_type) == 'vinewtonrsls' .or. & - trim(snes_type) == 'vinewtonssls') then - call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) - call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) - call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) - call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) - call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. - call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) - call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) - endif - -!-------------------------------------------------------------------------------------------------- -! init fields - call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) - CHKERRQ(ierr) - xend = xstart + xend - 1 - yend = ystart + yend - 1 - zend = zstart + zend - 1 - call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr) - allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) - allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) - allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) - -!-------------------------------------------------------------------------------------------------- -! damage reference diffusion update - cell = 0_pInt - D_ref = 0.0_pReal - mobility_ref = 0.0_pReal - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) - mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) - enddo; enddo; enddo - D_ref = D_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - mobility_ref = mobility_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - -end subroutine spectral_damage_init - -!-------------------------------------------------------------------------------------------------- -!> @brief solution for the spectral damage scheme with internal iterations -!-------------------------------------------------------------------------------------------------- -type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadCaseTime) - use numerics, only: & - itmax, & - err_damage_tolAbs, & - err_damage_tolRel - use mesh, only: & - grid, & - grid3 - use damage_nonlocal, only: & - damage_nonlocal_putNonLocalDamage - - implicit none - 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 - - integer :: i, j, k, cell - PetscInt ::position - PetscReal :: minDamage, maxDamage, stagNorm, solnNorm - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - spectral_damage_solution%converged =.false. - -!-------------------------------------------------------------------------------------------------- -! set module wide availabe data - params%timeinc = timeinc - params%timeincOld = timeinc_old - - call SNESSolve(damage_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) - call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) - - if (reason < 1) then - spectral_damage_solution%converged = .false. - spectral_damage_solution%iterationsNeeded = itmax - else - spectral_damage_solution%converged = .true. - spectral_damage_solution%iterationsNeeded = totalIter - endif - stagNorm = maxval(abs(damage_current - damage_stagInc)) - solnNorm = maxval(abs(damage_current)) - call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) - call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) - damage_stagInc = damage_current - spectral_damage_solution%stagConverged = stagNorm < err_damage_tolAbs & - .or. stagNorm < err_damage_tolRel*solnNorm - -!-------------------------------------------------------------------------------------------------- -! updating damage state - cell = 0 - do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 - call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) - enddo; enddo; enddo - - call VecMin(solution,position,minDamage,ierr); CHKERRQ(ierr) - call VecMax(solution,position,maxDamage,ierr); CHKERRQ(ierr) - if (spectral_damage_solution%converged) & - write(6,'(/,a)') ' ... nonlocal damage converged .....................................' - write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& - minDamage, maxDamage, stagNorm - write(6,'(/,a)') ' ===========================================================================' - flush(6) - -end function spectral_damage_solution - - -!-------------------------------------------------------------------------------------------------- -!> @brief forms the spectral damage residual vector -!-------------------------------------------------------------------------------------------------- -subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr) - use numerics, only: & - residualStiffness - use mesh, only: & - grid, & - grid3 - use math, only: & - math_mul33x3 - use spectral_utilities, only: & - scalarField_real, & - vectorField_real, & - utilities_FFTvectorForward, & - utilities_FFTvectorBackward, & - utilities_FFTscalarForward, & - utilities_FFTscalarBackward, & - utilities_fourierGreenConvolution, & - utilities_fourierScalarGradient, & - utilities_fourierVectorDivergence - use damage_nonlocal, only: & - damage_nonlocal_getSourceAndItsTangent,& - damage_nonlocal_getDiffusion33, & - damage_nonlocal_getMobility - - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & - in - PetscScalar, dimension( & - XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: & - x_scal - PetscScalar, dimension( & - X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & - f_scal - PetscObject :: dummy - PetscErrorCode :: ierr - integer(pInt) :: i, j, k, cell - real(pReal) :: phiDot, dPhiDot_dPhi, mobility - - damage_current = x_scal -!-------------------------------------------------------------------------------------------------- -! evaluate polarization field - scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current - call utilities_FFTscalarForward() - call utilities_fourierScalarGradient() !< calculate gradient of damage field - call utilities_FFTvectorBackward() - cell = 0_pInt - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt - vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & - vectorField_real(1:3,i,j,k)) - enddo; enddo; enddo - call utilities_FFTvectorForward() - call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field - call utilities_FFTscalarBackward() - cell = 0_pInt - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt - call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell) - mobility = damage_nonlocal_getMobility(1,cell) - scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & - params%timeinc*phiDot + & - mobility*damage_lastInc(i,j,k) - & - mobility*damage_current(i,j,k) + & - mobility_ref*damage_current(i,j,k) - enddo; enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! convolution of damage field with green operator - call utilities_FFTscalarForward() - call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) - call utilities_FFTscalarBackward() - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc - where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & - scalarField_real(1:grid(1),1:grid(2),1:grid3) = residualStiffness - -!-------------------------------------------------------------------------------------------------- -! constructing residual - f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current - -end subroutine spectral_damage_formResidual - -!-------------------------------------------------------------------------------------------------- -!> @brief spectral damage forwarding routine -!-------------------------------------------------------------------------------------------------- -subroutine spectral_damage_forward() - use mesh, only: & - grid, & - grid3 - use spectral_utilities, only: & - cutBack, & - wgt - use damage_nonlocal, only: & - damage_nonlocal_putNonLocalDamage, & - damage_nonlocal_getDiffusion33, & - damage_nonlocal_getMobility - - implicit none - integer(pInt) :: i, j, k, cell - DM :: dm_local - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr - - if (cutBack) then - damage_current = damage_lastInc - damage_stagInc = damage_lastInc -!-------------------------------------------------------------------------------------------------- -! reverting damage field state - cell = 0_pInt - call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with - x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current - call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt - call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) - enddo; enddo; enddo - else -!-------------------------------------------------------------------------------------------------- -! update rate and forward last inc - damage_lastInc = damage_current - cell = 0_pInt - D_ref = 0.0_pReal - mobility_ref = 0.0_pReal - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt - D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) - mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) - enddo; enddo; enddo - D_ref = D_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - mobility_ref = mobility_ref*wgt - call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - endif - -end subroutine spectral_damage_forward - -end module spectral_damage diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 deleted file mode 100644 index 4f32b098d..000000000 --- a/src/spectral_mech_Basic.f90 +++ /dev/null @@ -1,541 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @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 -!-------------------------------------------------------------------------------------------------- -module spectral_mech_basic -#include -#include - use PETScdmda - use PETScsnes - use prec, only: & - pInt, & - pReal - use math, only: & - math_I3 - use spectral_utilities, only: & - tSolutionState, & - tSolutionParams - - implicit none - private - - character (len=*), parameter, public :: & - DAMASK_spectral_SolverBasic_label = 'basic' - -!-------------------------------------------------------------------------------------------------- -! derived types - type(tSolutionParams), private :: params - -!-------------------------------------------------------------------------------------------------- -! PETSc data - DM, private :: da - SNES, private :: snes - Vec, private :: solution_vec - -!-------------------------------------------------------------------------------------------------- -! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot - -!-------------------------------------------------------------------------------------------------- -! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient - F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - - character(len=1024), private :: incInfo !< time and increment information - - real(pReal), private, dimension(3,3,3,3) :: & - C_volAvg = 0.0_pReal, & !< current volume average stiffness - C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness - C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness - C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness - S = 0.0_pReal !< current compliance (filled up with zeros) - - real(pReal), private :: & - err_BC, & !< deviation from stress BC - err_div !< RMS of div of P - - integer(pInt), private :: & - totalIter = 0_pInt !< total iteration in current increment - - public :: & - basic_init, & - basic_solution, & - basic_forward - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all necessary fields and fills them with data, potentially from restart info -!-------------------------------------------------------------------------------------------------- -subroutine basic_init - use IO, only: & - IO_intOut, & - IO_error, & - IO_open_jobFile_binary - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRestart - use FEsolving, only: & - restartInc - use numerics, only: & - worldrank, & - worldsize - use homogenization, only: & - materialpoint_F0 - use DAMASK_interface, only: & - getSolverJobName - use spectral_utilities, only: & - Utilities_constitutiveResponse, & - Utilities_updateGamma, & - Utilities_updateIPcoords, & - wgt - use mesh, only: & - grid, & - grid3 - use math, only: & - math_invSym3333 - - implicit none - real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - - PetscErrorCode :: ierr - PetscScalar, pointer, dimension(:,:,:,:) :: F - PetscInt, dimension(worldsize) :: localK - integer :: fileUnit - character(len=1024) :: rankStr - - write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' - - write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012' - - write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' - write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' - -!-------------------------------------------------------------------------------------------------- -! allocate global fields - allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) - -!-------------------------------------------------------------------------------------------------- -! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 - call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) - call DMDACreate3d(PETSC_COMM_WORLD, & - DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary - DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point - grid(1),grid(2),grid(3), & ! global grid - 1 , 1, worldsize, & - 9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - [grid(1)],[grid(2)],localK, & ! local grid - da,ierr) ! handle, error - CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da - call DMsetFromOptions(da,ierr); CHKERRQ(ierr) - call DMsetUp(da,ierr); CHKERRQ(ierr) - call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) - call DMDASNESsetFunctionLocal(da,INSERT_VALUES,Basic_formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,Basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" - CHKERRQ(ierr) - call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments - -!-------------------------------------------------------------------------------------------------- -! init fields - call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with - - restart: if (restartInc > 0_pInt) then - if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & - 'reading values of increment ', restartInc, ' from file' - flush(6) - endif - - fileUnit = IO_open_jobFile_binary('F_aimDot') - read(fileUnit) F_aimDot; close(fileUnit) - - write(rankStr,'(a1,i0)')'_',worldrank - - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr)) - read(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr)) - read(fileUnit) F_lastInc; close (fileUnit) - - F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F - call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim') - F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc - call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) - if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc') - elseif (restartInc == 0_pInt) then restart - F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity - F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) - endif restart - - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 - reshape(F,shape(F_lastInc)), & ! target F - 0.0_pReal, & ! time increment - math_I3) ! no rotation of boundary condition - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc - ! QUESTION: why not writing back right after reading (l.189)? - - restartRead: if (restartInc > 0_pInt) then ! QUESTION: are those values not calc'ed by constitutiveResponse? why reading from file? - if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 .and. worldrank == 0_pInt) & - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & - 'reading more values of increment ', restartInc, ' from file' - flush(6) - fileUnit = IO_open_jobFile_binary('C_volAvg') - read(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') - read(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_ref') - read(fileUnit) C_minMaxAvg; close(fileUnit) - endif restartRead - - call Utilities_updateGamma(C_minMaxAvg,.true.) - -end subroutine basic_init - -!-------------------------------------------------------------------------------------------------- -!> @brief solution for the Basic scheme with internal iterations -!-------------------------------------------------------------------------------------------------- -type(tSolutionState) function basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) - use IO, only: & - IO_error - use numerics, only: & - update_gamma - use spectral_utilities, only: & - tBoundaryCondition, & - Utilities_maskedCompliance, & - Utilities_updateGamma - use FEsolving, only: & - restartWrite, & - terminallyIll - - implicit none - -!-------------------------------------------------------------------------------------------------- -! input data for solution - character(len=*), intent(in) :: & - incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< increment time for current solution - timeinc_old !< increment time of last successful increment - type(tBoundaryCondition), intent(in) :: & - stress_BC - real(pReal), dimension(3,3), intent(in) :: rotation_BC - -!-------------------------------------------------------------------------------------------------- -! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - incInfo = incInfoIn - -!-------------------------------------------------------------------------------------------------- -! update stiffness (and gamma operator) - S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) - if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite) - - -!-------------------------------------------------------------------------------------------------- -! set module wide availabe data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - -!-------------------------------------------------------------------------------------------------- -! solve BVP - call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) - -!-------------------------------------------------------------------------------------------------- -! check convergence - call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) - - basic_solution%converged = reason > 0 - basic_solution%iterationsNeeded = totalIter - basic_solution%termIll = terminallyIll - terminallyIll = .false. - if (reason == -4) call IO_error(893_pInt) ! MPI error - -end function basic_solution - - -!-------------------------------------------------------------------------------------------------- -!> @brief forms the basic residual vector -!-------------------------------------------------------------------------------------------------- -subroutine Basic_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) - F, & ! defgrad field on grid - residuum, & ! residuum field on grid - dummy, & - ierr) - use numerics, only: & - itmax, & - itmin - use mesh, only: & - grid, & - grid3 - use math, only: & - math_rotate_backward33, & - math_mul3333xx33 - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRotation - use spectral_utilities, only: & - tensorField_real, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & - utilities_FFTtensorBackward, & - Utilities_constitutiveResponse, & - Utilities_divergenceRMS - use IO, only: & - IO_intOut - use FEsolving, only: & - terminallyIll - - implicit none - DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in - PetscScalar, & - dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: F - PetscScalar, & - dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: residuum - PetscInt :: & - PETScIter, & - nfuncs - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal), dimension(3,3) :: & - deltaF_aim - - call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) - - if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment -!-------------------------------------------------------------------------------------------------- -! begin of new iteration - newIteration: if (totalIter <= PETScIter) then - totalIter = totalIter + 1_pInt - write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & - trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax - if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) - write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & - ' deformation gradient aim =', transpose(F_aim) - flush(6) - endif newIteration - -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F,params%timeinc,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) - -!-------------------------------------------------------------------------------------------------- -! stress BC handling - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc - -!-------------------------------------------------------------------------------------------------- -! updated deformation gradient using fix point algorithm of basic scheme - tensorField_real = 0.0_pReal - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform - call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real" - err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use - call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg - call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier - -!-------------------------------------------------------------------------------------------------- -! constructing residual - residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too - -end subroutine Basic_formResidual - - -!-------------------------------------------------------------------------------------------------- -!> @brief convergence check -!-------------------------------------------------------------------------------------------------- -subroutine Basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use FEsolving, only: & - terminallyIll - - implicit none - SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - divTol, & - BCTol - - divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs) - - converged: if ((totalIter >= itmin .and. & - all([ err_div/divTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then converged - reason = -1 - else converged - reason = 0 - endif converged - -!-------------------------------------------------------------------------------------------------- -! 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_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' -flush(6) - -end subroutine Basic_converged - -!-------------------------------------------------------------------------------------------------- -!> @brief forwarding routine -!> @details find new boundary conditions and best F estimate for end of current timestep -!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates -!-------------------------------------------------------------------------------------------------- -subroutine Basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) - use math, only: & - math_mul33x33 ,& - math_rotate_backward33 - use numerics, only: & - worldrank - use homogenization, only: & - materialpoint_F0 - use mesh, only: & - grid, & - grid3 - use CPFEM2, only: & - CPFEM_age - use spectral_utilities, only: & - Utilities_calculateRate, & - Utilities_forwardField, & - Utilities_updateIPcoords, & - tBoundaryCondition, & - cutBack - use IO, only: & - IO_open_jobFile_binary - use FEsolving, only: & - restartWrite - - implicit none - logical, intent(in) :: & - guess - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - type(tBoundaryCondition), intent(in) :: & - stress_BC, & - deformation_BC - real(pReal), dimension(3,3), intent(in) :: & - rotation_BC - PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: F - - integer :: fileUnit - character(len=32) :: rankStr - - call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - - if (cutBack) then - C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? - C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? - else - !-------------------------------------------------------------------------------------------------- - ! restart information for spectral solver - if (restartWrite) then ! QUESTION: where is this logical properly set? - write(6,'(/,a)') ' writing converged results for restart' - flush(6) - - if (worldrank == 0) then - fileUnit = IO_open_jobFile_binary('C_volAvg','w') - write(fileUnit) C_volAvg; close(fileUnit) - fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w') - write(fileUnit) C_volAvgLastInc; close(fileUnit) - fileUnit = IO_open_jobFile_binary('F_aimDot','w') - write(fileUnit) F_aimDot; close(fileUnit) - endif - - write(rankStr,'(a1,i0)')'_',worldrank - fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w') - write(fileUnit) F; close (fileUnit) - fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w') - write(fileUnit) F_lastInc; close (fileUnit) - endif - - call CPFEM_age() ! age state and kinematics - call utilities_updateIPcoords(F) - - C_volAvgLastInc = C_volAvg - C_minMaxAvgLastInc = C_minMaxAvg - - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) - F_aim_lastInc = F_aim - - !-------------------------------------------------------------------------------------------------- - ! calculate rate for aim - if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) - elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values - elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime - endif - - - Fdot = Utilities_calculateRate(guess, & - F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, & - math_rotate_backward33(F_aimDot,rotation_BC)) - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward - materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - endif - -!-------------------------------------------------------------------------------------------------- -! update average and local deformation gradients - F_aim = F_aim_lastInc + F_aimDot * timeinc - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average - math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3]) - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - -end subroutine Basic_forward - -end module spectral_mech_basic