diff --git a/src/spectral_damage.f90 b/src/spectral_damage.f90 index 1b741ee2d..4c581350e 100644 --- a/src/spectral_damage.f90 +++ b/src/spectral_damage.f90 @@ -11,14 +11,9 @@ module spectral_damage use prec, only: & pInt, & pReal - use math, only: & - math_I3 use spectral_utilities, only: & tSolutionState, & tSolutionParams - use numerics, only: & - worldrank, & - worldsize implicit none private @@ -42,7 +37,7 @@ module spectral_damage !-------------------------------------------------------------------------------------------------- ! reference diffusion tensor, mobility etc. - integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment + integer(pInt), private :: totalIter = 0 !< total iteration in current increment real(pReal), dimension(3,3), private :: D_ref real(pReal), private :: mobility_ref @@ -57,93 +52,94 @@ 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 - - implicit none - PetscInt, dimension(:), allocatable :: localK - integer(pInt) :: proc - integer(pInt) :: 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, volume in press, ' - write(6,'(a,/)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 ' + 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, volume in press, ' + write(6,'(a,/)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 ' !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 - do proc = 1, worldsize - call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) - enddo - 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 + 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) + 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_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) + 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 @@ -151,74 +147,69 @@ 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 + 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 - implicit none - -!-------------------------------------------------------------------------------------------------- -! input data for solution - real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old, & !< increment in time of last increment - loadCaseTime !< remaining time of current load case - integer(pInt) :: i, j, k, cell - PetscInt ::position - PetscReal :: minDamage, maxDamage, stagNorm, solnNorm - -!-------------------------------------------------------------------------------------------------- -! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason + 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 + 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_pInt !< material point = 0 - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt !< material point increase - 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) + 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 diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90 index 9508b6f0a..ee7116486 100644 --- a/src/spectral_mech_Basic.f90 +++ b/src/spectral_mech_Basic.f90 @@ -108,8 +108,8 @@ subroutine basic_init PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: F - PetscInt, dimension(:), allocatable :: localK - integer :: proc, fileUnit + PetscInt, dimension(worldsize) :: localK + integer :: fileUnit character(len=1024) :: rankStr write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' @@ -125,10 +125,9 @@ subroutine basic_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 - do proc = 1, worldsize !ToDo: there are smarter options in MPI - call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) - enddo + 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 diff --git a/src/spectral_mech_Polarisation.f90 b/src/spectral_mech_Polarisation.f90 index 3b613263f..c36f2b716 100644 --- a/src/spectral_mech_Polarisation.f90 +++ b/src/spectral_mech_Polarisation.f90 @@ -118,8 +118,8 @@ subroutine Polarisation_init FandF_tau, & ! overall pointer to solution data F, & ! specific (sub)pointer F_tau ! specific (sub)pointer - PetscInt, dimension(:), allocatable :: localK - integer :: proc, fileUnit + PetscInt, dimension(worldsize) :: localK + integer :: fileUnit character(len=1024) :: rankStr write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' @@ -137,10 +137,9 @@ subroutine Polarisation_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 - do proc = 1, worldsize !ToDo: there are smarter options in MPI - call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) - enddo + 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 diff --git a/src/spectral_thermal.f90 b/src/spectral_thermal.f90 index 2ad5c4eab..cb73d73ec 100644 --- a/src/spectral_thermal.f90 +++ b/src/spectral_thermal.f90 @@ -9,16 +9,10 @@ module spectral_thermal use PETScdmda use PETScsnes use prec, only: & - pInt, & pReal - use math, only: & - math_I3 use spectral_utilities, only: & tSolutionState, & tSolutionParams - use numerics, only: & - worldrank, & - worldsize implicit none private @@ -42,7 +36,7 @@ module spectral_thermal !-------------------------------------------------------------------------------------------------- ! reference diffusion tensor, mobility etc. - integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment + integer, private :: totalIter = 0 !< total iteration in current increment real(pReal), dimension(3,3), private :: D_ref real(pReal), private :: mobility_ref @@ -57,104 +51,96 @@ contains !> @brief allocates all neccessary fields and fills them with data, potentially from restart info !-------------------------------------------------------------------------------------------------- subroutine spectral_thermal_init -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use IO, only: & - IO_timeStamp - use spectral_utilities, only: & - wgt - use mesh, only: & - grid, & - grid3 - use thermal_conduction, only: & - thermal_conduction_getConductivity33, & - thermal_conduction_getMassDensity, & - thermal_conduction_getSpecificHeat - use material, only: & - mappingHomogenization, & - temperature, & - thermalMapping - - implicit none - integer(pInt), dimension(:), allocatable :: localK - integer(pInt) :: proc - integer(pInt) :: i, j, k, cell - DM :: thermal_grid - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr + use spectral_utilities, only: & + wgt + use mesh, only: & + grid, & + grid3 + use thermal_conduction, only: & + thermal_conduction_getConductivity33, & + thermal_conduction_getMassDensity, & + thermal_conduction_getSpecificHeat + use material, only: & + mappingHomogenization, & + temperature, & + thermalMapping + use numerics, only: & + worldrank, & + worldsize - write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' - write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,' - write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" + implicit none + integer, dimension(worldsize) :: localK + integer :: i, j, k, cell + DM :: thermal_grid + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr + + write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' + write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,' + write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018' !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) - allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 - do proc = 1, worldsize - call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) - enddo - 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 (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 - CHKERRQ(ierr) - call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments + call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) + call SNESSetOptionsPrefix(thermal_snes,'thermal_',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 (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 + CHKERRQ(ierr) + call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- ! init fields - call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr) - CHKERRQ(ierr) - xend = xstart + xend - 1 - yend = ystart + yend - 1 - zend = zstart + zend - 1 - allocate(temperature_current(grid(1),grid(2),grid3), source=0.0_pReal) - 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_pInt - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt - temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% & - p(thermalMapping(mappingHomogenization(2,1,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 - x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current - call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) + call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr) + CHKERRQ(ierr) + xend = xstart + xend - 1 + yend = ystart + yend - 1 + zend = zstart + zend - 1 + allocate(temperature_current(grid(1),grid(2),grid3), source=0.0_pReal) + 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) + cell = cell + 1 + temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% & + p(thermalMapping(mappingHomogenization(2,1,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 + x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current + call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! thermal reference diffusion update - 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 + 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) + 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) end subroutine spectral_thermal_init @@ -162,76 +148,72 @@ end subroutine spectral_thermal_init !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,loadCaseTime) - use numerics, only: & - itmax, & - err_thermal_tolAbs, & - err_thermal_tolRel - use mesh, only: & - grid, & - grid3 - use thermal_conduction, only: & - thermal_conduction_putTemperatureAndItsRate + use numerics, only: & + itmax, & + err_thermal_tolAbs, & + err_thermal_tolRel + use mesh, only: & + grid, & + grid3 + use thermal_conduction, only: & + thermal_conduction_putTemperatureAndItsRate + + implicit none - 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 :: minTemperature, maxTemperature, stagNorm, solnNorm -!-------------------------------------------------------------------------------------------------- -! input data for solution - real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old, & !< increment in time of last increment - loadCaseTime !< remaining time of current load case - integer(pInt) :: i, j, k, cell - PetscInt :: position - PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm + PetscErrorCode :: ierr + SNESConvergedReason :: reason -!-------------------------------------------------------------------------------------------------- -! PETSc Data - PetscErrorCode :: ierr - SNESConvergedReason :: reason - - spectral_thermal_solution%converged =.false. + spectral_thermal_solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data - params%timeinc = timeinc - params%timeincOld = timeinc_old + params%timeinc = timeinc + params%timeincOld = timeinc_old - call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) - call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) + call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,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 - else - spectral_thermal_solution%converged = .true. - spectral_thermal_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 + if (reason < 1) then + spectral_thermal_solution%converged = .false. + spectral_thermal_solution%iterationsNeeded = itmax + else + spectral_thermal_solution%converged = .true. + spectral_thermal_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 !-------------------------------------------------------------------------------------------------- ! updating thermal state - cell = 0_pInt !< material point = 0 - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt !< material point increase - 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 + cell = 0 + 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 - call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr) - call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr) - if (spectral_thermal_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) + call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr) + call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr) + if (spectral_thermal_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 @@ -272,7 +254,7 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) f_scal PetscObject :: dummy PetscErrorCode :: ierr - integer(pInt) :: i, j, k, cell + integer :: i, j, k, cell real(pReal) :: Tdot, dTdot_dT temperature_current = x_scal @@ -283,18 +265,18 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr) 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 + 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_pInt - do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt + 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 + & @@ -333,10 +315,10 @@ subroutine spectral_thermal_forward() thermal_conduction_getSpecificHeat implicit none - integer(pInt) :: i, j, k, cell - DM :: dm_local - PetscScalar, dimension(:,:,:), pointer :: x_scal - PetscErrorCode :: ierr + integer :: i, j, k, cell + DM :: dm_local + PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr if (cutBack) then temperature_current = temperature_lastInc @@ -344,13 +326,13 @@ subroutine spectral_thermal_forward() !-------------------------------------------------------------------------------------------------- ! reverting thermal field state - cell = 0_pInt !< material point = 0 + 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_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) - cell = cell + 1_pInt !< material point increase + 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, & @@ -360,11 +342,11 @@ subroutine spectral_thermal_forward() !-------------------------------------------------------------------------------------------------- ! update rate and forward last inc temperature_lastInc = temperature_current - cell = 0_pInt + cell = 0 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 + 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)