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

moved public data res,size and homog from mesh to DAMASK_spectral_utilities (as grid and geomSize)

This commit is contained in:
Martin Diehl 2013-05-08 15:52:29 +00:00
parent 7ef0ba688a
commit 85d4a37d95
9 changed files with 533 additions and 537 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
code/DAMASK_spectral_driver.f90
View File

@ -51,12 +51,13 @@ program DAMASK_spectral_Driver
IO_read_jobBinaryFile, &
IO_write_jobBinaryFile, &
IO_intOut, &
IO_warning
IO_warning, &
IO_timeStamp
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic
use math ! need to include the whole module for FFTW
use mesh, only : &
res, &
geomdim, &
mesh_NcpElems
use CPFEM, only: &
CPFEM_initAll
use FEsolving, only: &
@ -71,10 +72,10 @@ program DAMASK_spectral_Driver
materialpoint_sizeResults, &
materialpoint_results
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
tBoundaryCondition, &
tSolutionState, &
debugGeneral, &
debugDivergence, &
cutBack
use DAMASK_spectral_SolverBasic
#ifdef PETSc
@ -150,8 +151,9 @@ program DAMASK_spectral_Driver
!--------------------------------------------------------------------------------------------------
! init DAMASK (all modules)
call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(/,a)') ' <<<+- DAMASK_spectral_driver init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
!--------------------------------------------------------------------------------------------------
@ -327,7 +329,8 @@ program DAMASK_spectral_Driver
case (DAMASK_spectral_SolverBasicPETSc_label)
call basicPETSc_init(loadCases(1)%temperature)
case (DAMASK_spectral_SolverAL_label)
if(debugDivergence) call IO_warning(42_pInt, ext_msg='debug Divergence')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
call AL_init(loadCases(1)%temperature)
#endif
case default
@ -349,10 +352,10 @@ program DAMASK_spectral_Driver
write(resUnit) 'load', trim(loadCaseFile) ! ... and write header
write(resUnit) 'workingdir', trim(getSolverWorkingDirectoryName())
write(resUnit) 'geometry', trim(geometryFile)
!write(resUnit) 'grid', res
!write(resUnit) 'size', geomdim
write(resUnit) 'resolution', res
write(resUnit) 'dimension', geomdim
!write(resUnit) 'grid', grid
!write(resUnit) 'size', geomSize
write(resUnit) 'resolution', grid
write(resUnit) 'dimension', geomSize
write(resUnit) 'materialpoint_sizeResults', materialpoint_sizeResults
write(resUnit) 'loadcases', size(loadCases)
write(resUnit) 'frequencies', loadCases%outputfrequency ! one entry per currentLoadCase
@ -361,11 +364,12 @@ program DAMASK_spectral_Driver
write(resUnit) 'increments', loadCases%incs ! one entry per currentLoadCase
write(resUnit) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc
write(resUnit) 'eoh' ! end of header
write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results
write(resUnit) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:product(grid)) ! initial (non-deformed or read-in) results
open(newunit=statUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
'.sta',form='FORMATTED',status='REPLACE')
write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded' ! statistics file
if (debugGeneral) write(6,'(/,a)') ' header of result file written out'
if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
write(6,'(/,a)') ' header of result file written out'
flush(6)
endif

6
code/DAMASK_spectral_interface.f90
View File

@ -250,7 +250,7 @@ character(len=1024) function storeWorkingDirectory(workingDirectoryArg,geometryA
endif
if (storeWorkingDirectory(len(trim(storeWorkingDirectory)):len(trim(storeWorkingDirectory))) & ! if path seperator is not given, append it
/= pathSep) storeWorkingDirectory = trim(storeWorkingDirectory)//pathSep
!here check if exists and use chdir!
!> @ToDO here check if exists and use chdir!
else ! using path to geometry file as working dir
if (geometryArg(1:1) == pathSep) then ! absolute path given as command line argument
storeWorkingDirectory = geometryArg(1:scan(geometryArg,pathSep,back=.true.))
@ -383,8 +383,8 @@ function rectifyPath(path)
l = len_trim(path)
rectifyPath = path
do i = l,3,-1
if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) &
rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
if ( rectifyPath(i-1:i) == '.'//pathSep .and. rectifyPath(i-2:i-2) /= '.' ) &
rectifyPath(i-1:l) = rectifyPath(i+1:l)//' '
enddo
!remove ../ and corresponding directory from rectifyPath

127
code/DAMASK_spectral_solverAL.f90
View File

@ -94,22 +94,22 @@ module DAMASK_spectral_solverAL
AL_solution, &
AL_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASetLocalFunction, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASetLocalFunction, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort
contains
@ -119,9 +119,14 @@ contains
subroutine AL_init(temperature)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_intOut, &
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile, &
IO_timeStamp
use debug, only : &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use DAMASK_interface, only: &
@ -130,14 +135,12 @@ subroutine AL_init(temperature)
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
grid, &
geomSize, &
wgt
use numerics, only: &
petsc_options
use mesh, only: &
res, &
geomdim, &
wgt, &
mesh_NcpElems, &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
use math, only: &
@ -148,7 +151,7 @@ subroutine AL_init(temperature)
temperature
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P
real(pReal), dimension(:,:,:,:,:), allocatable :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
real(pReal), dimension(3,3,3,3) :: &
@ -164,18 +167,21 @@ subroutine AL_init(temperature)
write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
write(6,'(a16,a)') ' Current time : ',IO_timeStamp()
#include "compilation_info.f90"
allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) !< @Todo sourced allocation allocate(Fdot,source = F_lastInc)
allocate (F_tau_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_tauDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal) !< @Todo sourced allocation allocate(Fdot,source = F_lastInc)
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)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
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)
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
@ -191,13 +197,14 @@ subroutine AL_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,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity
F_tau_lastInc = F_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
F_tau = F
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(/,a,i6,a)') ' reading values of increment ',&
restartInc - 1_pInt,' from file'
elseif (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
flush(6)
call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F))
@ -230,8 +237,8 @@ subroutine AL_init(temperature)
read (777,rec=1) C_minmaxAvg
close (777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(&
F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,&
temp33_Real,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
@ -263,14 +270,13 @@ type(tSolutionState) function &
math_rotate_backward33, &
math_invSym3333
use mesh, only: &
res, &
geomdim, &
mesh_NcpElems, &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
tBoundaryCondition, &
Utilities_forwardField, &
Utilities_calculateRate, &
@ -342,8 +348,8 @@ use mesh, only: &
if ( cutBack) then
F_aim = F_aim_lastInc
F_tau= reshape(F_tau_lastInc,[9,res(1),res(2),res(3)])
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)])
F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)])
C = C_lastInc
else
C_lastInc = C
@ -359,23 +365,23 @@ use mesh, only: &
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(&
F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
F_tauDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,res(1),res(2),res(3)]))
timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc = reshape(F, [3,3,res(1),res(2),res(3)])
F_tau_lastInc = reshape(F_tau,[3,3,res(1),res(2),res(3)])
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)])
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,res(1),res(2),res(3)])
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,res(1),res(2),res(3)]) ! does not have any average value as boundary condition
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
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
CHKERRQ(ierr)
@ -422,22 +428,25 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
itmin, &
polarAlpha, &
polarBeta
use IO, only: &
IO_intOut
use math, only: &
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33, &
math_invSym3333
use mesh, only: &
res, &
wgt
use DAMASK_spectral_Utilities, only: &
grid, &
wgt, &
field_real, &
Utilities_FFTforward, &
Utilities_fourierConvolution, &
Utilities_FFTbackward, &
Utilities_constitutiveResponse, &
debugRotation
use IO, only: IO_intOut
Utilities_constitutiveResponse
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use homogenization, only: &
materialpoint_P, &
materialpoint_dPdF
@ -490,7 +499,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',reportIter, '≤', itmax
if (debugRotation) &
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
@ -503,7 +512,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
!
field_real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
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) = math_mul3333xx33(C_scale,(polarAlpha + polarBeta)*F(1:3,1:3,i,j,k) - &
(polarAlpha)*F_tau(1:3,1:3,i,j,k))
enddo; enddo; enddo
@ -516,8 +525,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = polarBeta*F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),&
[3,3,res(1),res(2),res(3)],order=[3,4,5,1,2])
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])
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -534,7 +543,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
! constructing residual
n_ele = 0_pInt
err_p = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
n_ele = n_ele + 1_pInt
err_p = err_p + sum((math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) - &
(F_tau(1:3,1:3,i,j,k) - &

61
code/DAMASK_spectral_solverBasic.f90
View File

@ -73,6 +73,10 @@ subroutine basic_init(temperature)
IO_write_JobBinaryFile, &
IO_intOut, &
IO_timeStamp
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use DAMASK_interface, only: &
@ -81,20 +85,17 @@ subroutine basic_init(temperature)
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
use mesh, only: &
res, &
grid, &
wgt, &
geomdim, &
scaledDim, &
geomSize
use mesh, only: &
mesh_ipCoordinates, &
mesh_NcpElems, &
mesh_deformedCoordsFFT
implicit none
real(pReal), intent(inout) :: &
temperature
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
real(pReal), dimension(:,:,:,:,:), allocatable :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
real(pReal), dimension(3,3,3,3) :: &
@ -105,22 +106,23 @@ subroutine basic_init(temperature)
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim
allocate (P (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F (3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lastInc (3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal)
allocate (F_lastInc (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal)
allocate (Fdot (3,3,grid(1), grid(2),grid(3)), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! init fields and average quantities
if (restartInc == 1_pInt) then ! no deformation (no restart)
F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid(3)) ! initialize to identity
F_lastInc = F
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
flush(6)
call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F))
@ -147,8 +149,9 @@ subroutine basic_init(temperature)
read (777,rec=1) temp3333_Real
close (777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,F),[3,1,mesh_NcpElems])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,&
temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C_minmaxAvg
endif
@ -173,15 +176,15 @@ type(tSolutionState) function &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
res,&
geomdim, &
wgt, &
mesh_ipCoordinates,&
mesh_NcpElems, &
mesh_deformedCoordsFFT
use IO, only: &
IO_write_JobBinaryFile, &
IO_intOut
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
field_real, &
@ -194,7 +197,9 @@ type(tSolutionState) function &
Utilities_updateGamma, &
Utilities_constitutiveResponse, &
Utilities_calculateRate, &
debugRotation
grid,&
geomSize, &
wgt
use FEsolving, only: &
restartWrite, &
restartRead, &
@ -224,7 +229,7 @@ type(tSolutionState) function &
real(pReal), dimension(3,3) :: &
F_aim_lastIter, & !< aim of last iteration
P_av
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: P
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal) :: err_div, err_stress
@ -261,7 +266,7 @@ type(tSolutionState) function &
C = C_lastInc
else
C_lastInc = C
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,F),[3,1,mesh_NcpElems])
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,F),[3,1,product(grid)])
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
@ -298,7 +303,7 @@ type(tSolutionState) function &
! report begin of new iteration
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',iter, '≤', itmax
if (debugRotation) &
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab)=', &
math_transpose33(math_rotate_backward33(F_aim,rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
@ -322,13 +327,13 @@ type(tSolutionState) function &
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(P,[res(1),res(2),res(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(P,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,rotation_BC))
call Utilities_FFTbackward()
F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
F = F - reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av)
write(6,'(/,a)') ' =========================================================================='
flush(6)
@ -368,7 +373,7 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
err_stress_tol, &
pAvgDivL2
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html)
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html)
err_stress_tol = max(maxval(abs(pAvgStress))*err_stress_tolrel,err_stress_tolabs)
basic_Converged = all([ err_div/pAvgDivL2/err_div_tol,&

119
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -84,22 +84,22 @@ module DAMASK_spectral_SolverBasicPETSc
basicPETSc_solution ,&
basicPETSc_destroy
external :: &
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASetLocalFunction, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort
VecDestroy, &
DMDestroy, &
DMDACreate3D, &
DMCreateGlobalVector, &
DMDASetLocalFunction, &
PETScFinalize, &
SNESDestroy, &
SNESGetNumberFunctionEvals, &
SNESGetIterationNumber, &
SNESSolve, &
SNESSetDM, &
SNESGetConvergedReason, &
SNESSetConvergenceTest, &
SNESSetFromOptions, &
SNESCreate, &
MPI_Abort
contains
@ -111,7 +111,12 @@ subroutine basicPETSc_init(temperature)
use IO, only: &
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile, &
IO_intOut, &
IO_timeStamp
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use DAMASK_interface, only: &
@ -120,15 +125,12 @@ subroutine basicPETSc_init(temperature)
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
grid, &
wgt, &
geomSize
use numerics, only: &
petsc_options
use mesh, only: &
res, &
wgt, &
geomdim, &
scaledDim, &
mesh_NcpElems, &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
use math, only: &
@ -139,7 +141,7 @@ subroutine basicPETSc_init(temperature)
temperature
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
real(pReal), dimension(:,:,:,:,:), allocatable :: P
PetscScalar, dimension(:,:,:,:), pointer :: F
PetscErrorCode :: ierr
PetscObject :: dummy
@ -153,19 +155,19 @@ subroutine basicPETSc_init(temperature)
write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
write(6,'(a,3(f12.5)/)') ' scaledDim x y z:', scaledDim
allocate (P (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc(3,3,res(1),res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3,res(1),res(2),res(3)), source = 0.0_pReal)
allocate (F_lastInc(3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid(3)),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
DMDA_STENCIL_BOX,grid(1),grid(2),grid(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
9,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
@ -180,11 +182,12 @@ 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,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
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)])
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(/,a,i6,a)') ' reading values of increment ',&
restartInc - 1_pInt,' from file'
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
flush(6)
call IO_read_jobBinaryFile(777,'F',&
trim(getSolverJobName()),size(F))
@ -210,12 +213,12 @@ subroutine basicPETSc_init(temperature)
read (777,rec=1) temp3333_Real
close (777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(&
F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call Utilities_constitutiveResponse(&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3)
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)]),&
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)]),&
temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C_minmaxAvg
@ -237,14 +240,13 @@ type(tSolutionState) function &
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
res,&
geomdim,&
mesh_ipCoordinates,&
mesh_NcpElems, &
mesh_deformedCoordsFFT
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
grid, &
geomSize, &
tBoundaryCondition, &
Utilities_calculateRate, &
Utilities_forwardField, &
@ -295,16 +297,16 @@ type(tSolutionState) function &
write (777,rec=1) C_lastInc
close(777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(&
F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
if ( cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)])
C = C_lastInc
else
C_lastInc = C
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(&
F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
@ -319,13 +321,13 @@ type(tSolutionState) function &
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)]))
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)])
endif
F_aim = F_aim + f_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, &
rotation_BC)),[9,res(1),res(2),res(3)])
rotation_BC)),[9,grid(1),grid(2),grid(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
@ -370,23 +372,26 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
res, &
wgt
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use DAMASK_spectral_Utilities, only: &
grid, &
wgt, &
field_real, &
Utilities_FFTforward, &
Utilities_FFTbackward, &
Utilities_fourierConvolution, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS, &
debugRotation
use IO, only : IO_intOut
Utilities_divergenceRMS
use IO, only: &
IO_intOut
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: &
PetscScalar, dimension(3,3,grid(1),grid(2),grid(3)) :: &
x_scal, &
f_scal
PetscInt :: &
@ -408,7 +413,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ Iteration ', itmin, '≤',reportIter, '≤', itmax
if (debugRotation) &
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab)=', &
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
@ -433,7 +438,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(f_scal,[res(1),res(2),res(3),3,3],&
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()
err_div = Utilities_divergenceRMS()
@ -442,7 +447,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2])
f_scal = reshape(field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2])
end subroutine BasicPETSc_formResidual

295
code/DAMASK_spectral_utilities.f90
View File

@ -36,14 +36,21 @@ module DAMASK_spectral_utilities
#include <finclude/petscsys.h>
#endif
logical, public :: cutBack =.false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
!--------------------------------------------------------------------------------------------------
! grid related information information
integer(pInt), public, dimension(3) :: grid !< grid points as specified in geometry file
real(pReal), public :: wgt !< weighting factor 1/Nelems
real(pReal), public, dimension(3) :: geomSize !< size of geometry as specified in geometry file
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
integer(pInt), public :: grid1Red !< grid(1)/2
real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier
complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates
real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator
real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness
real(pReal), private, dimension(3) :: scaledGeomSize !< scaled geometry size for calculation of divergence (Basic, Basic PETSc)
!--------------------------------------------------------------------------------------------------
! debug fftw
@ -66,10 +73,9 @@ module DAMASK_spectral_utilities
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, public :: &
logical, private :: &
debugGeneral, & !< general debugging of spectral solver
debugDivergence, & !< debugging of divergence calculation (comparison to function used for post processing)
debugRestart, & !< debbuging of restart features
debugFFTW, & !< doing additional FFT on scalar field and compare to results of strided 3D FFT
debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
@ -118,22 +124,25 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine utilities_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use DAMASK_interface, only: &
geometryFile
use IO, only: &
IO_error, &
IO_warning, &
IO_timeStamp
IO_timeStamp, &
IO_open_file
use numerics, only: &
DAMASK_NumThreadsInt, &
fftw_planner_flag, &
fftw_timelimit, &
memory_efficient, &
petsc_options
petsc_options, &
divergence_correction
use debug, only: &
debug_level, &
debug_spectral, &
debug_levelBasic, &
debug_spectralDivergence, &
debug_spectralRestart, &
debug_spectralFFTW, &
debug_spectralPETSc, &
debug_spectralRotation
@ -142,11 +151,10 @@ subroutine utilities_init()
debug_spectralPETSc, &
PETScDebug
#endif
use mesh, only: &
res, &
res1_red, &
scaledDim
use math ! must use the whole module for use of FFTW
use mesh, only: &
mesh_spectral_getSize, &
mesh_spectral_getGrid
implicit none
#ifdef PETSc
@ -157,6 +165,7 @@ subroutine utilities_init()
PetscErrorCode :: ierr
#endif
integer(pInt) :: i, j, k
integer(pInt), parameter :: fileUnit = 228_pInt
integer(pInt), dimension(3) :: k_s
type(C_PTR) :: &
tensorField, & !< field cotaining data for FFTW in real and fourier space (in place)
@ -172,7 +181,6 @@ subroutine utilities_init()
! set debugging parameters
debugGeneral = iand(debug_level(debug_spectral),debug_levelBasic) /= 0
debugDivergence = iand(debug_level(debug_spectral),debug_spectralDivergence) /= 0
debugRestart = iand(debug_level(debug_spectral),debug_spectralRestart) /= 0
debugFFTW = iand(debug_level(debug_spectral),debug_spectralFFTW) /= 0
debugRotation = iand(debug_level(debug_spectral),debug_spectralRotation) /= 0
debugPETSc = iand(debug_level(debug_spectral),debug_spectralPETSc) /= 0
@ -187,12 +195,41 @@ subroutine utilities_init()
call IO_warning(41_pInt, ext_msg='debug PETSc')
#endif
call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file...
grid = mesh_spectral_getGrid(fileUnit)
grid1Red = grid(1)/2_pInt + 1_pInt
wgt = 1.0/real(product(grid),pReal)
geomSize = mesh_spectral_getSize(fileUnit)
close(fileUnit)
write(6,'(a,3(i12 ))') ' grid a b c: ', grid
write(6,'(a,3(f12.5))') ' size x y z: ', geomSize
!--------------------------------------------------------------------------------------------------
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso-
! lution-independent divergence
if (divergence_correction == 1_pInt) then
do j = 1_pInt, 3_pInt
if (j /= minloc(geomSize,1) .and. j /= maxloc(geomSize,1)) &
scaledGeomSize = geomSize/geomSize(j)
enddo
elseif (divergence_correction == 2_pInt) then
do j = 1_pInt, 3_pInt
if (j /= minloc(geomSize/grid,1) .and. j /= maxloc(geomSize/grid,1)) &
scaledGeomSize = geomSize/geomSize(j)*grid(j)
enddo
else
scaledGeomSize = geomSize
endif
!--------------------------------------------------------------------------------------------------
! allocation
allocate (xi(3,res1_red,res(2),res(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension
tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
call c_f_pointer(tensorField, field_real, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, field_fourier,[res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField
allocate (xi(3,grid1Red,grid(2),grid(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension
tensorField = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
call c_f_pointer(tensorField, field_real, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3,3])! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, field_fourier,[grid1Red, grid(2),grid(3),3,3])! place a pointer for a complex representation on tensorField
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
@ -206,41 +243,41 @@ subroutine utilities_init()
!--------------------------------------------------------------------------------------------------
! creating plans for the convolution
planForth = fftw_plan_many_dft_r2c(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! input data, physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
field_fourier,[res(3),res(2) ,res1_red], & ! output data, physical length in each dimension in reversed order
1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision
planForth = fftw_plan_many_dft_r2c(3,[grid(3),grid(2) ,grid(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! input data, physical length in each dimension in reversed order
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
field_fourier,[grid(3),grid(2) ,grid1Red], & ! output data, physical length in each dimension in reversed order
1, grid(3)*grid(2)* grid1Red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision
planBack = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_fourier,[res(3),res(2) ,res1_red], & ! input data, physical length in each dimension in reversed order
1, res(3)*res(2)* res1_red, & ! striding, product of physical length in the 3 dimensions
field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! output data, physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
planBack = fftw_plan_many_dft_c2r(3,[grid(3),grid(2) ,grid(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_fourier,[grid(3),grid(2) ,grid1Red], & ! input data, physical length in each dimension in reversed order
1, grid(3)*grid(2)* grid1Red, & ! striding, product of physical length in the 3 dimensions
field_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], & ! output data, physical length in each dimension in reversed order
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
fftw_planner_flag) ! planner precision
!--------------------------------------------------------------------------------------------------
! depending on debug options, allocate more memory and create additional plans
if (debugDivergence) then
div = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T))
call c_f_pointer(div,divReal, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3])
call c_f_pointer(div,divFourier,[res1_red, res(2),res(3),3])
planDiv = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)],3,&
divFourier,[res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,&
divReal,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], &
1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), &
div = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*3_pInt,C_SIZE_T))
call c_f_pointer(div,divReal, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3])
call c_f_pointer(div,divFourier,[grid1Red, grid(2),grid(3),3])
planDiv = fftw_plan_many_dft_c2r(3,[grid(3),grid(2) ,grid(1)],3,&
divFourier,[grid(3),grid(2) ,grid1Red],&
1, grid(3)*grid(2)* grid1Red,&
divReal,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)], &
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), &
fftw_planner_flag)
endif
if (debugFFTW) then
scalarField_realC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for real representation (no in place transform)
scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for fourier representation (no in place transform)
call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) ! place a pointer for a real representation
call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) ! place a pointer for a fourier representation
planDebugForth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style)
scalarField_realC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for real representation (no in place transform)
scalarField_fourierC = fftw_alloc_complex(int(product(grid),C_SIZE_T)) ! allocate data for fourier representation (no in place transform)
call c_f_pointer(scalarField_realC, scalarField_real, grid) ! place a pointer for a real representation
call c_f_pointer(scalarField_fourierC, scalarField_fourier, grid) ! place a pointer for a fourier representation
planDebugForth = fftw_plan_dft_3d(grid(3),grid(2),grid(1),& ! reversed order (C style)
scalarField_real,scalarField_fourier,-1,fftw_planner_flag) ! input, output, forward FFT(-1), planner precision
planDebugBack = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style)
planDebugBack = fftw_plan_dft_3d(grid(3),grid(2),grid(1),& ! reversed order (C style)
scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision
endif
@ -249,21 +286,21 @@ subroutine utilities_init()
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
do k = 1_pInt, res(3)
do k = 1_pInt, grid(3)
k_s(3) = k - 1_pInt
if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1_pInt, res(2)
if(k > grid(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - grid(3) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do j = 1_pInt, grid(2)
k_s(2) = j - 1_pInt
if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do i = 1_pInt, res1_red
if(j > grid(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - grid(2) ! running from 0,1,...,N/2,N/2+1,-N/2,-N/2+1,...,-1
do i = 1_pInt, grid1Red
k_s(1) = i - 1_pInt ! symmetry, junst running from 0,1,...,N/2,N/2+1
xi(1:3,i,j,k) = real(k_s, pReal)/scaledDim ! if divergence_correction is set, frequencies are calculated on unit length
xi(1:3,i,j,k) = real(k_s, pReal)/scaledGeomSize ! if divergence_correction is set, frequencies are calculated on unit length
enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = 0.0_pReal)
else ! precalculation of gamma_hat field
allocate (gamma_hat(3,3,3,3,res1_red ,res(2),res(3)), source =0.0_pReal)
allocate (gamma_hat(3,3,3,3,grid1Red ,grid(2),grid(3)), source = 0.0_pReal)
endif
end subroutine utilities_init
@ -284,9 +321,6 @@ subroutine utilities_updateGamma(C,saveReference)
memory_efficient
use math, only: &
math_inv33
use mesh, only: &
res, &
res1_red
implicit none
real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness
@ -307,7 +341,7 @@ subroutine utilities_updateGamma(C,saveReference)
endif
if(.not. memory_efficient) then
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k)
@ -333,10 +367,6 @@ end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTforward() !< @ToDo make row and column between randomly between 1 and 3
use math
use mesh, only : &
scaledDim, &
res, &
res1_red
implicit none
integer(pInt) :: row, column ! if debug FFTW, compare 3D array field of row and column
@ -349,7 +379,7 @@ subroutine utilities_FFTforward()
call random_number(myRand) ! two numbers: 0 <= x < 1
row = nint(myRand(1)*2_pReal + 1_pReal,pInt)
column = nint(myRand(2)*2_pReal + 1_pReal,pInt)
scalarField_real = cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) ! store the selected component
scalarField_real = cmplx(field_real(1:grid(1),1:grid(2),1:grid(3),row,column),0.0_pReal,pReal) ! store the selected component
endif
!--------------------------------------------------------------------------------------------------
@ -360,34 +390,34 @@ subroutine utilities_FFTforward()
! comparing 1 and 3x3 FT results
if (debugFFTW) then
call fftw_execute_dft(planDebugForth,scalarField_real,scalarField_fourier)
where(abs(scalarField_fourier(1:res1_red,1:res(2),1:res(3))) > tiny(1.0_pReal)) ! avoid division by zero
scalarField_fourier(1:res1_red,1:res(2),1:res(3)) = &
(scalarField_fourier(1:res1_red,1:res(2),1:res(3))-&
field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3))
where(abs(scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3))) > tiny(1.0_pReal)) ! avoid division by zero
scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)) = &
(scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3))-&
field_fourier(1:grid1Red,1:grid(2),1:grid(3),row,column))/&
scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3))
else where
scalarField_real = cmplx(0.0,0.0,pReal)
end where
write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..'
write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately
maxval(real (scalarField_fourier(1:res1_red,1:res(2),1:res(3)))),&
maxval(aimag(scalarField_fourier(1:res1_red,1:res(2),1:res(3))))
maxval(real (scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)))),&
maxval(aimag(scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3))))
flush(6)
endif
!--------------------------------------------------------------------------------------------------
! removing highest frequencies
Nyquist(2,1:2) = [res(2)/2_pInt + 1_pInt, res(2)/2_pInt + 1_pInt + mod(res(2),2_pInt)]
Nyquist(3,1:2) = [res(3)/2_pInt + 1_pInt, res(3)/2_pInt + 1_pInt + mod(res(3),2_pInt)]
Nyquist(2,1:2) = [grid(2)/2_pInt + 1_pInt, grid(2)/2_pInt + 1_pInt + mod(grid(2),2_pInt)]
Nyquist(3,1:2) = [grid(3)/2_pInt + 1_pInt, grid(3)/2_pInt + 1_pInt + mod(grid(3),2_pInt)]
if(res(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (res1_red, 1:res(2), 1:res(3), 1:3,1:3) &
if(grid(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (grid1Red, 1:grid(2), 1:grid(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (1:res1_red,Nyquist(2,1):Nyquist(2,2),1:res(3), 1:3,1:3) &
if(grid(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (1:grid1Red,Nyquist(2,1):Nyquist(2,2),1:grid(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (1:res1_red,1:res(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) &
if(grid(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (1:grid1Red,1:grid(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
end subroutine utilities_FFTforward
@ -401,10 +431,6 @@ end subroutine utilities_FFTforward
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTbackward()
use math !< must use the whole module for use of FFTW
use mesh, only: &
wgt, &
res, &
res1_red
implicit none
integer(pInt) :: row, column !< if debug FFTW, compare 3D array field of row and column
@ -417,18 +443,18 @@ subroutine utilities_FFTbackward()
call random_number(myRand) ! two numbers: 0 <= x < 1
row = nint(myRand(1)*2_pReal + 1_pReal,pInt)
column = nint(myRand(2)*2_pReal + 1_pReal,pInt)
scalarField_fourier(1:res1_red,1:res(2),1:res(3)) &
= field_fourier(1:res1_red,1:res(2),1:res(3),row,column)
do i = 0_pInt, res(1)/2_pInt-2_pInt + mod(res(1),2_pInt)
scalarField_fourier(1:grid1Red,1:grid(2),1:grid(3)) &
= field_fourier(1:grid1Red,1:grid(2),1:grid(3),row,column)
do i = 0_pInt, grid(1)/2_pInt-2_pInt + mod(grid(1),2_pInt)
m = 1_pInt
do k = 1_pInt, res(3)
do k = 1_pInt, grid(3)
n = 1_pInt
do j = 1_pInt, res(2)
scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m))
if(n == 1_pInt) n = res(2) + 1_pInt
do j = 1_pInt, grid(2)
scalarField_fourier(grid(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m))
if(n == 1_pInt) n = grid(2) + 1_pInt
n = n-1_pInt
enddo
if(m == 1_pInt) m = res(3) + 1_pInt
if(m == 1_pInt) m = grid(3) + 1_pInt
m = m -1_pInt
enddo; enddo
endif
@ -443,7 +469,7 @@ subroutine utilities_FFTbackward()
call fftw_execute_dft(planDebugBack,scalarField_fourier,scalarField_real)
where(abs(real(scalarField_real,pReal)) > tiny(1.0_pReal)) ! avoid division by zero
scalarField_real = (scalarField_real &
- cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column), 0.0, pReal))/ &
- cmplx(field_real(1:grid(1),1:grid(2),1:grid(3),row,column), 0.0, pReal))/ &
scalarField_real
else where
scalarField_real = cmplx(0.0,0.0,pReal)
@ -466,10 +492,6 @@ subroutine utilities_fourierConvolution(fieldAim)
memory_efficient
use math, only: &
math_inv33
use mesh, only: &
mesh_NcpElems, &
res, &
res1_red
implicit none
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
@ -486,7 +508,7 @@ subroutine utilities_fourierConvolution(fieldAim)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat
do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red
if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k)
@ -502,13 +524,13 @@ subroutine utilities_fourierConvolution(fieldAim)
endif
enddo; enddo; enddo
else ! use precalculated gamma-operator
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3, i,j,k) * field_fourier(i,j,k,1:3,1:3))
field_fourier(i,j,k, 1:3,1:3) = temp33_Complex
enddo; enddo; enddo
endif
field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(mesh_NcpElems,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(product(grid),pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
end subroutine utilities_fourierConvolution
@ -518,10 +540,6 @@ end subroutine utilities_fourierConvolution
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_divergenceRMS()
use math !< must use the whole module for use of FFTW
use mesh, only: &
wgt, &
res, &
res1_red
implicit none
integer(pInt) :: i, j, k
@ -537,32 +555,32 @@ real(pReal) function utilities_divergenceRMS()
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2)
do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2)
do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
enddo
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if res(1) /= 1)
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
+ sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)
+ sum( real(math_mul33x3_complex(field_fourier(grid1Red,j,k,1:3,1:3),&
xi(1:3,grid1Red,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(field_fourier(grid1Red,j,k,1:3,1:3),&
xi(1:3,grid1Red,j,k))*TWOPIIMG)**2.0_pReal)
enddo; enddo
if(res(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of res(1) == 1
if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
!--------------------------------------------------------------------------------------------------
! calculate additional divergence criteria and report
if (debugDivergence) then ! calculate divergence again
err_div_max = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier
xi(1:3,i,j,k))*TWOPIIMG
err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal))
@ -590,10 +608,6 @@ end function utilities_divergenceRMS
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_curlRMS()
use math !< must use the whole module for use of FFTW
use mesh, only: &
res, &
res1_red, &
wgt
implicit none
integer(pInt) :: i, j, k, l
@ -605,8 +619,8 @@ real(pReal) function utilities_curlRMS()
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2);
do i = 2_pInt, res1_red - 1_pInt
do k = 1_pInt, grid(3); do j = 1_pInt, grid(2);
do i = 2_pInt, grid1Red - 1_pInt
do l = 1_pInt, 3_pInt
curl_fourier(l,1) = (field_fourier(i,j,k,l,3)*xi(2,i,j,k)&
- field_fourier(i,j,k,l,2)*xi(3,i,j,k))*TWOPIIMG
@ -629,12 +643,12 @@ real(pReal) function utilities_curlRMS()
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)
do l = 1_pInt, 3_pInt
curl_fourier = (field_fourier(res1_red,j,k,l,3)*xi(2,res1_red,j,k)&
- field_fourier(res1_red,j,k,l,2)*xi(3,res1_red,j,k))*TWOPIIMG
curl_fourier = (-field_fourier(res1_red,j,k,l,3)*xi(1,res1_red,j,k)&
+field_fourier(res1_red,j,k,l,1)*xi(3,res1_red,j,k) )*TWOPIIMG
curl_fourier = ( field_fourier(res1_red,j,k,l,2)*xi(1,res1_red,j,k)&
-field_fourier(res1_red,j,k,l,1)*xi(2,res1_red,j,k) )*TWOPIIMG
curl_fourier = ( field_fourier(grid1Red,j,k,l,3)*xi(2,grid1Red,j,k)&
-field_fourier(grid1Red,j,k,l,2)*xi(3,grid1Red,j,k))*TWOPIIMG
curl_fourier = (-field_fourier(grid1Red,j,k,l,3)*xi(1,grid1Red,j,k)&
+field_fourier(grid1Red,j,k,l,1)*xi(3,grid1Red,j,k))*TWOPIIMG
curl_fourier = ( field_fourier(grid1Red,j,k,l,2)*xi(1,grid1Red,j,k)&
-field_fourier(grid1Red,j,k,l,1)*xi(2,grid1Red,j,k))*TWOPIIMG
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)
@ -758,10 +772,6 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
math_det33
use FEsolving, only: &
restartWrite
use mesh, only: &
res, &
wgt, &
mesh_NcpElems
use CPFEM, only: &
CPFEM_general, &
CPFEM_COLLECT, &
@ -778,16 +788,16 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
implicit none
real(pReal), intent(inout) :: temperature !< temperature (no field)
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: &
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid(3)) :: &
F_lastInc, & !< target deformation gradient
F !< previous deformation gradient
real(pReal), intent(in) :: timeinc !< loading time
logical, intent(in) :: forwardData !< age results
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,res(1),res(2),res(3)) :: P !< PK stress
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid(3)) :: P !< PK stress
integer(pInt) :: &
calcMode, & !< CPFEM mode for calculation
@ -812,8 +822,8 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
call CPFEM_general(collectMode,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration
temperature,timeinc,1_pInt,1_pInt)
materialpoint_F0 = reshape(F_lastInc, [3,3,1,mesh_NcpElems])
materialpoint_F = reshape(F, [3,3,1,mesh_NcpElems])
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)])
materialpoint_F = reshape(F, [3,3,1,product(grid)])
materialpoint_Temperature = temperature
call debug_reset()
@ -823,7 +833,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
if(debugGeneral) then
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do j = 1_pInt, mesh_NcpElems
do j = 1_pInt, product(grid)
defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
@ -840,7 +850,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
max_dPdF_norm = 0.0_pReal
min_dPdF = huge(1.0_pReal)
min_dPdF_norm = huge(1.0_pReal)
do k = 1_pInt, mesh_NcpElems
do k = 1_pInt, product(grid)
if (max_dPdF_norm < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then
max_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)
max_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)
@ -851,7 +861,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
endif
end do
P = reshape(materialpoint_P, [3,3,res(1),res(2),res(3)])
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid(3)])
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
@ -875,24 +885,22 @@ end subroutine utilities_constitutiveResponse
!> @brief calculates forward rate, either guessing or just add delta/timeinc
!--------------------------------------------------------------------------------------------------
pure function utilities_calculateRate(avRate,timeinc_old,guess,field_lastInc,field)
use mesh, only: &
res
implicit none
real(pReal), intent(in), dimension(3,3) :: avRate !< homogeneous addon
real(pReal), intent(in) :: &
timeinc_old !< timeinc of last step
logical, intent(in) :: &
guess !< guess along former trajectory
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: &
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid(3)) :: &
field_lastInc, & !< data of previous step
field !< data of current step
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_calculateRate
real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: utilities_calculateRate
if(guess) then
utilities_calculateRate = (field-field_lastInc) / timeinc_old
else
utilities_calculateRate = spread(spread(spread(avRate,3,res(1)),4,res(2)),5,res(3))
utilities_calculateRate = spread(spread(spread(avRate,3,grid(1)),4,grid(2)),5,grid(3))
endif
end function utilities_calculateRate
@ -903,26 +911,23 @@ end function utilities_calculateRate
!> ensures that the average matches the aim
!--------------------------------------------------------------------------------------------------
pure function utilities_forwardField(timeinc,field_lastInc,rate,aim)
use mesh, only: &
res, &
wgt
implicit none
real(pReal), intent(in) :: &
timeinc !< timeinc of current step
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: &
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid(3)) :: &
field_lastInc, & !< initial field
rate !< rate by which to forward
real(pReal), intent(in), optional, dimension(3,3) :: &
aim !< average field value aim
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_forwardField
real(pReal), dimension(3,3,grid(1),grid(2),grid(3)) :: utilities_forwardField
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
utilities_forwardField = field_lastInc + rate*timeinc
if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim
utilities_forwardField = utilities_forwardField - &
spread(spread(spread(fieldDiff,3,res(1)),4,res(2)),5,res(3))
spread(spread(spread(fieldDiff,3,grid(1)),4,grid(2)),5,grid(3))
endif
end function utilities_forwardField
@ -936,8 +941,6 @@ real(pReal) function utilities_getFilter(k)
IO_error
use numerics, only: &
myfilter
use mesh, only: &
res
use math, only: &
PI
@ -949,9 +952,9 @@ real(pReal) function utilities_getFilter(k)
select case (myfilter)
case ('none') !< default is already nothing (1.0_pReal)
case ('cosine') !< cosine curve with 1 for avg and zero for highest freq
utilities_getFilter = (1.0_pReal + cos(PI*k(3)/res(3))) &
*(1.0_pReal + cos(PI*k(2)/res(2))) &
*(1.0_pReal + cos(PI*k(1)/res(1)))/8.0_pReal
utilities_getFilter = (1.0_pReal + cos(PI*k(3)/grid(3))) &
*(1.0_pReal + cos(PI*k(2)/grid(2))) &
*(1.0_pReal + cos(PI*k(1)/grid(1)))/8.0_pReal
case default
call IO_error(error_ID = 892_pInt, ext_msg = trim(myfilter))
end select

3
code/lattice.f90
View File

@ -1018,7 +1018,8 @@ pure function lattice_symmetrizeC66(structName,C66)
! TwinTwinInteraction
!--------------------------------------------------------------------------------------------------
function lattice_configNchunks(struct)
use prec, only: pReal,pInt
use prec, only: &
pInt
implicit none
integer(pInt), dimension(6) :: lattice_configNchunks

10
code/math.f90
View File

@ -727,14 +727,14 @@ pure function math_exp33(A,n)
order = 5
if (present(n)) order = n
B = math_identity2nd(3) ! init
invfac = 1.0_pReal ! 0!
math_exp33 = B ! A^0 = eye2
B = math_identity2nd(3) ! init
invfac = 1.0_pReal ! 0!
math_exp33 = B ! A^0 = eye2
do i = 1_pInt,n
invfac = invfac/real(i) ! invfac = 1/i!
invfac = invfac/real(i) ! invfac = 1/i!
B = math_mul33x33(B,A)
math_exp33 = math_exp33 + invfac*B ! exp = SUM (A^i)/i!
math_exp33 = math_exp33 + invfac*B ! exp = SUM (A^i)/i!
enddo
end function math_exp33

413
code/mesh.f90
View File

@ -121,16 +121,6 @@ module mesh
#ifdef Spectral
include 'fftw3.f03'
real(pReal), dimension(3), public, protected :: &
geomdim, & !< physical dimension of volume element per direction
scaledDim !< scaled dimension of volume element, depending on selected divergence calculation
integer(pInt), dimension(3), public, protected :: &
res !< resolution, e.g. number of Fourier points in each direction
real(pReal), public, protected :: &
wgt
integer(pInt), public, protected :: &
res1_red, &
homog
#endif
! These definitions should actually reside in the FE-solver specific part (different for MARC/ABAQUS)
@ -428,6 +418,8 @@ module mesh
mesh_get_nodeAtIP
#ifdef Spectral
public :: &
mesh_spectral_getGrid, &
mesh_spectral_getSize, &
mesh_regrid, &
mesh_nodesAroundCentres, &
mesh_deformedCoordsFFT, &
@ -438,13 +430,9 @@ module mesh
private :: &
#ifdef Spectral
mesh_spectral_getGrid, &
mesh_spectral_getSize, &
mesh_spectral_getHomogenization, &
mesh_spectral_count_nodesAndElements, &
mesh_spectral_count_cpElements, &
mesh_spectral_map_elements, &
mesh_spectral_map_nodes, &
mesh_spectral_count, &
mesh_spectral_mapNodesAndElems, &
mesh_spectral_count_cpSizes, &
mesh_spectral_build_nodes, &
mesh_spectral_build_elements, &
@ -554,51 +542,23 @@ subroutine mesh_init(ip,el)
if (allocated(FE_ipNeighbor)) deallocate(FE_ipNeighbor)
if (allocated(FE_cellnodeParentnodeWeights)) deallocate(FE_cellnodeParentnodeWeights)
if (allocated(FE_subNodeOnIPFace)) deallocate(FE_subNodeOnIPFace)
call mesh_build_FEdata ! get properties of the different types of elements
mesh_unitlength = numerics_unitlength ! set physical extent of a length unit in mesh
#ifdef Spectral
call IO_open_file(fileUnit,geometryFile) ! parse info from geometry file...
res = mesh_spectral_getGrid(fileUnit)
res1_red = res(1)/2_pInt + 1_pInt
wgt = 1.0/real(product(res),pReal)
geomdim = mesh_spectral_getSize(fileUnit)
homog = mesh_spectral_getHomogenization(fileUnit)
!--------------------------------------------------------------------------------------------------
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and reso-
! lution-independent divergence
if (divergence_correction == 1_pInt) then
do j = 1_pInt, 3_pInt
if (j /= minloc(geomdim,1) .and. j /= maxloc(geomdim,1)) scaledDim = geomdim/geomdim(j)
enddo
elseif (divergence_correction == 2_pInt) then
do j = 1_pInt, 3_pInt
if (j /= minloc(geomdim/res,1) .and. j /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(j)*res(j)
enddo
else
scaledDim = geomdim
endif
write(6,'(a,3(i12 ))') ' grid a b c: ', res
write(6,'(a,3(f12.5))') ' size x y z: ', geomdim
write(6,'(a,i5,/)') ' homogenization: ', homog
call mesh_spectral_count_nodesAndElements
call mesh_spectral_count_cpElements
call mesh_spectral_map_elements
call mesh_spectral_map_nodes
call mesh_spectral_count(fileUnit)
call mesh_spectral_mapNodesAndElems
call mesh_spectral_count_cpSizes
call mesh_spectral_build_nodes
call mesh_spectral_build_nodes(fileUnit)
call mesh_spectral_build_elements(fileUnit)
call mesh_get_damaskOptions(fileUnit)
close (fileUnit)
call mesh_build_cellconnectivity
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes)
call mesh_build_ipCoordinates
call mesh_build_ipVolumes
call mesh_build_ipAreas
call mesh_spectral_build_ipNeighborhood
call mesh_spectral_build_ipNeighborhood(fileUnit)
#endif
#ifdef Marc4DAMASK
call IO_open_inputFile(fileUnit,modelName) ! parse info from input file...
@ -613,7 +573,6 @@ subroutine mesh_init(ip,el)
call mesh_marc_count_cpSizes(fileunit)
call mesh_marc_build_elements(fileUnit)
call mesh_get_damaskOptions(fileUnit)
close (fileUnit)
call mesh_build_cellconnectivity
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes)
call mesh_build_ipCoordinates
@ -638,7 +597,6 @@ subroutine mesh_init(ip,el)
call mesh_abaqus_count_cpSizes(fileunit)
call mesh_abaqus_build_elements(fileUnit)
call mesh_get_damaskOptions(fileUnit)
close (fileUnit)
call mesh_build_cellconnectivity
mesh_cellnode = mesh_build_cellnodes(mesh_node,mesh_Ncellnodes)
call mesh_build_ipCoordinates
@ -649,6 +607,7 @@ subroutine mesh_init(ip,el)
call mesh_build_ipNeighborhood
#endif
close (fileUnit)
call mesh_tell_statistics
call mesh_write_meshfile
call mesh_write_cellGeom
@ -950,25 +909,24 @@ end subroutine mesh_build_ipCoordinates
!--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(ip,el)
implicit none
implicit none
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
integer(pInt), intent(in) :: el, & !< element number
ip !< integration point number
real(pReal), dimension(3) :: mesh_cellCenterCoordinates !< x,y,z coordinates of the cell center of the requested IP cell
integer(pInt) :: t,g,c,n
integer(pInt) :: t,g,c,n
t = mesh_element(2_pInt,el) ! get element type
g = FE_geomtype(t) ! get geometry type
c = FE_celltype(g) ! get cell type
mesh_cellCenterCoordinates = 0.0_pReal
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
enddo
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c)
t = mesh_element(2_pInt,el) ! get element type
g = FE_geomtype(t) ! get geometry type
c = FE_celltype(g) ! get cell type
mesh_cellCenterCoordinates = 0.0_pReal
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
enddo
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c)
endfunction mesh_cellCenterCoordinates
end function mesh_cellCenterCoordinates
#ifdef Spectral
@ -1189,61 +1147,39 @@ end function mesh_spectral_getHomogenization
!--------------------------------------------------------------------------------------------------
!> @brief Count overall number of nodes and elements in mesh and stores them in
!! 'mesh_Nelems' and 'mesh_Nnodes'
!! 'mesh_Nelems', 'mesh_Nnodes' and 'mesh_NcpElems'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_count_nodesAndElements()
subroutine mesh_spectral_count(myUnit)
implicit none
mesh_Nelems = product(res)
mesh_Nnodes = product(res+1_pInt)
end subroutine mesh_spectral_count_nodesAndElements
!--------------------------------------------------------------------------------------------------
!> @brief Count overall number of CP elements in mesh and stores them in 'mesh_NcpElems'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_count_cpElements
implicit none
integer(pInt), intent(in) :: myUnit
integer(pInt), dimension(3) :: grid
grid = mesh_spectral_getGrid(myUnit)
mesh_Nelems = product(grid)
mesh_NcpElems = mesh_Nelems
end subroutine mesh_spectral_count_cpElements
mesh_Nnodes = product(grid+1_pInt)
end subroutine mesh_spectral_count
!--------------------------------------------------------------------------------------------------
!> @brief Maps elements from FE ID to internal (consecutive) representation.
!! Allocates global array 'mesh_mapFEtoCPelem'
!> @brief fake map node from FE ID to internal (consecutive) representation for node and element
!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_map_elements
subroutine mesh_spectral_mapNodesAndElems
use math, only: &
math_range
implicit none
integer(pInt) :: i
allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems)) ; mesh_mapFEtoCPelem = 0_pInt
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt)
allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt)
forall (i = 1_pInt:mesh_NcpElems) &
mesh_mapFEtoCPelem(1:2,i) = i
mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2)
mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2)
end subroutine mesh_spectral_map_elements
!--------------------------------------------------------------------------------------------------
!> @brief Maps node from FE ID to internal (consecutive) representation.
!! Allocates global array 'mesh_mapFEtoCPnode'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_map_nodes
implicit none
integer(pInt) :: i
allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes)) ; mesh_mapFEtoCPnode = 0_pInt
forall (i = 1_pInt:mesh_Nnodes) &
mesh_mapFEtoCPnode(1:2,i) = i
end subroutine mesh_spectral_map_nodes
end subroutine mesh_spectral_mapNodesAndElems
!--------------------------------------------------------------------------------------------------
@ -1272,24 +1208,30 @@ end subroutine mesh_spectral_count_cpSizes
!> @brief Store x,y,z coordinates of all nodes in mesh.
!! Allocates global arrays 'mesh_node0' and 'mesh_node'
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_nodes()
subroutine mesh_spectral_build_nodes(myUnit)
implicit none
integer(pInt) :: n
integer(pInt), intent(in) :: myUnit
integer(pInt) :: n
integer(pInt), dimension(3) :: grid
real(pReal), dimension(3) :: geomSize
allocate ( mesh_node0 (3,mesh_Nnodes) ); mesh_node0 = 0.0_pReal
allocate ( mesh_node (3,mesh_Nnodes) ); mesh_node = 0.0_pReal
allocate (mesh_node0 (3,mesh_Nnodes), source = 0.0_pReal)
allocate (mesh_node (3,mesh_Nnodes), source = 0.0_pReal)
grid = mesh_spectral_getGrid(myUnit)
geomSize = mesh_spectral_getSize(myUnit)
forall (n = 0_pInt:mesh_Nnodes-1_pInt)
mesh_node0(1,n+1_pInt) = mesh_unitlength * &
geomdim(1) * real(mod(n,(res(1)+1_pInt) ),pReal) &
/ real(res(1),pReal)
geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) &
/ real(grid(1),pReal)
mesh_node0(2,n+1_pInt) = mesh_unitlength * &
geomdim(2) * real(mod(n/(res(1)+1_pInt),(res(2)+1_pInt)),pReal) &
/ real(res(2),pReal)
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 * &
geomdim(3) * real(mod(n/(res(1)+1_pInt)/(res(2)+1_pInt),(res(3)+1_pInt)),pReal) &
/ real(res(3),pReal)
geomSize(3)*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid(3)+1_pInt)),pReal) &
/ real(grid(3),pReal)
end forall
mesh_node = mesh_node0
@ -1314,15 +1256,29 @@ subroutine mesh_spectral_build_elements(myUnit)
IO_countContinuousIntValues
implicit none
integer(pInt), intent(in) :: myUnit
integer(pInt), dimension (1_pInt+7_pInt*2_pInt) :: myPos
integer(pInt) :: e, i, headerLength = 0_pInt, maxIntCount
integer(pInt), dimension(:), allocatable :: microstructures
integer(pInt), dimension(1,1) :: dummySet = 0_pInt
character(len=65536) :: line,keyword
character(len=64), dimension(1) :: dummyName = ''
integer(pInt), intent(in) :: &
myUnit
integer(pInt), dimension(1_pInt+7_pInt*2_pInt) :: &
myPos
integer(pInt) :: &
e, i, &
headerLength = 0_pInt, &
maxIntCount, &
homog
integer(pInt), dimension(:), allocatable :: &
microstructures
integer(pInt), dimension(3) :: &
grid
integer(pInt), dimension(1,1) :: &
dummySet = 0_pInt
character(len=65536) :: &
line, &
keyword
character(len=64), dimension(1) :: &
dummyName = ''
grid = mesh_spectral_getGrid(myUnit)
homog = mesh_spectral_getHomogenization(myUnit)
call IO_checkAndRewind(myUnit)
read(myUnit,'(a65536)') line
@ -1364,15 +1320,15 @@ subroutine mesh_spectral_build_elements(myUnit)
mesh_element( 2,e) = FE_mapElemtype('C3D8R') ! 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)/res(1) + &
((e-1_pInt)/(res(1)*res(2)))*(res(1)+1_pInt) ! base node
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( 6,e) = mesh_element(5,e) + 1_pInt
mesh_element( 7,e) = mesh_element(5,e) + res(1) + 2_pInt
mesh_element( 8,e) = mesh_element(5,e) + res(1) + 1_pInt
mesh_element( 9,e) = mesh_element(5,e) +(res(1) + 1_pInt) * (res(2) + 1_pInt) ! second floor base node
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(10,e) = mesh_element(9,e) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + res(1) + 2_pInt
mesh_element(12,e) = mesh_element(9,e) + res(1) + 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_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
@ -1388,56 +1344,59 @@ end subroutine mesh_spectral_build_elements
!> @brief build neighborhood relations for spectral
!> @details assign globals: mesh_ipNeighborhood
!--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_ipNeighborhood
subroutine mesh_spectral_build_ipNeighborhood(myUnit)
implicit none
integer(pInt) x,y,z, &
e
allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems))
mesh_ipNeighborhood = 0_pInt
e = 0_pInt
do z = 0_pInt,res(3)-1_pInt
do y = 0_pInt,res(2)-1_pInt
do x = 0_pInt,res(1)-1_pInt
e = e + 1_pInt
mesh_ipNeighborhood(1,1,1,e) = z * res(1) * res(2) &
+ y * res(1) &
+ modulo(x+1_pInt,res(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,2,1,e) = z * res(1) * res(2) &
+ y * res(1) &
+ modulo(x-1_pInt,res(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,3,1,e) = z * res(1) * res(2) &
+ modulo(y+1_pInt,res(2)) * res(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,4,1,e) = z * res(1) * res(2) &
+ modulo(y-1_pInt,res(2)) * res(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,res(3)) * res(1) * res(2) &
+ y * res(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,res(3)) * res(1) * res(2) &
+ y * res(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt
mesh_ipNeighborhood(3,1,1,e) = 2_pInt
mesh_ipNeighborhood(3,2,1,e) = 1_pInt
mesh_ipNeighborhood(3,3,1,e) = 4_pInt
mesh_ipNeighborhood(3,4,1,e) = 3_pInt
mesh_ipNeighborhood(3,5,1,e) = 6_pInt
mesh_ipNeighborhood(3,6,1,e) = 5_pInt
enddo
enddo
enddo
implicit none
integer(pInt), intent(in) :: &
myUnit
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(myUnit)
e = 0_pInt
do z = 0_pInt,grid(1)-1_pInt
do y = 0_pInt,grid(2)-1_pInt
do x = 0_pInt,grid(3)-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)) &
+ 1_pInt
mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x-1_pInt,grid(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) &
+ modulo(y+1_pInt,grid(2)) * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) &
+ modulo(y-1_pInt,grid(2)) * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid(3)) * grid(1) * grid(2) &
+ y * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt
mesh_ipNeighborhood(3,1,1,e) = 2_pInt
mesh_ipNeighborhood(3,2,1,e) = 1_pInt
mesh_ipNeighborhood(3,3,1,e) = 4_pInt
mesh_ipNeighborhood(3,4,1,e) = 3_pInt
mesh_ipNeighborhood(3,5,1,e) = 6_pInt
mesh_ipNeighborhood(3,6,1,e) = 5_pInt
enddo
enddo
enddo
end subroutine mesh_spectral_build_ipNeighborhood
@ -1454,6 +1413,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
getSolverJobName, &
GeometryFile
use IO, only: &
IO_open_file, &
IO_read_jobBinaryFile ,&
IO_read_jobBinaryIntFile ,&
IO_write_jobBinaryFile, &
@ -1467,16 +1427,17 @@ function mesh_regrid(adaptive,resNewInput,minRes)
math_mul33x3
implicit none
character(len=1024):: formatString, N_Digits
logical, intent(in) :: adaptive ! if true, choose adaptive grid based on resNewInput, otherwise keep it constant
integer(pInt), dimension(3), optional, intent(in) :: resNewInput ! f2py cannot handle optional arguments correctly (they are always present)
integer(pInt), dimension(3), optional, intent(in) :: minRes
integer(pInt), dimension(3) :: mesh_regrid, ratio
integer(pInt), dimension(3) :: mesh_regrid, ratio, grid
integer(pInt), parameter :: myUnit = 777_pInt
integer(pInt), dimension(3,2) :: possibleResNew
integer(pInt):: maxsize, i, j, k, ielem, NpointsNew, spatialDim
integer(pInt):: maxsize, i, j, k, ielem, NpointsNew, spatialDim, Nelems
integer(pInt), dimension(3) :: resNew
integer(pInt), dimension(:), allocatable :: indices
real(pReal), dimension(3) :: geomdimNew
real(pReal) :: wgt
real(pReal), dimension(3) :: geomSizeNew, geomSize
real(pReal), dimension(3,3) :: Favg, Favg_LastInc, &
FavgNew, Favg_LastIncNew, &
deltaF, deltaF_lastInc
@ -1497,13 +1458,21 @@ function mesh_regrid(adaptive,resNewInput,minRes)
F_lastIncNew
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: &
dPdF, dPdFNew
character(len=1024):: formatString, N_Digits
integer(pInt), dimension(:,:), allocatable :: &
sizeStateHomog
integer(pInt), dimension(:,:,:), allocatable :: &
material_phase, material_phaseNew, &
sizeStateConst
call IO_open_file(myUnit,trim(geometryFile))
grid = mesh_spectral_getGrid(myUnit)
geomSize = mesh_spectral_getsize(myUnit)
close(myUnit)
Nelems = product(grid)
wgt = 1.0_pReal/real(Nelems,pReal)
write(6,'(a)') 'Regridding geometry'
if (adaptive) then
write(6,'(a)') 'adaptive resolution determination'
@ -1519,53 +1488,54 @@ function mesh_regrid(adaptive,resNewInput,minRes)
endif
endif
allocate(coordinates(3,mesh_NcpElems))
allocate(coordinates(3,Nelems))
!--------------------------------------------------------------------------------------------------
! read in deformation gradient to calculate coordinates, shape depend of selected solver
select case(myspectralsolver)
case('basic')
allocate(spectralF33(3,3,res(1),res(2),res(3)))
allocate(spectralF33(3,3,grid(1),grid(2),grid(3)))
call IO_read_jobBinaryFile(777,'F',trim(getSolverJobName()),size(spectralF33))
read (777,rec=1) spectralF33
close (777)
Favg = sum(sum(sum(spectralF33,dim=5),dim=4),dim=3) * wgt
coordinates = reshape(mesh_deformedCoordsFFT(geomdim,spectralF33),[3,mesh_NcpElems])
coordinates = reshape(mesh_deformedCoordsFFT(geomSize,spectralF33),[3,mesh_NcpElems])
case('basicpetsc','al')
allocate(spectralF9(9,res(1),res(2),res(3)))
allocate(spectralF9(9,grid(1),grid(2),grid(3)))
call IO_read_jobBinaryFile(777,'F',trim(getSolverJobName()),size(spectralF9))
read (777,rec=1) spectralF9
close (777)
Favg = reshape(sum(sum(sum(spectralF9,dim=4),dim=3),dim=2) * wgt, [3,3])
coordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(spectralF9, &
[3,3,res(1),res(2),res(3)])),[3,mesh_NcpElems])
coordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(spectralF9, &
[3,3,grid(1),grid(2),grid(3)])),[3,mesh_NcpElems])
end select
!--------------------------------------------------------------------------------------------------
! sanity check 2D/3D case
if (res(3)== 1_pInt) then
if (grid(3)== 1_pInt) then
spatialDim = 2_pInt
if (present (minRes)) then
if (minRes(1) > 0_pInt .or. minRes(2) > 0_pInt) then
if (minRes(3) /= 1_pInt .or. &
mod(minRes(1),2_pInt) /= 0_pInt .or. &
mod(minRes(2),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '2D minRes') ! as f2py has problems with present, use pyf file for initialization to -1
mod(minRes(2),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '2D minRes') ! as f2py has problems with present, use pyf file for initialization to -1
endif; endif
else
spatialDim = 3_pInt
if (present (minRes)) then
if (any(minRes > 0_pInt)) then
if (mod(minRes(1),2_pInt) /= 0_pInt.or. &
if (mod(minRes(1),2_pInt) /= 0_pInt .or. &
mod(minRes(2),2_pInt) /= 0_pInt .or. &
mod(minRes(3),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '3D minRes') ! as f2py has problems with present, use pyf file for initialization to -1
mod(minRes(3),2_pInt) /= 0_pInt) call IO_error(890_pInt, ext_msg = '3D minRes') ! as f2py has problems with present, use pyf file for initialization to -1
endif; endif
endif
!--------------------------------------------------------------------------------------------------
! Automatic detection based on current geom
geomdimNew = math_mul33x3(Favg,geomdim)
geomSizeNew = math_mul33x3(Favg,geomSize)
if (adaptive) then
ratio = floor(real(resNewInput,pReal) * (geomdimNew/geomdim), pInt)
ratio = floor(real(resNewInput,pReal) * (geomSizeNew/geomSize), pInt)
possibleResNew = 1_pInt
do i = 1_pInt, spatialDim
@ -1574,15 +1544,14 @@ function mesh_regrid(adaptive,resNewInput,minRes)
else
possibleResNew(i,1:2) = [ratio(i)-1_pInt, ratio(i) + 1_pInt]
endif
if (.not.present(minRes)) then ! calling from fortran, optional argument not given
if (.not.present(minRes)) then ! calling from fortran, optional argument not given
possibleResNew = possibleResNew
else ! optional argument is there
else ! optional argument is there
if (any(minRes<1_pInt)) then
possibleResNew = possibleResNew ! f2py calling, but without specification (or choosing invalid values), standard from pyf = -1
else ! given useful values
do k = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
possibleResNew(j,k) = max(possibleResNew(j,k), minRes(j))
enddo; enddo
possibleResNew = possibleResNew ! f2py calling, but without specification (or choosing invalid values), standard from pyf = -1
else ! given useful values
forall(k = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
possibleResNew(j,k) = max(possibleResNew(j,k), minRes(j))
endif
endif
enddo
@ -1602,11 +1571,11 @@ function mesh_regrid(adaptive,resNewInput,minRes)
endif
enddo
else
resNew = res
resNew = grid
endif
mesh_regrid = resNew
NpointsNew = resNew(1)*resNew(2)*resNew(3)
NpointsNew = product(resNew)
!--------------------------------------------------------------------------------------------------
! Calculate regular new coordinates
@ -1614,14 +1583,14 @@ function mesh_regrid(adaptive,resNewInput,minRes)
ielem = 0_pInt
do k=1_pInt,resNew(3); do j=1_pInt, resNew(2); do i=1_pInt, resNew(1)
ielem = ielem + 1_pInt
coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomdim/real(resNew,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*resNew,pReal))
coordinatesNew(1:3,ielem) = math_mul33x3(Favg, geomSize/real(resNew,pReal)*real([i,j,k],pReal) &
- geomSize/real(2_pInt*resNew,pReal))
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! Nearest neighbour search
allocate(indices(NpointsNew))
indices = math_periodicNearestNeighbor(geomdim, Favg, coordinatesNew, coordinates)
indices = math_periodicNearestNeighbor(geomSize, Favg, coordinatesNew, coordinates)
deallocate(coordinates)
!--------------------------------------------------------------------------------------------------
@ -1665,8 +1634,8 @@ function mesh_regrid(adaptive,resNewInput,minRes)
formatString = '(I'//trim(N_Digits)//'.'//trim(N_Digits)//',a)'
open(777,file=trim(getSolverWorkingDirectoryName())//trim(GeometryFile),status='REPLACE')
write(777, '(A)') '3 header'
write(777, '(3(A, I8))') 'resolution a ', resNew(1), ' b ', resNew(2), ' c ', resNew(3)
write(777, '(3(A, g17.10))') 'dimension x ', geomdim(1), ' y ', geomdim(2), ' z ', geomdim(3)
write(777, '(3(A, I8))') 'grid a ', resNew(1), ' b ', resNew(2), ' c ', resNew(3)
write(777, '(3(A, g17.10))') 'size x ', geomSize(1), ' y ', geomSize(2), ' z ', geomSize(3)
write(777, '(A)') 'homogenization 1'
do i = 1_pInt, NpointsNew
write(777,trim(formatString),advance='no') mesh_element(4,indices(i)), ' '
@ -2763,11 +2732,12 @@ end subroutine mesh_marc_map_nodes
!--------------------------------------------------------------------------------------------------
subroutine mesh_marc_build_nodes(myUnit)
use IO, only: IO_lc, &
IO_stringValue, &
IO_stringPos, &
IO_fixedIntValue, &
IO_fixedNoEFloatValue
use IO, only: &
IO_lc, &
IO_stringValue, &
IO_stringPos, &
IO_fixedIntValue, &
IO_fixedNoEFloatValue
implicit none
integer(pInt), intent(in) :: myUnit
@ -3464,7 +3434,7 @@ subroutine mesh_abaqus_build_nodes(myUnit)
myPos = IO_stringPos(line,maxNchunks)
m = mesh_FEasCP('node',IO_intValue(line,myPos,1_pInt))
do j=1_pInt, 3_pInt
mesh_node0(j,m) = mesh_unitlength * IO_floatValue(line,myPos,j+1_pInt)
mesh_node0(j,m) = numerics_unitlength * IO_floatValue(line,myPos,j+1_pInt)
enddo
enddo
endif
@ -3705,9 +3675,10 @@ use IO, only: &
endselect
endif
enddo
#endif
610 FORMAT(A300)
#endif
620 end subroutine mesh_get_damaskOptions
@ -3715,19 +3686,17 @@ use IO, only: &
!> @brief calculation of IP interface areas, allocate globals '_ipArea', and '_ipAreaNormal'
!--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipAreas
use math, only: &
math_norm3, &
math_vectorproduct
implicit none
integer(pInt) :: e,t,g,c,i,f,n,m
real(pReal), dimension (3,FE_maxNcellnodesPerCellface) :: nodePos, normals
real(pReal), dimension(3) :: normal
allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipArea = 0.0_pReal
allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)) ; mesh_ipAreaNormal = 0.0_pReal
allocate(mesh_ipArea(mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)); mesh_ipArea = 0.0_pReal
allocate(mesh_ipAreaNormal(3_pInt,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems)); mesh_ipAreaNormal = 0.0_pReal
!$OMP PARALLEL DO PRIVATE(t,g,c,nodePos,normal,normals)
do e = 1_pInt,mesh_NcpElems ! loop over cpElems