223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Spectral solver now fully parallel (parallel IO, domain decomposition, FFTs and restart). Working but not extensively tested so please report bugs to me

This commit is contained in:
Pratheek Shanthraj 2015-03-25 16:06:19 +00:00
parent 37a7364a3e
commit d44fce4a76
9 changed files with 1238 additions and 1038 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
code/CPFEM.f90
View File

@ -159,6 +159,7 @@ subroutine CPFEM_init
implicit none
integer(pInt) :: k,l,m,ph,homog
character(len=1024) :: rankStr
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
@ -181,39 +182,41 @@ subroutine CPFEM_init
flush(6)
endif
call IO_read_intFile(777,'recordedPhase',modelName,size(material_phase))
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_intFile(777,'recordedPhase'//trim(rankStr),modelName,size(material_phase))
read (777,rec=1) material_phase
close (777)
call IO_read_realFile(777,'convergedF',modelName,size(crystallite_F0))
call IO_read_realFile(777,'convergedF'//trim(rankStr),modelName,size(crystallite_F0))
read (777,rec=1) crystallite_F0
close (777)
call IO_read_realFile(777,'convergedFp',modelName,size(crystallite_Fp0))
call IO_read_realFile(777,'convergedFp'//trim(rankStr),modelName,size(crystallite_Fp0))
read (777,rec=1) crystallite_Fp0
close (777)
call IO_read_realFile(777,'convergedFi',modelName,size(crystallite_Fi0))
call IO_read_realFile(777,'convergedFi'//trim(rankStr),modelName,size(crystallite_Fi0))
read (777,rec=1) crystallite_Fi0
close (777)
call IO_read_realFile(777,'convergedLp',modelName,size(crystallite_Lp0))
call IO_read_realFile(777,'convergedLp'//trim(rankStr),modelName,size(crystallite_Lp0))
read (777,rec=1) crystallite_Lp0
close (777)
call IO_read_realFile(777,'convergedLi',modelName,size(crystallite_Li0))
call IO_read_realFile(777,'convergedLi'//trim(rankStr),modelName,size(crystallite_Li0))
read (777,rec=1) crystallite_Li0
close (777)
call IO_read_realFile(777,'convergeddPdF',modelName,size(crystallite_dPdF0))
call IO_read_realFile(777,'convergeddPdF'//trim(rankStr),modelName,size(crystallite_dPdF0))
read (777,rec=1) crystallite_dPdF0
close (777)
call IO_read_realFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v))
call IO_read_realFile(777,'convergedTstar'//trim(rankStr),modelName,size(crystallite_Tstar0_v))
read (777,rec=1) crystallite_Tstar0_v
close (777)
call IO_read_realFile(777,'convergedStateConst',modelName)
call IO_read_realFile(777,'convergedStateConst'//trim(rankStr),modelName)
m = 0_pInt
readPlasticityInstances: do ph = 1_pInt, size(phase_plasticity)
do k = 1_pInt, plasticState(ph)%sizeState
@ -224,7 +227,7 @@ subroutine CPFEM_init
enddo readPlasticityInstances
close (777)
call IO_read_realFile(777,'convergedStateHomog',modelName)
call IO_read_realFile(777,'convergedStateHomog'//trim(rankStr),modelName)
m = 0_pInt
readHomogInstances: do homog = 1_pInt, material_Nhomogenization
do k = 1_pInt, homogState(homog)%sizeState
@ -265,7 +268,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
#endif
use numerics, only: &
defgradTolerance, &
iJacoStiffness
iJacoStiffness, &
worldrank
use debug, only: &
debug_level, &
debug_CPFEM, &
@ -376,6 +380,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog
logical updateJaco ! flag indicating if JAcobian has to be updated
character(len=1024) :: rankStr
#if defined(Marc4DAMASK) || defined(Abaqus)
elCP = mesh_FEasCP('elem',elFE)
@ -443,39 +448,41 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) &
write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files'
call IO_write_jobRealFile(777,'recordedPhase',size(material_phase))
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'recordedPhase'//trim(rankStr),size(material_phase))
write (777,rec=1) material_phase
close (777)
call IO_write_jobRealFile(777,'convergedF',size(crystallite_F0))
call IO_write_jobRealFile(777,'convergedF'//trim(rankStr),size(crystallite_F0))
write (777,rec=1) crystallite_F0
close (777)
call IO_write_jobRealFile(777,'convergedFp',size(crystallite_Fp0))
call IO_write_jobRealFile(777,'convergedFp'//trim(rankStr),size(crystallite_Fp0))
write (777,rec=1) crystallite_Fp0
close (777)
call IO_write_jobRealFile(777,'convergedFi',size(crystallite_Fi0))
call IO_write_jobRealFile(777,'convergedFi'//trim(rankStr),size(crystallite_Fi0))
write (777,rec=1) crystallite_Fi0
close (777)
call IO_write_jobRealFile(777,'convergedLp',size(crystallite_Lp0))
call IO_write_jobRealFile(777,'convergedLp'//trim(rankStr),size(crystallite_Lp0))
write (777,rec=1) crystallite_Lp0
close (777)
call IO_write_jobRealFile(777,'convergedLi',size(crystallite_Li0))
call IO_write_jobRealFile(777,'convergedLi'//trim(rankStr),size(crystallite_Li0))
write (777,rec=1) crystallite_Li0
close (777)
call IO_write_jobRealFile(777,'convergeddPdF',size(crystallite_dPdF0))
call IO_write_jobRealFile(777,'convergeddPdF'//trim(rankStr),size(crystallite_dPdF0))
write (777,rec=1) crystallite_dPdF0
close (777)
call IO_write_jobRealFile(777,'convergedTstar',size(crystallite_Tstar0_v))
call IO_write_jobRealFile(777,'convergedTstar'//trim(rankStr),size(crystallite_Tstar0_v))
write (777,rec=1) crystallite_Tstar0_v
close (777)
call IO_write_jobRealFile(777,'convergedStateConst')
call IO_write_jobRealFile(777,'convergedStateConst'//trim(rankStr))
m = 0_pInt
writePlasticityInstances: do ph = 1_pInt, size(phase_plasticity)
do k = 1_pInt, plasticState(ph)%sizeState
@ -486,7 +493,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip)
enddo writePlasticityInstances
close (777)
call IO_write_jobRealFile(777,'convergedStateHomog')
call IO_write_jobRealFile(777,'convergedStateHomog'//trim(rankStr))
m = 0_pInt
writeHomogInstances: do homog = 1_pInt, material_Nhomogenization
do k = 1_pInt, homogState(homog)%sizeState

129
code/DAMASK_spectral_driver.f90
View File

@ -41,12 +41,17 @@ program DAMASK_spectral_Driver
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW
use mesh, only: &
gridGlobal, &
geomSizeGlobal
use CPFEM, only: &
CPFEM_initAll
use FEsolving, only: &
restartWrite, &
restartInc
use numerics, only: &
worldrank, &
worldsize, &
maxCutBack, &
spectral_solver, &
continueCalculation
@ -54,17 +59,17 @@ program DAMASK_spectral_Driver
materialpoint_sizeResults, &
materialpoint_results
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
tBoundaryCondition, &
tSolutionState, &
cutBack
use DAMASK_spectral_SolverBasic
use DAMASK_spectral_SolverBasicPETSC
use DAMASK_spectral_SolverAL
use DAMASK_spectral_SolverPolarisation
implicit none
#include <petsc-finclude/petscsys.h>
type tLoadCase
real(pReal), dimension (3,3) :: rotation = math_I3 !< rotation of BC
type(tBoundaryCondition) :: P, & !< stress BC
@ -129,14 +134,19 @@ program DAMASK_spectral_Driver
character(len=1024) :: incInfo !< string parsed to solution with information about current load case
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState) solres
integer(kind=MPI_OFFSET_KIND) :: my_offset
integer, dimension(:), allocatable :: outputSize
PetscErrorCode :: ierr
external :: quit
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll(temperature = 300.0_pReal, el = 1_pInt, ip = 1_pInt)
mainProcess: if (worldrank == 0) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
!--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases
@ -256,6 +266,7 @@ program DAMASK_spectral_Driver
! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0_pInt
if (worldrank == 0) then
checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
write (loadcase_string, '(i6)' ) currentLoadCase
write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
@ -308,20 +319,19 @@ program DAMASK_spectral_Driver
loadCases(currentLoadCase)%restartfrequency
if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
enddo checkLoadcases
endif
!--------------------------------------------------------------------------------------------------
! doing initialization depending on selected solver
select case (spectral_solver)
case (DAMASK_spectral_SolverBasic_label)
call basic_init(loadCases(1)%temperature)
case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init(loadCases(1)%temperature)
case (DAMASK_spectral_SolverAL_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call AL_init(loadCases(1)%temperature)
case (DAMASK_spectral_SolverPolarisation_label)
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call Polarisation_init(loadCases(1)%temperature)
case default
@ -330,19 +340,15 @@ program DAMASK_spectral_Driver
!--------------------------------------------------------------------------------------------------
! write header of output file
if (appendToOutFile) then ! after restart, append to existing results file
open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED', position='APPEND', status='OLD')
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD')
else ! open new files ...
if (.not. appendToOutFile) then ! after restart, append to existing results file
if (worldrank == 0) then
open(newunit=resUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.spectralOut',form='UNFORMATTED',status='REPLACE')
write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header
write(resUnit) 'workingdir:', trim(getSolverWorkingDirectoryName())
write(resUnit) 'geometry:', trim(geometryFile)
write(resUnit) 'grid:', grid
write(resUnit) 'size:', geomSize
write(resUnit) 'grid:', gridGlobal
write(resUnit) 'size:', geomSizeGlobal
write(resUnit) 'materialpoint_sizeResults:', materialpoint_sizeResults
write(resUnit) 'loadcases:', size(loadCases)
write(resUnit) 'frequencies:', loadCases%outputfrequency ! one entry per currentLoadCase
@ -350,14 +356,36 @@ program DAMASK_spectral_Driver
write(resUnit) 'logscales:', loadCases%logscale
write(resUnit) 'increments:', loadCases%incs ! one entry per currentLoadCase
write(resUnit) 'startingIncrement:', restartInc - 1_pInt ! start with writing out the previous inc
write(resUnit) 'eoh' ! end of header
write(resUnit) materialpoint_results ! initial (non-deformed or read-in) results
write(resUnit) 'eoh'
close(resUnit) ! end of header
endif
endif
allocate(outputSize(worldsize), source = 0_pInt); outputSize(worldrank+1) = sizeof(materialpoint_results)
call MPI_Allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
call MPI_File_open(PETSC_COMM_WORLD, &
trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//'.spectralOut', &
MPI_MODE_WRONLY + MPI_MODE_APPEND, &
MPI_INFO_NULL, &
resUnit, &
ierr)
call MPI_File_get_position(resUnit,my_offset,ierr)
my_offset = my_offset + sum(outputSize(1:worldrank))
call MPI_File_seek (resUnit,my_offset,MPI_SEEK_SET,ierr)
call MPI_File_write(resUnit, materialpoint_results, size(materialpoint_results), &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a)') ' header of result file written out'
flush(6)
if (worldrank == 0) then
if (appendToOutFile) then ! after restart, append to existing results file
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED', position='APPEND', status='OLD')
else ! open new files ...
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
write(6,'(/,a)') ' header of result file written out'
flush(6)
endif
endif
!--------------------------------------------------------------------------------------------------
@ -405,7 +433,6 @@ program DAMASK_spectral_Driver
!--------------------------------------------------------------------------------------------------
! forward solution
select case(spectral_solver)
case (DAMASK_spectral_SolverBasic_label)
case (DAMASK_spectral_SolverBasicPETSC_label)
call BasicPETSC_forward (&
guess,timeinc,timeIncOld,remainingLoadCaseTime, &
@ -430,6 +457,7 @@ program DAMASK_spectral_Driver
!--------------------------------------------------------------------------------------------------
! report begin of new increment
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,a,es12.5'//&
',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
@ -444,17 +472,11 @@ program DAMASK_spectral_Driver
',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
'-',stepFraction, '/', subStepFactor**cutBackLevel
select case(spectral_solver)
endif
!--------------------------------------------------------------------------------------------------
! calculate solution
case (DAMASK_spectral_SolverBasic_label)
solres = basic_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
P_BC = loadCases(currentLoadCase)%P, &
F_BC = loadCases(currentLoadCase)%deformation, &
temperature_bc = loadCases(currentLoadCase)%temperature, &
rotation_BC = loadCases(currentLoadCase)%rotation)
select case(spectral_solver)
case (DAMASK_spectral_SolverBasicPETSC_label)
solres = BasicPETSC_solution (&
incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
@ -486,7 +508,7 @@ program DAMASK_spectral_Driver
cutBack = .False.
if(solres%termIll .or. .not. solres%converged) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back
write(6,'(/,a)') ' cut back detected'
if (worldrank == 0) write(6,'(/,a)') ' cut back detected'
cutBack = .True.
stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1_pInt
@ -505,28 +527,36 @@ program DAMASK_spectral_Driver
else
guess = .true. ! start guessing after first converged (sub)inc
endif
if (.not. cutBack) &
if (.not. cutBack) then
if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
flush(statUnit)
endif
endif
enddo subIncLooping
cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half number of subincs next inc
if(solres%converged) then ! report converged inc
convergedCounter = convergedCounter + 1_pInt
if (worldrank == 0) &
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
' increment ', totalIncsCounter, ' converged'
else
if (worldrank == 0) &
write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') & ! report non-converged inc
' increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt
endif; flush(6)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency
if (worldrank == 0) &
write(6,'(1/,a)') ' ... writing results to file ......................................'
write(resUnit) materialpoint_results ! write result to file
flush(resUnit)
my_offset = my_offset + sum(outputSize)
call MPI_File_seek (resUnit,my_offset,MPI_SEEK_SET,ierr)
call MPI_File_write(resUnit, materialpoint_results, size(materialpoint_results), &
MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
endif
if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. & ! at frequency of writing restart information set restart parameter for FEsolving
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! ToDo first call to CPFEM_general will write?
mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then ! first call to CPFEM_general will write?
restartWrite = .true.
lastRestartWritten = inc
endif
@ -538,9 +568,20 @@ program DAMASK_spectral_Driver
enddo incLooping
enddo loadCaseLooping
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
if (worldrank == 0) then
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
endif
call MPI_file_close(resUnit,ierr)
close(statUnit)
select case (spectral_solver)
case (DAMASK_spectral_SolverBasic_label)
call basic_destroy()
case (DAMASK_spectral_SolverBasicPETSC_label)
call BasicPETSC_destroy()
case (DAMASK_spectral_SolverAL_label)
@ -549,16 +590,6 @@ program DAMASK_spectral_Driver
call Polarisation_destroy()
end select
!--------------------------------------------------------------------------------------------------
! report summary of whole calculation
write(6,'(/,a)') ' ###########################################################################'
write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
notConvergedCounter + convergedCounter, ' (', &
real(convergedCounter, pReal)/&
real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
' %) increments converged!'
close(resUnit)
close(statUnit)
if (notConvergedCounter > 0_pInt) call quit(3_pInt) ! error if some are not converged
call quit(0_pInt) ! no complains ;)
@ -576,11 +607,14 @@ end program DAMASK_spectral_Driver
subroutine quit(stop_id)
use prec, only: &
pInt
use numerics, only: &
worldrank
implicit none
integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
if (worldrank == 0_pInt) then
call date_and_time(values = dateAndTime)
write(6,'(/,a)') 'DAMASK terminated on:'
write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
@ -589,8 +623,11 @@ subroutine quit(stop_id)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
endif
if (stop_id == 0_pInt) stop 0 ! normal termination
if (stop_id < 0_pInt) then ! trigger regridding
if (worldrank == 0_pInt) &
write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
stop 2
endif

185
code/DAMASK_spectral_solverAL.f90
View File

@ -45,11 +45,9 @@ module DAMASK_spectral_solverAL
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients
F_lastInc2, & !< field of 2nd previous compatible deformation gradients
F_lambda_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient
F_lambdaDot !< field of assumed rate of incopatible deformation gradient
complex(pReal),private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
@ -118,19 +116,21 @@ subroutine AL_init(temperature)
debug_spectralRestart
use FEsolving, only: &
restartInc
use numerics, only: &
worldrank, &
worldsize
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
grid, &
Utilities_updateIPcoords, &
grid1Red, &
geomSize, &
wgt
use mesh, only: &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
gridLocal, &
gridGlobal
use math, only: &
math_invSym3333
@ -144,30 +144,41 @@ subroutine AL_init(temperature)
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
character(len=1024) :: rankStr
call Utilities_init()
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lastInc2 (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lambda_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lambdaDot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lambda_lastInc(3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lambdaDot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
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, &
DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr)
@ -183,33 +194,31 @@ subroutine AL_init(temperature)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity
F_lastInc2 = F_lastInc
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity
F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_lambda = F
F_lambda_lastInc = F_lastInc
elseif (restartInc > 1_pInt) then ! read in F to calculate coordinates and initialize CPFEM general
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6)
call IO_read_realFile(777,'F', trim(getSolverJobName()),size(F))
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr), trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_realFile(777,'F_lastInc', trim(getSolverJobName()),size(F_lastInc))
call IO_read_realFile(777,'F_lastInc'//trim(rankStr), trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
call IO_read_realFile(777,'F_lastInc2', trim(getSolverJobName()),size(F_lastInc2))
read (777,rec=1) F_lastInc2
close (777)
call IO_read_realFile(777,'F_lambda',trim(getSolverJobName()),size(F_lambda))
call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda))
read (777,rec=1) F_lambda
close (777)
call IO_read_realFile(777,'F_lambda_lastInc',&
call IO_read_realFile(777,'F_lambda_lastInc'//trim(rankStr),&
trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F_lambda_lastInc
close (777)
if (worldrank == 0_pInt) then
call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim))
read (777,rec=1) F_aim
close (777)
@ -220,16 +229,16 @@ subroutine AL_init(temperature)
read (777,rec=1) f_aimDot
close (777)
endif
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call Utilities_constitutiveResponse(F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]),&
call utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc,F,&
temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
nullify(F)
nullify(F_lambda)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
if (restartInc > 1_pInt) then ! using old values from files
if (restartInc > 1_pInt .and. worldrank == 0_pInt) then ! using old values from files
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -256,7 +265,8 @@ end subroutine AL_init
!> @brief solution for the AL scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function &
AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, &
rotation_BC,density)
use numerics, only: &
update_gamma
use math, only: &
@ -342,7 +352,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
itmax, &
itmin, &
polarAlpha, &
polarBeta
polarBeta, &
worldrank
use mesh, only: &
gridLocal
use IO, only: &
IO_intOut
use math, only: &
@ -353,11 +366,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_mul33x33, &
PI
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
wgt, &
field_real, &
field_fourier, &
field_realMPI, &
field_fourierMPI, &
Utilities_FFTforward, &
Utilities_fourierConvolution, &
Utilities_inverseLaplace, &
@ -371,6 +382,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
debug_spectralRotation
use homogenization, only: &
materialpoint_dPdF
use FEsolving, only: &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
@ -410,12 +423,14 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
if(totalIter <= PETScIter) then ! new iteration
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
if (worldrank == 0_pInt) then
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
@ -423,30 +438,15 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
endif
flush(6)
endif
!--------------------------------------------------------------------------------------------------
! evaluate inertia
dynamic: if (params%density > 0.0_pReal) then
residual_F = ((F - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/&
((params%timeinc + params%timeincOld)/2.0_pReal)
residual_F = params%density*product(geomSize/grid)*residual_F
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
call Utilities_inverseLaplace()
inertiaField_fourier = field_fourier
else dynamic
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
endif dynamic
!--------------------------------------------------------------------------------------------------
!
field_real = 0.0_pReal
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
field_real(i,j,k,1:3,1:3) = &
field_realMPI = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
field_realMPI(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))
@ -456,26 +456,25 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call Utilities_FFTforward()
field_fourier = field_fourier + polarAlpha*inertiaField_fourier
call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_lambda = polarBeta*F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),&
[3,3,grid(1),grid(2),grid(3)],order=[3,4,5,1,2])
residual_F_lambda = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
P_avLastEval = P_av
call Utilities_constitutiveResponse(F_lastInc,F - residual_F_lambda/polarBeta,params%temperature,params%timeinc, &
residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .False.
!--------------------------------------------------------------------------------------------------
! calculate divergence
field_real = 0.0_pReal
field_real = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3])
field_realMPI = 0.0_pReal
field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_FFTbackward()
@ -483,7 +482,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0_pInt
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
e = e + 1_pInt
residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - &
@ -494,8 +493,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! calculating curl
field_real = 0.0_pReal
field_real = reshape(F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3])
field_realMPI = 0.0_pReal
field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
call Utilities_FFTforward()
err_curl = Utilities_curlRMS()
call Utilities_FFTbackward()
@ -515,7 +514,8 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
err_curl_tolRel, &
err_curl_tolAbs, &
err_stress_tolAbs, &
err_stress_tolRel
err_stress_tolRel, &
worldrank
use math, only: &
math_mul3333xx33
use FEsolving, only: &
@ -562,6 +562,7 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
!--------------------------------------------------------------------------------------------------
! report
if (worldrank == 0_pInt) then
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')'
@ -571,6 +572,7 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr
err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
endif
end subroutine AL_converged
@ -583,16 +585,16 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
math_mul3333xx33, &
math_transpose33, &
math_rotate_backward33
use numerics, only: &
worldrank
use mesh, only: &
gridLocal
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use mesh, only: &
mesh_ipCoordinates,&
mesh_deformedCoordsFFT
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
@ -613,6 +615,7 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda
integer(pInt) :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
character(len=1024) :: rankStr
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
@ -620,23 +623,24 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if (restartWrite) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
call IO_write_jobRealFile(777,'F',size(F)) ! writing deformation gradient field to file
endif
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobRealFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
call IO_write_jobRealFile(777,'F_lastInc2',size(F_lastInc2)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc2
close (777)
call IO_write_jobRealFile(777,'F_lambda',size(F_lambda)) ! writing deformation gradient field to file
call IO_write_jobRealFile(777,'F_lambda'//trim(rankStr),size(F_lambda)) ! writing deformation gradient field to file
write (777,rec=1) F_lambda
close (777)
call IO_write_jobRealFile(777,'F_lambda_lastInc',size(F_lambda_lastInc)) ! writing F_lastInc field to file
call IO_write_jobRealFile(777,'F_lambda_lastInc'//trim(rankStr),size(F_lambda_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lambda_lastInc
close (777)
if (worldrank == 0_pInt) then
call IO_write_jobRealFile(777,'F_aim',size(F_aim))
write (777,rec=1) F_aim
close(777)
@ -652,15 +656,15 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
call IO_write_jobRealFile(777,'C_volAvgLastInc',size(C_volAvgLastInc))
write (777,rec=1) C_volAvgLastInc
close(777)
endif
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid(3)])
F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)])
F_lambda = reshape(F_lambda_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)])
C_volAvg = C_volAvgLastInc
else
ForwardData = .True.
@ -679,25 +683,24 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
timeinc_old,guess,F_lastInc,reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc2 = F_lastInc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])
timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,gridLocal(1), &
gridLocal(2),gridLocal(3)]))
F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
endif
F_aim = F_aim + f_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)), &
[9,grid(1),grid(2),grid(3)])
[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), &
[9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition
[9,gridLocal(1),gridLocal(2),gridLocal(3)]) ! does not have any average value as boundary condition
if (.not. guess) then ! large strain forwarding
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,&

148
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -42,8 +42,7 @@ module DAMASK_spectral_SolverBasicPETSc
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot, F_lastInc2
complex(pReal), private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
@ -105,6 +104,9 @@ subroutine basicPETSc_init(temperature)
debug_spectralRestart
use FEsolving, only: &
restartInc
use numerics, only: &
worldrank, &
worldsize
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
@ -112,10 +114,12 @@ subroutine basicPETSc_init(temperature)
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
utilities_updateIPcoords, &
grid, &
grid1Red, &
wgt, &
geomSize
wgt
use mesh, only: &
gridLocal, &
gridGlobal, &
mesh_ipCoordinates
use math, only: &
math_invSym3333
@ -128,31 +132,38 @@ subroutine basicPETSc_init(temperature)
PetscObject :: dummy
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
character(len=1024) :: rankStr
call Utilities_init()
mainProcess: if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif mainProcess
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lastInc2(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
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), & ! overall grid
PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & ! domain decomposition strategy (or local (per core) grid)
9,1, & ! #dof (F tensor), ghost boundary width (domain overlap)
PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! todo
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
1, 1, worldsize, &
9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
@ -168,33 +179,31 @@ subroutine basicPETSc_init(temperature)
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_lastInc2 = F_lastInc
F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity
F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
elseif (restartInc > 1_pInt) then ! using old values from file
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file'
flush(6)
call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F))
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_realFile(777,'F_lastInc',trim(getSolverJobName()),size(F_lastInc))
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
call IO_read_realFile(777,'F_lastInc2',trim(getSolverJobName()),size(F_lastInc2))
read (777,rec=1) F_lastInc2
close (777)
if (worldrank == 0_pInt) then
call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
read (777,rec=1) f_aimDot
close (777)
endif
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
endif
call utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc, &
reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]), &
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc, F, &
temperature, &
0.0_pReal, &
P, &
@ -204,7 +213,7 @@ subroutine basicPETSc_init(temperature)
math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
if (restartInc > 1_pInt) then ! using old values from files
if (restartInc > 1_pInt .and. worldrank == 0_pInt) then ! using old values from files
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -231,7 +240,8 @@ type(tSolutionState) function basicPETSc_solution( &
incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
use numerics, only: &
update_gamma, &
itmax
itmax, &
worldrank
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
Utilities_maskedCompliance, &
@ -305,6 +315,10 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
itmax, &
itmin
use numerics, only: &
worldrank
use mesh, only: &
gridLocal
use math, only: &
math_rotate_backward33, &
math_transpose33, &
@ -314,11 +328,9 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
debug_spectral, &
debug_spectralRotation
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
wgt, &
field_real, &
field_fourier, &
field_realMPI, &
field_fourierMPI, &
Utilities_FFTforward, &
Utilities_FFTbackward, &
Utilities_fourierConvolution, &
@ -327,12 +339,17 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
Utilities_divergenceRMS
use IO, only: &
IO_intOut
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: &
x_scal, &
PetscScalar, dimension(3,3, &
XG_RANGE,YG_RANGE,ZG_RANGE) :: &
x_scal
PetscScalar, dimension(3,3, &
X_RANGE,Y_RANGE,Z_RANGE) :: &
f_scal
PetscInt :: &
PETScIter, &
@ -348,6 +365,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
if (worldrank == 0_pInt) then
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
@ -355,29 +373,15 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
flush(6)
endif
!--------------------------------------------------------------------------------------------------
! evaluate inertia
if (params%density > 0.0_pReal) then
f_scal = ((x_scal - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/&
((params%timeinc + params%timeincOld)/2.0_pReal)
f_scal = params%density*product(geomSize/grid)*f_scal
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
call Utilities_inverseLaplace()
inertiaField_fourier = field_fourier
else
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
flush(6)
endif
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, &
f_scal,C_volAvg,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
@ -388,18 +392,16 @@ if (params%density > 0.0_pReal) then
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(f_scal,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
field_realMPI = 0.0_pReal
field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = f_scal
call Utilities_FFTforward()
field_fourier = field_fourier + inertiaField_fourier
err_divDummy = Utilities_divergenceRMS()
err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC))
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(f_scal),order=[3,4,5,1,2])
f_scal = field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
end subroutine BasicPETSc_formResidual
@ -414,7 +416,8 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
err_div_tolRel, &
err_div_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
err_stress_tolAbs, &
worldrank
use FEsolving, only: &
terminallyIll
@ -434,7 +437,6 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs)
stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
err_divPrev = err_div; err_div = err_divDummy
converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
@ -449,6 +451,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
!--------------------------------------------------------------------------------------------------
! report
if (worldrank == 0_pInt) then
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,')'
@ -456,6 +459,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du
err_stress/stressTol, ' (',err_stress, ' Pa, tol =',stressTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
endif
end subroutine BasicPETSc_converged
@ -466,9 +470,9 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
gridLocal
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
Utilities_calculateRate, &
Utilities_forwardField, &
utilities_updateIPcoords, &
@ -478,6 +482,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
use numerics, only: &
worldrank
implicit none
real(pReal), intent(in) :: &
@ -492,22 +498,24 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
guess
PetscScalar, pointer :: F(:,:,:,:)
PetscErrorCode :: ierr
character(len=1024) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
call IO_write_jobRealFile(777,'F',size(F)) ! writing deformation gradient field to file
endif
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobRealFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
call IO_write_jobRealFile(777,'F_lastInc2',size(F_lastInc2)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc2
close (777)
if (worldrank == 0_pInt) then
call IO_write_jobRealFile(777,'F_aimDot',size(F_aimDot))
write (777,rec=1) F_aimDot
close(777)
@ -518,11 +526,12 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
write (777,rec=1) C_volAvgLastInc
close(777)
endif
endif
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)])
F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)])
C_volAvg = C_volAvgLastInc
else
ForwardData = .True.
@ -543,16 +552,15 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
! update coordinates and rate and forward last inc
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc2 = F_lastInc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)])
timeinc_old,guess,F_lastInc,reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)])
math_rotate_backward33(F_aim,rotation_BC)),[9,gridLocal(1),gridLocal(2),gridLocal(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine BasicPETSc_forward

189
code/DAMASK_spectral_solverPolarisation.f90
View File

@ -45,11 +45,9 @@ module DAMASK_spectral_solverPolarisation
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients
F_lastInc2, & !< field of 2nd previous compatible deformation gradients
F_tau_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient
F_tauDot !< field of assumed rate of incopatible deformation gradient
complex(pReal),private, dimension(:,:,:,:,:), allocatable :: inertiaField_fourier
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
@ -118,19 +116,21 @@ subroutine Polarisation_init(temperature)
debug_spectralRestart
use FEsolving, only: &
restartInc
use numerics, only: &
worldrank, &
worldsize
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
grid, &
Utilities_updateIPcoords, &
grid1Red, &
geomSize, &
wgt
use mesh, only: &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
gridLocal, &
gridGlobal
use math, only: &
math_invSym3333
@ -144,30 +144,41 @@ subroutine Polarisation_init(temperature)
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_tau
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
character(len=1024) :: rankStr
call Utilities_init()
if (worldrank == 0_pInt) then
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
endif
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_lastInc2 (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (F_tauDot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (inertiaField_fourier (grid1Red,grid(2),grid(3),3,3),source = cmplx(0.0_pReal,0.0_pReal,pReal))
allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_tau_lastInc(3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_tauDot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
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, &
DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,dummy,ierr)
@ -183,32 +194,30 @@ subroutine Polarisation_init(temperature)
F => xx_psc(0:8,:,:,:)
F_tau => xx_psc(9:17,:,:,:)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity
F_lastInc2 = F_lastInc
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity
F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_tau = 2.0_pReal* F
F_tau_lastInc = 2.0_pReal*F_lastInc
elseif (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
flush(6)
call IO_read_realFile(777,'F',trim(getSolverJobName()),size(F))
write(rankStr,'(a1,i0)')'_',worldrank
call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F))
read (777,rec=1) F
close (777)
call IO_read_realFile(777,'F_lastInc',trim(getSolverJobName()),size(F_lastInc))
call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc
close (777)
call IO_read_realFile(777,'F_lastInc2',trim(getSolverJobName()),size(F_lastInc2))
read (777,rec=1) F_lastInc2
close (777)
call IO_read_realFile(777,'F_tau',trim(getSolverJobName()),size(F_tau))
call IO_read_realFile(777,'F_tau'//trim(rankStr),trim(getSolverJobName()),size(F_tau))
read (777,rec=1) F_tau
close (777)
call IO_read_realFile(777,'F_tau_lastInc',&
call IO_read_realFile(777,'F_tau_lastInc'//trim(rankStr),&
trim(getSolverJobName()),size(F_tau_lastInc))
read (777,rec=1) F_tau_lastInc
close (777)
if (worldrank == 0_pInt) then
call IO_read_realFile(777,'F_aim', trim(getSolverJobName()),size(F_aim))
read (777,rec=1) F_aim
close (777)
@ -219,16 +228,16 @@ subroutine Polarisation_init(temperature)
read (777,rec=1) f_aimDot
close (777)
endif
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call Utilities_constitutiveResponse(F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]),&
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc,F,&
temperature,0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
nullify(F)
nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
if (restartInc > 1_pInt) then ! using old values from files
if (restartInc > 1_pInt .and. worldrank == 0_pInt) then ! using old values from files
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -255,7 +264,8 @@ end subroutine Polarisation_init
!> @brief solution for the Polarisation scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function &
Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC,density)
Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, &
rotation_BC,density)
use numerics, only: &
update_gamma
use math, only: &
@ -267,6 +277,8 @@ type(tSolutionState) function &
use FEsolving, only: &
restartWrite, &
terminallyIll
use numerics, only: &
worldrank
implicit none
@ -340,9 +352,12 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
itmax, &
itmin, &
polarAlpha, &
polarBeta
polarBeta, &
worldrank
use IO, only: &
IO_intOut
use mesh, only: &
gridLocal
use math, only: &
math_rotate_backward33, &
math_transpose33, &
@ -351,11 +366,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_mul33x33, &
PI
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
wgt, &
field_real, &
field_fourier, &
field_realMPI, &
field_fourierMPI, &
Utilities_FFTforward, &
Utilities_fourierConvolution, &
Utilities_inverseLaplace, &
@ -369,6 +382,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
debug_spectralRotation
use homogenization, only: &
materialpoint_dPdF
use FEsolving, only: &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
@ -408,12 +423,14 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
if (totalIter <= PETScIter) then ! new iteration
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
if (worldrank == 0_pInt) then
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
@ -421,30 +438,15 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
endif
flush(6)
endif
!--------------------------------------------------------------------------------------------------
! evaluate inertia
dynamic: if (params%density > 0.0_pReal) then
residual_F = ((F - F_lastInc)/params%timeinc - (F_lastInc - F_lastInc2)/params%timeincOld)/&
((params%timeinc + params%timeincOld)/2.0_pReal)
residual_F = params%density*product(geomSize/grid)*residual_F
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
call Utilities_inverseLaplace()
inertiaField_fourier = field_fourier
else dynamic
inertiaField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
endif dynamic
!--------------------------------------------------------------------------------------------------
!
field_real = 0.0_pReal
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
field_real(i,j,k,1:3,1:3) = &
field_realMPI = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
field_realMPI(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
@ -453,26 +455,25 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call Utilities_FFTforward()
field_fourier = field_fourier + polarAlpha*inertiaField_fourier
call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = polarBeta*F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),&
[3,3,grid(1),grid(2),grid(3)],order=[3,4,5,1,2])
residual_F_tau = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
P_avLastEval = P_av
call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%temperature,params%timeinc, &
residual_F,C_volAvg,C_minMaxAvg,P_av,ForwardData,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .False.
!--------------------------------------------------------------------------------------------------
! calculate divergence
field_real = 0.0_pReal
field_real = reshape(residual_F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3])
field_realMPI = 0.0_pReal
field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_FFTbackward()
@ -480,7 +481,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0_pInt
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
e = e + 1_pInt
residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
@ -491,8 +492,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! calculating curl
field_real = 0.0_pReal
field_real = reshape(F,[grid(1),grid(2),grid(3),3,3],order=[4,5,1,2,3])
field_realMPI = 0.0_pReal
field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
call Utilities_FFTforward()
err_curl = Utilities_curlRMS()
call Utilities_FFTbackward()
@ -512,7 +513,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
err_curl_tolRel, &
err_curl_tolAbs, &
err_stress_tolabs, &
err_stress_tolrel
err_stress_tolrel, &
worldrank
use math, only: &
math_mul3333xx33
use FEsolving, only: &
@ -559,6 +561,7 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
!--------------------------------------------------------------------------------------------------
! report
if (worldrank == 0_pInt) then
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol =',curlTol,')'
@ -568,6 +571,7 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
err_BC/BC_tol, ' (',err_BC, ' Pa, tol =',BC_tol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
endif
end subroutine Polarisation_converged
@ -581,19 +585,19 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
math_transpose33, &
math_rotate_backward33
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use mesh, only: &
mesh_ipCoordinates,&
mesh_deformedCoordsFFT
gridLocal
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
use numerics, only: &
worldrank
implicit none
real(pReal), intent(in) :: &
@ -610,6 +614,7 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau
integer(pInt) :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
character(len=1024) :: rankStr
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
@ -617,23 +622,22 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
F => xx_psc(0:8,:,:,:)
F_tau => xx_psc(9:17,:,:,:)
if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart'
if (worldrank == 0_pInt) write(6,'(/,a)') ' writing converged results for restart'
flush(6)
call IO_write_jobRealFile(777,'F',size(F)) ! writing deformation gradient field to file
write(rankStr,'(a1,i0)')'_',worldrank
call IO_write_jobRealFile(777,'F'//trim(rankStr),size(F)) ! writing deformation gradient field to file
write (777,rec=1) F
close (777)
call IO_write_jobRealFile(777,'F_lastInc',size(F_lastInc)) ! writing F_lastInc field to file
call IO_write_jobRealFile(777,'F_lastInc'//trim(rankStr),size(F_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc
close (777)
call IO_write_jobRealFile(777,'F_lastInc2',size(F_lastInc2)) ! writing F_lastInc field to file
write (777,rec=1) F_lastInc2
close (777)
call IO_write_jobRealFile(777,'F_tau',size(F_tau)) ! writing deformation gradient field to file
call IO_write_jobRealFile(777,'F_tau'//trim(rankStr),size(F_tau)) ! writing deformation gradient field to file
write (777,rec=1) F_tau
close (777)
call IO_write_jobRealFile(777,'F_tau_lastInc',size(F_tau_lastInc)) ! writing F_lastInc field to file
call IO_write_jobRealFile(777,'F_tau_lastInc'//trim(rankStr),size(F_tau_lastInc)) ! writing F_lastInc field to file
write (777,rec=1) F_tau_lastInc
close (777)
if (worldrank == 0_pInt) then
call IO_write_jobRealFile(777,'F_aim',size(F_aim))
write (777,rec=1) F_aim
close(777)
@ -650,13 +654,14 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
write (777,rec=1) C_volAvgLastInc
close(777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
endif
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)])
F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)])
F_tau= reshape(F_tau_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)])
C_volAvg = C_volAvgLastInc
else
ForwardData = .True.
@ -675,25 +680,27 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
timeinc_old,guess,F_lastInc, &
reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), &
timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc2 = F_lastInc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])
timeinc_old,guess,F_tau_lastInc, &
reshape(F_tau,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
F_tau_lastInc = reshape(F_tau,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)])
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition
math_rotate_backward33(F_aim,rotation_BC)), &
[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & ! does not have any average value as boundary condition
[9,gridLocal(1),gridLocal(2),gridLocal(3)])
if (.not. guess) then ! large strain forwarding
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,&

672
code/DAMASK_spectral_utilities.f90
View File

File diff suppressed because it is too large Load Diff

11
code/Makefile
View File

@ -336,7 +336,7 @@ SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o
FEsolving.o mesh.o material.o lattice.o \
$(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(PLASTIC_FILES) \
crystallite.o $(HOMOGENIZATION_FILES) CPFEM.o \
DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o \
DAMASK_spectral_utilities.o \
DAMASK_spectral_solverAL.o DAMASK_spectral_solverBasicPETSc.o DAMASK_spectral_solverPolarisation.o
DAMASK_spectral.exe: DAMASK_spectral_driver.o
@ -345,7 +345,7 @@ DAMASK_spectral.exe: DAMASK_spectral_driver.o
$(SPECTRAL_FILES) $(LIBRARIES) $(SUFFIX)
DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 DAMASK_spectral_solverBasic.o \
DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 \
DAMASK_spectral_solverAL.o \
DAMASK_spectral_solverBasicPETSc.o \
DAMASK_spectral_solverPolarisation.o
@ -354,9 +354,6 @@ DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 DAMASK_spectral_solverBasic
DAMASK_spectral_solverAL.o: DAMASK_spectral_solverAL.f90 \
DAMASK_spectral_utilities.o
DAMASK_spectral_solverBasic.o: DAMASK_spectral_solverBasic.f90 \
DAMASK_spectral_utilities.o
DAMASK_spectral_solverPolarisation.o: DAMASK_spectral_solverPolarisation.f90 \
DAMASK_spectral_utilities.o
@ -370,8 +367,8 @@ DAMASK_spectral_utilities.o: DAMASK_spectral_utilities.f90 \
# FEM Solver
#####################
VPATH := ../private/FEM/code
DAMASK_FEM.exe: COMPILE += -DFEM -DmultiphysicsOut
DAMASK_FEM.exe: COMPILE_MAXOPTI += -DFEM -DmultiphysicsOut
DAMASK_FEM.exe: COMPILE += -DFEM
DAMASK_FEM.exe: COMPILE_MAXOPTI += -DFEM
DAMASK_FEM.exe: MESHNAME := ../private/FEM/code/meshFEM.f90
DAMASK_FEM.exe: INTERFACENAME := ../private/FEM/code/DAMASK_FEM_interface.f90
DAMASK_FEM.exe: INCLUDE_DIRS += -I./

163
code/mesh.f90
View File

@ -15,7 +15,7 @@ module mesh
implicit none
private
integer(pInt), public, protected :: &
mesh_NcpElems, & !< total number of CP elements in mesh
mesh_NcpElems, & !< total number of CP elements in local mesh
mesh_NelemSets, &
mesh_maxNelemInSet, &
mesh_Nmaterials, &
@ -29,6 +29,18 @@ module mesh
mesh_maxNcellnodes, & !< max number of cell nodes in any CP element
mesh_Nelems !< total number of elements in mesh
#ifdef Spectral
integer(pInt), public, protected :: &
mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh
gridGlobal(3), &
gridLocal (3), &
gridOffset
real(pReal), public, protected :: &
geomSizeGlobal(3), &
geomSizeLocal (3), &
geomSizeOffset
#endif
integer(pInt), dimension(:,:), allocatable, public, protected :: &
mesh_element, & !< FEid, type(internal representation), material, texture, node indices as CP IDs
mesh_sharedElem, & !< entryCount and list of elements containing node
@ -103,10 +115,15 @@ module mesh
logical, private :: noPart !< for cases where the ABAQUS input file does not use part/assembly information
#endif
#ifdef Spectral
include 'fftw3.f03'
#ifdef PETSc
#include <petsc-finclude/petscsys.h>
#endif
#ifdef Spectral
include 'fftw3-mpi.f03'
#endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
! Hence, I suggest to prefix with "FE_"
@ -470,6 +487,9 @@ contains
!! Order and routines strongly depend on type of solver
!--------------------------------------------------------------------------------------------------
subroutine mesh_init(ip,el)
#ifdef Spectral
use, intrinsic :: iso_c_binding
#endif
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: &
@ -503,6 +523,9 @@ subroutine mesh_init(ip,el)
modelName
implicit none
#ifdef Spectral
integer(C_INTPTR_T) :: gridMPI(3), alloc_local, local_K, local_K_offset
#endif
integer(pInt), parameter :: FILEUNIT = 222_pInt
integer(pInt), intent(in) :: el, ip
integer(pInt) :: j
@ -540,8 +563,25 @@ subroutine mesh_init(ip,el)
myDebug = (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt)
#ifdef Spectral
call fftw_mpi_init()
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
gridGlobal = mesh_spectral_getGrid(fileUnit)
gridMPI = gridGlobal
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, &
MPI_COMM_WORLD, local_K, local_K_offset)
gridLocal(1) = gridGlobal(1)
gridLocal(2) = gridGlobal(2)
gridLocal(3) = local_K
gridOffset = local_K_offset
geomSizeGlobal = mesh_spectral_getSize(fileUnit)
geomSizeLocal(1) = geomSizeGlobal(1)
geomSizeLocal(2) = geomSizeGlobal(2)
geomSizeLocal(3) = geomSizeGlobal(3)*real(gridLocal(3))/real(gridGlobal(3))
geomSizeOffset = geomSizeGlobal(3)*real(gridOffset) /real(gridGlobal(3))
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
call mesh_spectral_count(FILEUNIT)
if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
call mesh_spectral_mapNodesAndElems
@ -656,10 +696,12 @@ subroutine mesh_init(ip,el)
#endif
close (FILEUNIT)
if (worldrank == 0_pInt) then
call mesh_tell_statistics
call mesh_write_meshfile
call mesh_write_cellGeom
call mesh_write_elemGeom
endif
if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
call IO_error(600_pInt) ! ping-pong must be disabled when having non-DAMASK elements
@ -1204,12 +1246,12 @@ subroutine mesh_spectral_count(fileUnit)
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), dimension(3) :: grid
grid = mesh_spectral_getGrid(fileUnit)
mesh_Nelems = product(grid)
mesh_NcpElems = mesh_Nelems
mesh_Nnodes = product(grid+1_pInt)
mesh_Nelems = product(gridLocal)
mesh_NcpElems= mesh_Nelems
mesh_Nnodes = product(gridLocal + 1_pInt)
mesh_NcpElemsGlobal = product(gridGlobal)
end subroutine mesh_spectral_count
@ -1262,27 +1304,18 @@ subroutine mesh_spectral_build_nodes(fileUnit)
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt) :: n
integer(pInt), dimension(3) :: grid
real(pReal), dimension(3) :: geomSize
integer(pInt) :: n, i, j, k
allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal)
allocate (mesh_node (3,mesh_Nnodes), source = 0.0_pReal)
grid = mesh_spectral_getGrid(fileUnit)
geomSize = mesh_spectral_getSize(fileUnit)
forall (n = 0_pInt:mesh_Nnodes-1_pInt)
mesh_node0(1,n+1_pInt) = mesh_unitlength * &
geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) &
/ real(grid(1),pReal)
mesh_node0(2,n+1_pInt) = mesh_unitlength * &
geomSize(2)*real(mod(n/(grid(1)+1_pInt),(grid(2)+1_pInt)),pReal) &
/ real(grid(2),pReal)
mesh_node0(3,n+1_pInt) = mesh_unitlength * &
geomSize(3)*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid(3)+1_pInt)),pReal) &
/ real(grid(3),pReal)
end forall
n = 0_pInt
do k = 1, gridLocal(3); do j = 1, gridLocal(2); do i = 1, gridLocal(1)
n = n + 1_pInt
mesh_node0(1,n) = real(i)*geomSizeLocal(1)/real(gridLocal(1))
mesh_node0(2,n) = real(j)*geomSizeLocal(2)/real(gridLocal(2))
mesh_node0(3,n) = real(k)*geomSizeLocal(3)/real(gridLocal(3)) + geomSizeOffset
enddo; enddo; enddo
mesh_node = mesh_node0
@ -1315,11 +1348,11 @@ subroutine mesh_spectral_build_elements(fileUnit)
headerLength = 0_pInt, &
maxIntCount, &
homog, &
elemType
elemType, &
elemOffset
integer(pInt), dimension(:), allocatable :: &
microstructures
integer(pInt), dimension(3) :: &
grid
microstructures, &
mesh_microGlobal
integer(pInt), dimension(1,1) :: &
dummySet = 0_pInt
character(len=65536) :: &
@ -1328,7 +1361,6 @@ subroutine mesh_spectral_build_elements(fileUnit)
character(len=64), dimension(1) :: &
dummyName = ''
grid = mesh_spectral_getGrid(fileUnit)
homog = mesh_spectral_getHomogenization(fileUnit)
!--------------------------------------------------------------------------------------------------
@ -1359,6 +1391,7 @@ subroutine mesh_spectral_build_elements(fileUnit)
enddo
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems), source = 0_pInt)
allocate (microstructures (1_pInt+maxIntCount), source = 1_pInt)
allocate (mesh_microGlobal(mesh_NcpElemsGlobal), source = 1_pInt)
!--------------------------------------------------------------------------------------------------
! read in microstructures
@ -1367,31 +1400,39 @@ subroutine mesh_spectral_build_elements(fileUnit)
read(fileUnit,'(a65536)') line
enddo
elemType = FE_mapElemtype('C3D8R')
e = 0_pInt
do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!)
do while (e < mesh_NcpElemsGlobal .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!)
microstructures = IO_continuousIntValues(fileUnit,maxIntCount,dummyName,dummySet,0_pInt) ! get affected elements
do i = 1_pInt,microstructures(1_pInt)
e = e+1_pInt ! valid element entry
mesh_microGlobal(e) = microstructures(1_pInt+i)
enddo
enddo
elemType = FE_mapElemtype('C3D8R')
elemOffset = gridLocal(1)*gridLocal(2)*gridOffset
e = 0_pInt
do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!)
e = e+1_pInt ! valid element entry
mesh_element( 1,e) = e ! FE id
mesh_element( 2,e) = elemType ! elem type
mesh_element( 3,e) = homog ! homogenization
mesh_element( 4,e) = microstructures(1_pInt+i) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + &
((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node
mesh_element( 4,e) = mesh_microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/gridLocal(1) + &
((e-1_pInt)/(gridLocal(1)*gridLocal(2)))*(gridLocal(1)+1_pInt) ! base node
mesh_element( 6,e) = mesh_element(5,e) + 1_pInt
mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2_pInt
mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1_pInt
mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1_pInt) * (grid(2) + 1_pInt) ! second floor base node
mesh_element( 7,e) = mesh_element(5,e) + gridLocal(1) + 2_pInt
mesh_element( 8,e) = mesh_element(5,e) + gridLocal(1) + 1_pInt
mesh_element( 9,e) = mesh_element(5,e) +(gridLocal(1) + 1_pInt) * (gridLocal(2) + 1_pInt) ! second floor base node
mesh_element(10,e) = mesh_element(9,e) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2_pInt
mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + gridLocal(1) + 2_pInt
mesh_element(12,e) = mesh_element(9,e) + gridLocal(1) + 1_pInt
mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) ! needed for statistics
mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e))
enddo
enddo
deallocate(microstructures)
deallocate(mesh_microGlobal)
if (e /= mesh_NcpElems) call IO_error(880_pInt,e)
end subroutine mesh_spectral_build_elements
@ -1409,39 +1450,35 @@ subroutine mesh_spectral_build_ipNeighborhood(fileUnit)
integer(pInt) :: &
x,y,z, &
e
integer(pInt), dimension(3) :: &
grid
allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems),source=0_pInt)
grid = mesh_spectral_getGrid(fileUnit)
e = 0_pInt
do z = 0_pInt,grid(3)-1_pInt
do y = 0_pInt,grid(2)-1_pInt
do x = 0_pInt,grid(1)-1_pInt
do z = 0_pInt,gridLocal(3)-1_pInt
do y = 0_pInt,gridLocal(2)-1_pInt
do x = 0_pInt,gridLocal(1)-1_pInt
e = e + 1_pInt
mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x+1_pInt,grid(1)) &
mesh_ipNeighborhood(1,1,1,e) = z * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
+ modulo(x+1_pInt,gridLocal(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x-1_pInt,grid(1)) &
mesh_ipNeighborhood(1,2,1,e) = z * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
+ modulo(x-1_pInt,gridLocal(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) &
+ modulo(y+1_pInt,grid(2)) * grid(1) &
mesh_ipNeighborhood(1,3,1,e) = z * gridLocal(1) * gridLocal(2) &
+ modulo(y+1_pInt,gridLocal(2)) * gridLocal(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) &
+ modulo(y-1_pInt,grid(2)) * grid(1) &
mesh_ipNeighborhood(1,4,1,e) = z * gridLocal(1) * gridLocal(2) &
+ modulo(y-1_pInt,gridLocal(2)) * gridLocal(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,gridLocal(3)) * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,gridLocal(3)) * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt
@ -1841,7 +1878,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
deallocate(Tstar)
deallocate(TstarNew)
! for the state, we first have to know the size------------------------------------------------------------------
! for the state, we first have to know the size----------------------------------------------------
allocate(sizeStateConst(1,1,mesh_NcpElems))
call IO_read_intFile(FILEUNIT,'sizeStateConst',trim(getSolverJobName()),size(sizeStateConst))
read (FILEUNIT,rec=1) sizeStateConst