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

Merge remote-tracking branch 'remotes/origin/Grid-FEM-2' into development

This commit is contained in:
Franz Roters 2019-03-27 16:57:20 +01:00
parent df5cb3baf2 6070c0c2b1
commit 1e20e94792
19 changed files with 1810 additions and 1126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

96
.gitlab-ci.yml
View File

@ -4,8 +4,8 @@ stages:
- preprocessing
- postprocessing
- compilePETSc
- prepareSpectral
- spectral
- prepareGrid
- grid
- compileMarc
- marc
- compileAbaqus
@ -141,9 +141,9 @@ Pre_General:
- master
- release
Spectral_geometryPacking:
grid_geometryPacking:
stage: preprocessing
script: Spectral_geometryPacking/test.py
script: grid_geometryPacking/test.py
except:
- master
- release
@ -172,7 +172,7 @@ Post_General:
Post_GeometryReconstruction:
stage: postprocessing
script: Spectral_geometryReconstruction/test.py
script: spectral_geometryReconstruction/test.py
except:
- master
- release
@ -215,12 +215,12 @@ Post_OrientationConversion:
- release
###################################################################################################
Compile_Spectral_Intel:
grid_mech_compile_Intel:
stage: compilePETSc
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cp -r SpectralAll_compile SpectralAll_compile_Intel
- SpectralAll_compile_Intel/test.py
- cp -r grid_mech_compile grid_mech_compile_Intel
- grid_mech_compile_Intel/test.py
except:
- master
- release
@ -235,12 +235,12 @@ Compile_FEM_Intel:
- master
- release
Compile_Spectral_GNU:
grid_mech_compile_GNU:
stage: compilePETSc
script:
- module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
- cp -r SpectralAll_compile SpectralAll_compile_GNU
- SpectralAll_compile_GNU/test.py
- cp -r grid_mech_compile grid_mech_compile_GNU
- grid_mech_compile_GNU/test.py
except:
- master
- release
@ -257,134 +257,134 @@ Compile_FEM_GNU:
###################################################################################################
Compile_Intel_Prepare:
stage: prepareSpectral
stage: prepareGrid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd $DAMASKROOT
- make clean spectral processing
- make clean grid processing
except:
- master
- release
###################################################################################################
Thermal:
stage: spectral
stage: grid
script: Thermal/test.py
except:
- master
- release
Spectral_PackedGeometry:
stage: spectral
script: Spectral_PackedGeometry/test.py
grid_packedGeometry:
stage: grid
script: grid_packedGeometry/test.py
except:
- master
- release
Spectral_parsingArguments:
stage: spectral
script: Spectral_parsingArguments/test.py
grid_parsingArguments:
stage: grid
script: grid_parsingArguments/test.py
except:
- master
- release
StateIntegration_compareVariants:
stage: spectral
stage: grid
script: StateIntegration_compareVariants/test.py
except:
- master
- release
nonlocal_densityConservation:
stage: spectral
stage: grid
script: nonlocal_densityConservation/test.py
except:
- master
- release
Spectral_ipNeighborhood:
stage: spectral
stage: grid
script: Spectral_ipNeighborhood/test.py
except:
- master
- release
RGC_DetectChanges:
stage: spectral
stage: grid
script: RGC_DetectChanges/test.py
except:
- master
- release
Nonlocal_Damage_DetectChanges:
stage: spectral
stage: grid
script: Nonlocal_Damage_DetectChanges/test.py
except:
- master
- release
SpectralAll_restart:
stage: spectral
script: SpectralAll_restart/test.py
grid_all_restart:
stage: grid
script: grid_all_restart/test.py
except:
- master
- release
SpectralAll_parsingLoadCase:
stage: spectral
script: SpectralAll_parsingLoadCase/test.py
grid_parsingLoadCase:
stage: grid
script: grid_parsingLoadCase/test.py
except:
- master
- release
SpectralBasic_loadCaseRotation:
stage: spectral
script: SpectralBasic_loadCaseRotation/test.py
grid_all_loadCaseRotation:
stage: grid
script: grid_all_loadCaseRotation/test.py
except:
- master
- release
Spectral_MPI:
stage: spectral
grid_mech_MPI:
stage: grid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- Spectral_MPI/test.py
- grid_mech_MPI/test.py
except:
- master
- release
SpectralAll_restartMPI:
stage: spectral
grid_all_restartMPI:
stage: grid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- SpectralAll_restartMPI/test.py
- grid_all_restartMPI/test.py
except:
- master
- release
Plasticity_DetectChanges:
stage: spectral
stage: grid
script: Plasticity_DetectChanges/test.py
except:
- master
- release
Homogenization:
stage: spectral
stage: grid
script: Homogenization/test.py
except:
- master
- release
Phenopowerlaw_singleSlip:
stage: spectral
stage: grid
script: Phenopowerlaw_singleSlip/test.py
except:
- master
- release
TextureComponents:
stage: spectral
stage: grid
script: TextureComponents/test.py
except:
- master
@ -451,9 +451,9 @@ Abaqus_compile:
- release
###################################################################################################
SpectralExample:
grid_all_example:
stage: example
script: SpectralAll_example/test.py
script: grid_all_example/test.py
only:
- development
@ -463,7 +463,7 @@ SpectralRuntime:
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd $DAMASKROOT
- make clean spectral processing OPTIMIZATION=AGGRESSIVE
- make clean grid processing OPTIMIZATION=AGGRESSIVE
- cd $TESTROOT/performance # location of old results
- git checkout . # undo any changes (i.e. run time data from non-development branch)
- cd $DAMASKROOT/PRIVATE/testing
@ -501,7 +501,7 @@ Marc:
- master
- release
Spectral:
GridSolver:
stage: createDocumentation
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen

4
CMakeLists.txt
View File

@ -105,9 +105,9 @@ set (CMAKE_C_COMPILER "${PETSC_MPICC}")
# Now start to care about DAMASK
# DAMASK solver defines project to build
if (DAMASK_SOLVER STREQUAL "SPECTRAL")
if (DAMASK_SOLVER STREQUAL "GRID")
project (DAMASK_spectral Fortran C)
add_definitions (-DSpectral)
add_definitions (-DGrid)
message ("Building Spectral Solver\n")
elseif (DAMASK_SOLVER STREQUAL "FEM")
project (DAMASK_FEM Fortran C)

18
Makefile
View File

@ -3,20 +3,24 @@ SHELL = /bin/sh
# Makefile for the installation of DAMASK
########################################################################################
.PHONY: all
all: spectral FEM processing
all: grid FEM processing
.PHONY: grid
grid: build/grid
@(cd build/grid;make -j4 --no-print-directory -ws all install;)
.PHONY: spectral
spectral: build/spectral
@(cd build/spectral;make -j4 --no-print-directory -ws all install;)
spectral: build/grid
@(cd build/grid;make -j4 --no-print-directory -ws all install;)
.PHONY: FEM
FEM: build/FEM
@(cd build/FEM; make -j4 --no-print-directory -ws all install;)
.PHONY: build/spectral
build/spectral:
@mkdir -p build/spectral
@(cd build/spectral; cmake -Wno-dev -DDAMASK_SOLVER=SPECTRAL -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP} ../../;)
.PHONY: build/grid
build/grid:
@mkdir -p build/grid
@(cd build/grid; cmake -Wno-dev -DDAMASK_SOLVER=GRID -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DBUILDCMD_POST=${BUILDCMD_POST} -DBUILDCMD_PRE=${BUILDCMD_PRE} -DOPTIMIZATION=${OPTIMIZATION} -DOPENMP=${OPENMP} ../../;)
.PHONY: build/FEM
build/FEM:

2
PRIVATE

@ -1 +1 @@
Subproject commit d81a446bdfaa2bc3c939e802c50a5fd8f2fb38e3
Subproject commit 397d9265ef677966610831bbf4d1358d879a4ac2

10
installation/mods_Abaqus/abaqus_v6.env
View File

@ -11,11 +11,8 @@
# Compile_cpp and link_exe for Abaqus make utility.
#
import os, re, glob, driverUtils
from damask import version as DAMASKVERSION
from damask import Environment
myEnv = Environment()
if myEnv.options['DAMASK_HDF5'] == 'ON':
if false:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]
@ -44,7 +41,7 @@ compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " +
"-implicitnone -standard-semantics " +
"-assume nostd_mod_proc_name " +
"-real-size 64 " +
'-DDAMASKVERSION=\\\"%s\\\"'%DAMASKVERSION)
'-DDAMASKVERSION=\\\"n/a\\\"')
# Abaqus/CAE will generate an input file without parts and assemblies.
cae_no_parts_input_file=ON
@ -57,6 +54,3 @@ ask_delete=OFF
# Remove the temporary names from the namespace
del fortCmd
del Environment
del myEnv
del DAMASKVERSION

10
installation/mods_Abaqus/abaqus_v6_debug.env
View File

@ -11,11 +11,8 @@
# Compile_cpp and link_exe for Abaqus make utility.
#
import os, re, glob, driverUtils
from damask import version as DAMASKVERSION
from damask import Environment
myEnv = Environment()
if myEnv.options['DAMASK_HDF5'] == 'ON':
if false:
# use hdf5 compiler wrapper in $PATH
fortCmd = os.popen('h5fc -shlib -show').read().replace('\n','') # complicated way needed to pass in DAMASKVERSION string
link_sl += fortCmd.split()[1:]
@ -49,7 +46,7 @@ compile_fortran = (fortCmd + " -c -fPIC -auto -shared-intel " +
"-check bounds,format,output_conversion,uninit " +
"-ftrapuv -fpe-all0 " +
"-g -traceback -gen-interfaces -fp-stack-check -fp-model strict " +
'-DDAMASKVERSION=\\\"%s\\\"'%DAMASKVERSION)
'-DDAMASKVERSION=\\\"n/a\\\"')
# Abaqus/CAE will generate an input file without parts and assemblies.
cae_no_parts_input_file=ON
@ -62,6 +59,3 @@ ask_delete=OFF
# Remove the temporary names from the namespace
del fortCmd
del Environment
del myEnv
del DAMASKVERSION

1
src/CMakeLists.txt
View File

@ -176,6 +176,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
add_library(SPECTRAL_SOLVER OBJECT
"grid_thermal_spectral.f90"
"grid_damage_spectral.f90"
"grid_mech_FEM.f90"
"grid_mech_spectral_basic.f90"
"grid_mech_spectral_polarisation.f90")
add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)

19
src/DAMASK_grid.f90
View File

@ -31,6 +31,8 @@ program DAMASK_spectral
IO_lc, &
IO_intOut, &
IO_warning
use config, only: &
config_numerics
use debug, only: &
debug_level, &
debug_spectral, &
@ -50,7 +52,6 @@ program DAMASK_spectral
worldsize, &
stagItMax, &
maxCutBack, &
spectral_solver, &
continueCalculation
use homogenization, only: &
materialpoint_sizeResults, &
@ -73,6 +74,7 @@ program DAMASK_spectral
FIELD_DAMAGE_ID
use grid_mech_spectral_basic
use grid_mech_spectral_polarisation
use grid_mech_FEM
use grid_damage_spectral
use grid_thermal_spectral
use results
@ -165,21 +167,28 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! assign mechanics solver depending on selected type
select case (spectral_solver)
case (GRID_MECH_SPECTRAL_BASIC_LABEL)
select case (trim(config_numerics%getString('spectral_solver',defaultVal='basic')))
case ('basic')
mech_init => grid_mech_spectral_basic_init
mech_forward => grid_mech_spectral_basic_forward
mech_solution => grid_mech_spectral_basic_solution
case (GRID_MECH_SPECTRAL_POLARISATION_LABEL)
case ('polarisation')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
mech_init => grid_mech_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution
case ('fem')
if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
call IO_warning(42_pInt, ext_msg='debug Divergence')
mech_init => grid_mech_FEM_init
mech_forward => grid_mech_FEM_forward
mech_solution => grid_mech_FEM_solution
case default
call IO_error(error_ID = 891_pInt, ext_msg = trim(spectral_solver))
call IO_error(error_ID = 891_pInt, ext_msg = config_numerics%getString('spectral_solver'))
end select

2
src/FEsolving.f90
View File

@ -72,7 +72,7 @@ subroutine FE_init
modelName = getSolverJobName()
#if defined(Spectral) || defined(FEM)
#if defined(Grid) || defined(FEM)
restartInc = interface_RestartInc
if(restartInc < 0_pInt) then

2
src/IO.f90
View File

@ -823,8 +823,6 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'microstructure count mismatch'
case (846_pInt)
msg = 'rotation for load case rotation ill-defined (R:RT != I)'
case (847_pInt)
msg = 'update of gamma operator not possible when pre-calculated'
case (880_pInt)
msg = 'mismatch of microstructure count and a*b*c in geom file'
case (891_pInt)

4
src/config.f90
View File

@ -608,7 +608,7 @@ character(len=65536) function getString(this,key,defaultVal,raw)
implicit none
class(tPartitionedStringList), target, intent(in) :: this
character(len=*), intent(in) :: key
character(len=65536), intent(in), optional :: defaultVal
character(len=*), intent(in), optional :: defaultVal
logical, intent(in), optional :: raw
type(tPartitionedStringList), pointer :: item
logical :: found, &
@ -622,7 +622,7 @@ character(len=65536) function getString(this,key,defaultVal,raw)
found = present(defaultVal)
if (found) then
getString = trim(defaultVal)
if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0,ext_msg='getString')
!if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0,ext_msg='getString')
endif
item => this

6
src/grid_damage_spectral.f90
View File

@ -170,7 +170,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
loadCaseTime !< remaining time of current load case
integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt ::position
PetscInt :: devNull
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
PetscErrorCode :: ierr
@ -208,8 +208,8 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
call VecMin(solution_vec,position,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,position,maxDamage,ierr); CHKERRQ(ierr)
call VecMin(solution_vec,devNull,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution_vec,devNull,maxDamage,ierr); CHKERRQ(ierr)
if (solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&

737
src/grid_mech_FEM.f90 Normal file
View File

@ -0,0 +1,737 @@
!--------------------------------------------------------------------------------------------------
!> @author Arko Jyoti Bhattacharjee, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: FEM
!--------------------------------------------------------------------------------------------------
module grid_mech_FEM
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
use spectral_utilities, only: &
tSolutionState, &
tSolutionParams
implicit none
private
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: mech_grid
SNES, private :: mech_snes
Vec, private :: solution_current, solution_lastInc, solution_rate
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc
real(pReal), private :: detJ
real(pReal), private, dimension(3) :: delta
real(pReal), private, dimension(3,8) :: BMat
real(pReal), private, dimension(8,8) :: HGMat
PetscInt, private :: xstart,ystart,zstart,xend,yend,zend
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastIter = math_I3, &
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
real(pReal), private :: &
err_BC !< deviation from stress BC
integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment
public :: &
grid_mech_FEM_init, &
grid_mech_FEM_solution, &
grid_mech_FEM_forward
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_init
use IO, only: &
IO_intOut, &
IO_error, &
IO_open_jobFile_binary
use FEsolving, only: &
restartInc
use numerics, only: &
worldrank, &
worldsize, &
petsc_options
use homogenization, only: &
materialpoint_F0
use DAMASK_interface, only: &
getSolverJobName
use spectral_utilities, only: &
utilities_constitutiveResponse, &
utilities_updateIPcoords, &
wgt
use mesh, only: &
geomSize, &
grid, &
grid3
use math, only: &
math_invSym3333
implicit none
real(pReal) :: HGCoeff = 0e-2_pReal
PetscInt, dimension(:), allocatable :: localK
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
real(pReal), dimension(4,8) :: &
HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, &
1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, &
-1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, &
-1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, &
-1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, &
-1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
PetscErrorCode :: ierr
integer(pInt) :: rank
integer :: fileUnit
character(len=1024) :: rankStr
real(pReal), dimension(3,3,3,3) :: devNull
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
write(6,'(/,a)') ' <<<+- grid_mech_FEM init -+>>>'
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres &
&-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (P_current (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do rank = 1, worldsize
call MPI_Bcast(localK(rank),1,MPI_INTEGER,rank-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, &
DMDA_STENCIL_BOX, &
grid(1),grid(2),grid(3), &
1, 1, worldsize, &
3, 1, &
[grid(1)],[grid(2)],localK, &
mech_grid,ierr)
CHKERRQ(ierr)
call DMDASetUniformCoordinates(mech_grid,0.0,geomSize(1),0.0,geomSize(2),0.0,geomSize(3),ierr)
CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr)
call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr)
call DMsetUp(mech_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr)
call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) ! specify custom convergence check function "_converged"
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
!--------------------------------------------------------------------------------------------------
! init fields
call VecSet(solution_current,0.0,ierr);CHKERRQ(ierr)
call VecSet(solution_lastInc,0.0,ierr);CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent
CHKERRQ(ierr)
xend = xstart+xend-1
yend = ystart+yend-1
zend = zstart+zend-1
delta = geomSize/real(grid,pReal) ! grid spacing
detJ = product(delta) ! cell volume
BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
-1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), &
1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), &
-1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), &
1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), &
-1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), &
1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix
HGMat = matmul(transpose(HGcomp),HGcomp) &
* HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix
!--------------------------------------------------------------------------------------------------
! init fields
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('F_aim')
read(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
read(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u'//trim(rankStr))
read(fileUnit) u_current; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr))
read(fileUnit) u_lastInc; close (fileUnit)
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restart
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
restartRead: if (restartInc > 0_pInt) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
endif restartRead
end subroutine grid_mech_FEM_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
use IO, only: &
IO_error
use spectral_utilities, only: &
tBoundaryCondition, &
utilities_maskedCompliance
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
incInfo = incInfoIn
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr);CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(mech_snes,reason,ierr);CHKERRQ(ierr)
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
end function grid_mech_FEM_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use homogenization, only: &
materialpoint_F0
use mesh, only: &
grid, &
grid3
use CPFEM2, only: &
CPFEM_age
use spectral_utilities, only: &
utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
integer :: fileUnit
character(len=32) :: rankStr
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim','w')
write(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
write(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u'//trim(rankStr),'w')
write(fileUnit) u_current; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr),'w')
write(fileUnit) u_lastInc; close (fileUnit)
endif
call CPFEM_age() ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif
if (guess) then
call VecWAXPY(solution_rate,-1.0,solution_lastInc,solution_current,ierr)
CHKERRQ(ierr)
call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr)
else
call VecSet(solution_rate,0.0,ierr); CHKERRQ(ierr)
endif
call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr)
F_lastInc = F ! winding F forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
end subroutine grid_mech_FEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use mesh
use spectral_utilities
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
err_div, &
divTol, &
BCTol
err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ
divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then
reason = -1
else
reason = 0
endif
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(da_local,x_local,f_local,dummy,ierr)
use numerics, only: &
itmax, &
itmin
use numerics, only: &
worldrank
use mesh, only: &
grid
use math, only: &
math_rotate_backward33, &
math_mul3333xx33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: &
utilities_constitutiveResponse
use IO, only: &
IO_intOut
use FEsolving, only: &
terminallyIll
use homogenization, only: &
materialpoint_dPdF
implicit none
DM :: da_local
Vec :: x_local, f_local
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
PetscScalar, dimension(8,3) :: x_elem, f_elem
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
real(pReal), dimension(3,3) :: &
deltaF_aim
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal), dimension(3,3,3,3) :: devNull
call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1_pInt
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim)
flush(6)
endif newIteration
!--------------------------------------------------------------------------------------------------
! get deformation gradient
call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem))
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(P_current,&
P_av,C_volAvg,devNull, &
F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim_lastIter = F_aim
deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc
!--------------------------------------------------------------------------------------------------
! constructing residual
call VecSet(f_local,0.0,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
enddo; enddo; enddo
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + &
matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1
f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3)
enddo; enddo; enddo
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
if (zstart == 0) then
f_scal(0:2,xstart,ystart,zstart) = 0.0
f_scal(0:2,xend+1,ystart,zstart) = 0.0
f_scal(0:2,xstart,yend+1,zstart) = 0.0
f_scal(0:2,xend+1,yend+1,zstart) = 0.0
endif
if (zend + 1 == grid(3)) then
f_scal(0:2,xstart,ystart,zend+1) = 0.0
f_scal(0:2,xend+1,ystart,zend+1) = 0.0
f_scal(0:2,xstart,yend+1,zend+1) = 0.0
f_scal(0:2,xend+1,yend+1,zend+1) = 0.0
endif
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
end subroutine formResidual
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM stiffness matrix
!--------------------------------------------------------------------------------------------------
subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
use mesh, only: &
mesh_ipCoordinates
use homogenization, only: &
materialpoint_dPdF
implicit none
DM :: da_local
Vec :: x_local, coordinates
Mat :: Jac_pre, Jac
MatStencil,dimension(4,24) :: row, col
PetscScalar,pointer,dimension(:,:,:,:) :: x_scal
PetscScalar,dimension(24,24) :: K_ele
PetscScalar,dimension(9,24) :: BMatFull
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
PetscInt,dimension(3) :: rows
PetscScalar :: diag
PetscObject :: dummy
MatNullSpace :: matnull
PetscErrorCode :: ierr
BMatFull = 0.0
BMatFull(1:3,1 :8 ) = BMat
BMatFull(4:6,9 :16) = BMat
BMatFull(7:9,17:24) = BMat
call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr)
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1
col(MatStencil_i,ctr ) = i+ii
col(MatStencil_j,ctr ) = j+jj
col(MatStencil_k,ctr ) = k+kk
col(MatStencil_c,ctr ) = 0
col(MatStencil_i,ctr+8 ) = i+ii
col(MatStencil_j,ctr+8 ) = j+jj
col(MatStencil_k,ctr+8 ) = k+kk
col(MatStencil_c,ctr+8 ) = 1
col(MatStencil_i,ctr+16) = i+ii
col(MatStencil_j,ctr+16) = j+jj
col(MatStencil_k,ctr+16) = k+kk
col(MatStencil_c,ctr+16) = 2
enddo; enddo; enddo
row = col
ele = ele + 1
K_ele = 0.0
K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele = K_ele + &
matmul(transpose(BMatFull), &
matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), &
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr)
CHKERRQ(ierr)
enddo; enddo; enddo
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
rows = [0, 1, 2]
diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + &
C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + &
C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ
call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
CHKERRQ(ierr)
call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
ele = ele + 1
x_scal(0:2,i,j,k) = mesh_ipCoordinates(1:3,1,ele)
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates)
call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes
call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
end subroutine formJacobian
end module grid_mech_FEM

413
src/grid_mech_spectral_basic.f90
View File

@ -20,13 +20,17 @@ module grid_mech_spectral_basic
implicit none
private
character (len=*), parameter, public :: &
GRID_MECH_SPECTRAL_BASIC_LABEL = 'basic'
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
type, private :: tNumerics
logical :: &
update_gamma !< update gamma operator with current stiffness
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: da
@ -35,7 +39,9 @@ module grid_mech_spectral_basic
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, &
Fdot
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
@ -46,7 +52,6 @@ module grid_mech_spectral_basic
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
@ -65,6 +70,9 @@ module grid_mech_spectral_basic
grid_mech_spectral_basic_init, &
grid_mech_spectral_basic_solution, &
grid_mech_spectral_basic_forward
private :: &
converged, &
formResidual
contains
@ -76,12 +84,10 @@ subroutine grid_mech_spectral_basic_init
IO_intOut, &
IO_error, &
IO_open_jobFile_binary
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use config, only :&
config_numerics
use numerics, only: &
worldrank, &
worldsize, &
@ -107,7 +113,8 @@ subroutine grid_mech_spectral_basic_init
temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F
PetscScalar, pointer, dimension(:,:,:,:) :: &
F ! pointer to solution data
PetscInt, dimension(worldsize) :: localK
integer :: fileUnit
character(len=1024) :: rankStr
@ -120,6 +127,8 @@ subroutine grid_mech_spectral_basic_init
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
@ -152,23 +161,23 @@ subroutine grid_mech_spectral_basic_init
call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr)! residual vector of same shape as solution vector
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,grid_mech_spectral_basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "_converged"
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "converged"
CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
restart: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading values of increment ', restartInc, ' from file'
flush(6)
endif
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('F_aim')
read(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
read(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
@ -179,12 +188,6 @@ subroutine grid_mech_spectral_basic_init
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='F_aim')
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='F_aim_lastInc')
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
@ -196,14 +199,10 @@ subroutine grid_mech_spectral_basic_init
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
! QUESTION: why not writing back right after reading (l.189)?
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) &
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file'
flush(6)
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
@ -218,11 +217,9 @@ end subroutine grid_mech_spectral_basic_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic scheme with internal iterations
!> @brief solution for the basic scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
use numerics, only: &
update_gamma
use spectral_utilities, only: &
tBoundaryCondition, &
utilities_maskedCompliance, &
@ -245,7 +242,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
real(pReal), dimension(3,3), intent(in) :: rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
@ -256,7 +252,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
if (num%update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
!--------------------------------------------------------------------------------------------------
! set module wide available data
@ -279,14 +275,187 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
solution%termIll = terminallyIll
terminallyIll = .false.
end function grid_mech_spectral_basic_solution
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use homogenization, only: &
materialpoint_F0
use mesh, only: &
grid, &
grid3
use CPFEM2, only: &
CPFEM_age
use spectral_utilities, only: &
utilities_calculateRate, &
utilities_forwardField, &
utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer :: fileUnit
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim','w')
write(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
write(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
endif
call CPFEM_age ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_forward
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
divTol, &
BCTol
divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then
reason = -1
else
reason = 0
endif
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_formResidual(in, F, &
subroutine formResidual(in, F, &
residuum, dummy, ierr)
use numerics, only: &
itmax, &
@ -371,175 +540,7 @@ subroutine grid_mech_spectral_basic_formResidual(in, F, &
! constructing residual
residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
end subroutine grid_mech_spectral_basic_formResidual
end subroutine formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
divTol, &
BCTol
divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then converged
reason = -1
else converged
reason = 0
endif converged
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
write(6,'(a,f12.2,a,es8.2,a,es9.2,a)') ' error stress BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine grid_mech_spectral_basic_converged
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use numerics, only: &
worldrank
use homogenization, only: &
materialpoint_F0
use mesh, only: &
grid, &
grid3
use CPFEM2, only: &
CPFEM_age
use spectral_utilities, only: &
utilities_calculateRate, &
utilities_forwardField, &
utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
IO_open_jobFile_binary
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer :: fileUnit
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
endif
call CPFEM_age ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif
Fdot = Utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_forward
end module grid_mech_spectral_basic

549
src/grid_mech_spectral_polarisation.f90
View File

@ -2,7 +2,7 @@
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Polarisation scheme solver
!> @brief Grid solver for mechanics: Spectral Polarisation
!--------------------------------------------------------------------------------------------------
module grid_mech_spectral_polarisation
#include <petsc/finclude/petscsnes.h>
@ -10,7 +10,6 @@ module grid_mech_spectral_polarisation
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
@ -21,13 +20,17 @@ module grid_mech_spectral_polarisation
implicit none
private
character (len=*), parameter, public :: &
GRID_MECH_SPECTRAL_POLARISATION_LABEL = 'polarisation'
!--------------------------------------------------------------------------------------------------
! derived types
type(tSolutionParams), private :: params
type, private :: tNumerics
logical :: &
update_gamma !< update gamma operator with current stiffness
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: da
@ -49,8 +52,8 @@ module grid_mech_spectral_polarisation
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
P_avLastEval = 0.0_pReal !< average 1st Piola--Kirchhoff stress last call of CPFEM_general
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
@ -66,13 +69,16 @@ module grid_mech_spectral_polarisation
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment
integer, private :: &
totalIter = 0 !< total iteration in current increment
public :: &
grid_mech_spectral_polarisation_init, &
grid_mech_spectral_polarisation_solution, &
grid_mech_spectral_polarisation_forward
private :: &
converged, &
formResidual
contains
@ -84,12 +90,10 @@ subroutine grid_mech_spectral_polarisation_init
IO_intOut, &
IO_error, &
IO_open_jobFile_binary
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
use FEsolving, only: &
restartInc
use config, only :&
config_numerics
use numerics, only: &
worldrank, &
worldsize, &
@ -99,9 +103,9 @@ subroutine grid_mech_spectral_polarisation_init
use DAMASK_interface, only: &
getSolverJobName
use spectral_utilities, only: &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
Utilities_updateIPcoords, &
utilities_constitutiveResponse, &
utilities_updateGamma, &
utilities_updateIPcoords, &
wgt
use mesh, only: &
grid, &
@ -123,11 +127,13 @@ subroutine grid_mech_spectral_polarisation_init
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
@ -164,7 +170,7 @@ subroutine grid_mech_spectral_polarisation_init
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,grid_mech_spectral_polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged"
CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
@ -173,13 +179,14 @@ subroutine grid_mech_spectral_polarisation_init
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:)
restart: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading values of increment ', restartInc, ' from file'
flush(6)
endif
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('F_aim')
read(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
read(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
@ -189,19 +196,12 @@ subroutine grid_mech_spectral_polarisation_init
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr))
read(fileUnit) F_tau; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr))
read(fileUnit) F_tau_lastInc; close (fileUnit)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
call MPI_Allreduce(MPI_IN_PLACE,F_aim,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim')
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
call MPI_Allreduce(MPI_IN_PLACE,F_aim_lastInc,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='F_aim_lastInc')
elseif (restartInc == 0_pInt) then restart
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_tau = 2.0_pReal*F
@ -214,15 +214,10 @@ subroutine grid_mech_spectral_polarisation_init
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
nullify(F)
nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! write data back to PETSc
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead: if (restartInc > 0_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file'
flush(6)
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
@ -242,16 +237,12 @@ end subroutine grid_mech_spectral_polarisation_init
!> @brief solution for the Polarisation scheme with internal iterations
!--------------------------------------------------------------------------------------------------
function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
use IO, only: &
IO_error
use numerics, only: &
update_gamma
use math, only: &
math_invSym3333
use spectral_utilities, only: &
tBoundaryCondition, &
Utilities_maskedCompliance, &
Utilities_updateGamma
utilities_maskedCompliance, &
utilities_updateGamma
use FEsolving, only: &
restartWrite, &
terminallyIll
@ -280,14 +271,14 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (update_gamma) then
call Utilities_updateGamma(C_minMaxAvg,restartWrite)
if (num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg,restartWrite)
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
endif
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
! set module wide available data
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
@ -306,232 +297,10 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
if (reason == -4) call IO_error(893_pInt) ! MPI error
end function grid_mech_spectral_polarisation_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the Polarisation residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work)
FandF_tau, & ! defgrad fields on grid
residuum, & ! residuum fields on grid
dummy, &
ierr)
use numerics, only: &
itmax, &
itmin, &
polarAlpha, &
polarBeta
use mesh, only: &
grid, &
grid3
use IO, only: &
IO_intOut
use math, only: &
math_rotate_backward33, &
math_mul3333xx33, &
math_invSym3333, &
math_mul33x33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: &
wgt, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS, &
Utilities_curlRMS
use homogenization, only: &
materialpoint_dPdF
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in
PetscScalar, &
target, dimension(3,3,2, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: FandF_tau
PetscScalar, &
target, dimension(3,3,2, X_RANGE, Y_RANGE, Z_RANGE), intent(out) :: residuum
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, &
F_tau, &
residual_F, &
residual_F_tau
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
integer(pInt) :: &
i, j, k, e
F => FandF_tau(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => FandF_tau(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => residuum(1:3,1:3,1,&
X_RANGE, Y_RANGE, Z_RANGE)
residual_F_tau => residuum(1:3,1:3,2,&
X_RANGE, Y_RANGE, Z_RANGE)
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)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1_pInt
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim)
flush(6)
endif newIteration
!--------------------------------------------------------------------------------------------------
!
tensorField_real = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
tensorField_real(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))
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call utilities_FFTtensorForward()
call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call utilities_FFTtensorBackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
P_avLastEval = P_av
call Utilities_constitutiveResponse(residual_F,P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------
! calculate divergence
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise
call utilities_FFTtensorForward()
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0_pInt
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(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) - 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))) &
+ residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! calculating curl
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward()
err_curl = Utilities_curlRMS()
nullify(F)
nullify(F_tau)
nullify(residual_F)
nullify(residual_F_tau)
end subroutine formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_curl_tolRel, &
err_curl_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use math, only: &
math_mul3333xx33
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, &
snorm, &
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
curlTol, &
divTol, &
BCTol
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim-F_av) + &
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
!--------------------------------------------------------------------------------------------------
! error calculation
curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs)
divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs)
converged: if ((totalIter >= itmin .and. &
all([ err_div /divTol, &
err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then converged
reason = -1
else converged
reason = 0
endif converged
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine grid_mech_spectral_polarisation_converged
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep
@ -552,9 +321,9 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
use CPFEM2, only: &
CPFEM_age
use spectral_utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_updateIPcoords, &
utilities_calculateRate, &
utilities_forwardField, &
utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use IO, only: &
@ -576,14 +345,12 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
integer(pInt) :: i, j, k
integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
integer :: fileUnit
character(len=32) :: rankStr
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:)
@ -603,6 +370,10 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim','w')
write(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
write(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
@ -612,14 +383,13 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w')
write(fileUnit) F_tau; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w')
write(fileUnit) F_tau_lastInc; close (fileUnit)
endif
call CPFEM_age() ! age state and kinematics
call CPFEM_age ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
@ -642,10 +412,10 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
endif
Fdot = Utilities_calculateRate(guess, &
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC))
F_tauDot = Utilities_calculateRate(guess, &
F_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
@ -656,15 +426,14 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),&
[9,grid(1),grid(2),grid3])
if (guess) then
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
[9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
else
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(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,&
@ -675,10 +444,218 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
enddo; enddo; enddo
endif
nullify(F)
nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_forward
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tolRel, &
err_div_tolAbs, &
err_curl_tolRel, &
err_curl_tolAbs, &
err_stress_tolRel, &
err_stress_tolAbs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
real(pReal) :: &
curlTol, &
divTol, &
BCTol
curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs)
divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs)
if ((totalIter >= itmin .and. &
all([ err_div /divTol, &
err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then
reason = -1
else
reason = 0
endif
!--------------------------------------------------------------------------------------------------
! report
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ==========================================================================='
flush(6)
end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the polarisation residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, &
residuum, dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
polarAlpha, &
polarBeta
use mesh, only: &
grid, &
grid3
use math, only: &
math_rotate_forward33, &
math_rotate_backward33, &
math_mul3333xx33, &
math_invSym3333, &
math_mul33x33
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRotation
use spectral_utilities, only: &
wgt, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, &
utilities_constitutiveResponse, &
utilities_divergenceRMS, &
utilities_curlRMS
use IO, only: &
IO_intOut
use homogenization, only: &
materialpoint_dPdF
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), &
target, intent(in) :: FandF_tau
PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),&
target, intent(out) :: residuum !< residuum field
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, &
F_tau, &
residual_F, &
residual_F_tau
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
integer :: &
i, j, k, e
F => FandF_tau(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => FandF_tau(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => residuum(1:3,1:3,1,&
X_RANGE, Y_RANGE, Z_RANGE)
residual_F_tau => residuum(1:3,1:3,2,&
X_RANGE, Y_RANGE, Z_RANGE)
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)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim)
flush(6)
endif newIteration
!--------------------------------------------------------------------------------------------------
!
tensorField_real = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
tensorField_real(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))
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call utilities_FFTtensorForward
call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
-math_rotate_forward33(F_av,params%rotation_BC)) + &
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
! calculate divergence
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1
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) - 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))) &
+ residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! calculating curl
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward
err_curl = Utilities_curlRMS()
end subroutine formResidual
end module grid_mech_spectral_polarisation

81
src/numerics.f90
View File

@ -77,7 +77,7 @@ module numerics
!--------------------------------------------------------------------------------------------------
! spectral parameters:
#ifdef Spectral
#ifdef Grid
real(pReal), protected, public :: &
err_div_tolAbs = 1.0e-4_pReal, & !< absolute tolerance for equilibrium
err_div_tolRel = 5.0e-4_pReal, & !< relative tolerance for equilibrium
@ -85,27 +85,17 @@ module numerics
err_curl_tolRel = 5.0e-4_pReal, & !< relative tolerance for compatibility
err_stress_tolAbs = 1.0e3_pReal, & !< absolute tolerance for fullfillment of stress BC
err_stress_tolRel = 0.01_pReal, & !< relative tolerance for fullfillment of stress BC
fftw_timelimit = -1.0_pReal, & !< sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
rotation_tol = 1.0e-12_pReal, & !< tolerance of rotation specified in loadcase, Default 1.0e-12: first guess
polarAlpha = 1.0_pReal, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
polarBeta = 1.0_pReal !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
character(len=64), private :: &
fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
character(len=64), protected, public :: &
spectral_solver = 'basic', & !< spectral solution method
spectral_derivative = 'continuous' !< spectral spatial derivative method
character(len=1024), protected, public :: &
petsc_defaultOptions = '-mech_snes_type ngmres &
&-damage_snes_type ngmres &
&-thermal_snes_type ngmres ', &
petsc_options = ''
integer(pInt), protected, public :: &
fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw
divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: size scaled to 1, 2: size scaled to Npoints
logical, protected, public :: &
continueCalculation = .false., & !< false:exit if BVP solver does not converge, true: continue calculation despite BVP solver not converging
memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness
continueCalculation = .false. !< false:exit if BVP solver does not converge, true: continue calculation despite BVP solver not converging
#endif
!--------------------------------------------------------------------------------------------------
@ -319,7 +309,7 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Spectral
#ifdef Grid
case ('err_div_tolabs')
err_div_tolAbs = IO_floatValue(line,chunkPos,2_pInt)
case ('err_div_tolrel')
@ -330,22 +320,8 @@ subroutine numerics_init
err_stress_tolabs = IO_floatValue(line,chunkPos,2_pInt)
case ('continuecalculation')
continueCalculation = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('memory_efficient')
memory_efficient = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('fftw_timelimit')
fftw_timelimit = IO_floatValue(line,chunkPos,2_pInt)
case ('fftw_plan_mode')
fftw_plan_mode = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case ('spectralderivative')
spectral_derivative = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case ('divergence_correction')
divergence_correction = IO_intValue(line,chunkPos,2_pInt)
case ('update_gamma')
update_gamma = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
case ('petsc_options')
petsc_options = trim(line(chunkPos(4):))
case ('spectralsolver','myspectralsolver')
spectral_solver = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
case ('err_curl_tolabs')
err_curl_tolAbs = IO_floatValue(line,chunkPos,2_pInt)
case ('err_curl_tolrel')
@ -354,13 +330,6 @@ subroutine numerics_init
polarAlpha = IO_floatValue(line,chunkPos,2_pInt)
case ('polarbeta')
polarBeta = IO_floatValue(line,chunkPos,2_pInt)
#else
case ('err_div_tolabs','err_div_tolrel','err_stress_tolrel','err_stress_tolabs',& ! found spectral parameter for FEM build
'memory_efficient','fftw_timelimit','fftw_plan_mode', &
'divergence_correction','update_gamma','spectralfilter','myfilter', &
'err_curl_tolabs','err_curl_tolrel', &
'polaralpha','polarbeta')
call IO_warning(40_pInt,ext_msg=tag)
#endif
!--------------------------------------------------------------------------------------------------
@ -374,36 +343,16 @@ subroutine numerics_init
petsc_options = trim(line(chunkPos(4):))
case ('bbarstabilisation')
BBarStabilisation = IO_intValue(line,chunkPos,2_pInt) > 0_pInt
#else
case ('integrationorder','structorder','thermalorder', 'damageorder', &
'bbarstabilisation')
call IO_warning(40_pInt,ext_msg=tag)
#endif
case default ! found unknown keyword
call IO_error(300_pInt,ext_msg=tag)
end select
enddo
else fileExists
write(6,'(a,/)') ' using standard values'
flush(6)
endif fileExists
#ifdef Spectral
select case(IO_lc(fftw_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
fftw_planner_flag = 64_pInt
case('measure','fftw_measure')
fftw_planner_flag = 0_pInt
case('patient','fftw_patient')
fftw_planner_flag= 32_pInt
case('exhaustive','fftw_exhaustive')
fftw_planner_flag = 8_pInt
case default
call IO_warning(warning_ID=47_pInt,ext_msg=trim(IO_lc(fftw_plan_mode)))
fftw_planner_flag = 32_pInt
end select
#endif
!--------------------------------------------------------------------------------------------------
! writing parameters to output
@ -478,19 +427,8 @@ subroutine numerics_init
!--------------------------------------------------------------------------------------------------
! spectral parameters
#ifdef Spectral
#ifdef Grid
write(6,'(a24,1x,L8)') ' continueCalculation: ',continueCalculation
write(6,'(a24,1x,L8)') ' memory_efficient: ',memory_efficient
write(6,'(a24,1x,i8)') ' divergence_correction: ',divergence_correction
write(6,'(a24,1x,a)') ' spectral_derivative: ',trim(spectral_derivative)
if(fftw_timelimit<0.0_pReal) then
write(6,'(a24,1x,L8)') ' fftw_timelimit: ',.false.
else
write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit
endif
write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode)
write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,1x,L8,/)') ' update_gamma: ',update_gamma
write(6,'(a24,1x,es8.1)') ' err_stress_tolAbs: ',err_stress_tolAbs
write(6,'(a24,1x,es8.1)') ' err_stress_tolRel: ',err_stress_tolRel
write(6,'(a24,1x,es8.1)') ' err_div_tolAbs: ',err_div_tolAbs
@ -499,7 +437,6 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel
write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha
write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta
write(6,'(a24,1x,a)') ' spectral solver: ',trim(spectral_solver)
write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_defaultOptions)//' '//trim(petsc_options)
#endif
@ -564,11 +501,7 @@ subroutine numerics_init
if (err_thermal_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_thermal_tolrel')
if (err_damage_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_damage_tolabs')
if (err_damage_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_damage_tolrel')
#ifdef Spectral
if (divergence_correction < 0_pInt .or. &
divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction')
if (update_gamma .and. &
.not. memory_efficient) call IO_error(error_ID = 847_pInt)
#ifdef Grid
if (err_stress_tolrel <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolRel')
if (err_stress_tolabs <= 0.0_pReal) call IO_error(301_pInt,ext_msg='err_stress_tolAbs')
if (err_div_tolRel < 0.0_pReal) call IO_error(301_pInt,ext_msg='err_div_tolRel')

5
src/prec.f90
View File

@ -13,8 +13,11 @@ module prec
implicit none
private
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
#ifdef Abaqus
integer, parameter, public :: pReal = selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
#else
integer, parameter, public :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
#endif
#if(INT==8)
integer, parameter, public :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit)
#else

331
src/spectral_utilities.f90
View File

@ -9,7 +9,7 @@ module spectral_utilities
use PETScSys
use prec, only: &
pReal, &
pInt
pStringLen
use math, only: &
math_I3
@ -18,17 +18,17 @@ module spectral_utilities
include 'fftw3-mpi.f03'
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
integer(pInt), public, parameter :: maxPhaseFields = 2_pInt
integer(pInt), public :: nActiveFields = 0_pInt
integer, public, parameter :: maxPhaseFields = 2
integer, public :: nActiveFields = 0
!--------------------------------------------------------------------------------------------------
! field labels information
enum, bind(c)
enumerator :: FIELD_UNDEFINED_ID, &
enumerator :: &
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID, &
FIELD_THERMAL_ID, &
FIELD_DAMAGE_ID, &
FIELD_VACANCYDIFFUSION_ID
FIELD_DAMAGE_ID
end enum
!--------------------------------------------------------------------------------------------------
@ -37,7 +37,7 @@ module spectral_utilities
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
integer(pInt), public :: grid1Red !< grid(1)/2
integer, public :: grid1Red !< grid(1)/2
real (C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_real !< real representation (some stress or deformation) of field_fourier
complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:,:), pointer :: tensorField_fourier !< field on which the Fourier transform operates
real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field real representation for fftw
@ -70,15 +70,17 @@ module spectral_utilities
!--------------------------------------------------------------------------------------------------
! derived types
type, public :: tSolutionState !< return type of solution from spectral solver variants
logical :: converged = .true.
logical :: stagConverged = .true.
logical :: termIll = .false.
integer(pInt) :: iterationsNeeded = 0_pInt
integer :: &
iterationsNeeded = 0
logical :: &
converged = .true., &
stagConverged = .true., &
termIll = .false.
end type tSolutionState
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal
real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal
real(pReal), dimension(3,3) :: values = 0.0_pReal, &
maskFloat = 0.0_pReal
logical, dimension(3,3) :: maskLogical = .false.
character(len=64) :: myType = 'None'
end type tBoundaryCondition
@ -88,10 +90,10 @@ module spectral_utilities
type(tBoundaryCondition) :: stress, & !< stress BC
deformation !< deformation BC (Fdot or L)
real(pReal) :: time = 0.0_pReal !< length of increment
integer(pInt) :: incs = 0_pInt, & !< number of increments
outputfrequency = 1_pInt, & !< frequency of result writes
restartfrequency = 0_pInt, & !< frequency of restart writes
logscale = 0_pInt !< linear/logarithmic time inc flag
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
restartfrequency = 0, & !< frequency of restart writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
end type tLoadCase
@ -102,11 +104,29 @@ module spectral_utilities
real(pReal) :: timeincOld
end type tSolutionParams
type, private :: tNumerics
real(pReal) :: &
FFTW_timelimit !< timelimit for FFTW plan creation, see www.fftw.org
integer :: &
divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints]
logical :: &
memory_efficient !< calculate gamma operator on the fly
character(len=pStringLen) :: &
spectral_derivative, & !< approximation used for derivatives in Fourier space
FFTW_plan_mode, & !< FFTW plan mode, see www.fftw.org
PETSc_defaultOptions, &
PETSc_options
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
enum, bind(c)
enumerator :: DERIVATIVE_CONTINUOUS_ID, &
enumerator :: &
DERIVATIVE_CONTINUOUS_ID, &
DERIVATIVE_CENTRAL_DIFF_ID, &
DERIVATIVE_FWBW_DIFF_ID
end enum
integer(kind(DERIVATIVE_CONTINUOUS_ID)) :: &
spectral_derivative_ID
@ -149,18 +169,14 @@ contains
!> level chosen.
!> Initializes FFTW.
!--------------------------------------------------------------------------------------------------
subroutine utilities_init()
subroutine utilities_init
use IO, only: &
IO_error, &
IO_warning
IO_warning, &
IO_lc
use numerics, only: &
spectral_derivative, &
fftw_planner_flag, &
fftw_timelimit, &
memory_efficient, &
petsc_defaultOptions, &
petsc_options, &
divergence_correction
petsc_options
use debug, only: &
debug_level, &
debug_SPECTRAL, &
@ -169,6 +185,8 @@ subroutine utilities_init()
debug_SPECTRALFFTW, &
debug_SPECTRALPETSC, &
debug_SPECTRALROTATION
use config, only: &
config_numerics
use debug, only: &
PETSCDEBUG
use math
@ -180,8 +198,9 @@ subroutine utilities_init()
implicit none
PetscErrorCode :: ierr
integer(pInt) :: i, j, k
integer(pInt), dimension(3) :: k_s
integer :: i, j, k, &
FFTW_planner_flag
integer, dimension(3) :: k_s
type(C_PTR) :: &
tensorField, & !< field containing data for FFTW in real and fourier space (in place)
vectorField, & !< field containing data for FFTW in real space when debugging FFTW (no in place)
@ -195,6 +214,9 @@ subroutine utilities_init()
write(6,'(/,a)') ' <<<+- spectral_utilities init -+>>>'
write(6,'(/,a)') ' Diehl, Diploma Thesis TU München, 2010'
write(6,'(a)') ' https://doi.org/10.13140/2.1.3234.3840'
write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2012.09.012'
@ -230,7 +252,16 @@ subroutine utilities_init()
write(6,'(/,a,3(i12 ))') ' grid a b c: ', grid
write(6,'(a,3(es12.5))') ' size x y z: ', geomSize
select case (spectral_derivative)
num%memory_efficient = config_numerics%getInt ('memory_efficient', defaultVal=1) > 0
num%FFTW_timelimit = config_numerics%getFloat ('fftw_timelimit', defaultVal=-1.0)
num%divergence_correction = config_numerics%getInt ('divergence_correction', defaultVal=2)
num%spectral_derivative = config_numerics%getString('spectral_derivative', defaultVal='continuous')
num%FFTW_plan_mode = config_numerics%getString('fftw_plan_mode', defaultVal='FFTW_PATIENT')
if (num%divergence_correction < 0 .or. num%divergence_correction > 2) &
call IO_error(301,ext_msg='divergence_correction')
select case (num%spectral_derivative)
case ('continuous')
spectral_derivative_ID = DERIVATIVE_CONTINUOUS_ID
case ('central_difference')
@ -238,21 +269,21 @@ subroutine utilities_init()
case ('fwbw_difference')
spectral_derivative_ID = DERIVATIVE_FWBW_DIFF_ID
case default
call IO_error(892_pInt,ext_msg=trim(spectral_derivative))
call IO_error(892,ext_msg=trim(num%spectral_derivative))
end select
!--------------------------------------------------------------------------------------------------
! scale dimension to calculate either uncorrected, dimension-independent, or dimension- and
! resolution-independent divergence
if (divergence_correction == 1_pInt) then
do j = 1_pInt, 3_pInt
if (num%divergence_correction == 1) then
do j = 1, 3
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 /= int(minloc(geomSize/real(grid,pReal),1),pInt) &
.and. j /= int(maxloc(geomSize/real(grid,pReal),1),pInt)) &
elseif (num%divergence_correction == 2) then
do j = 1, 3
if ( j /= int(minloc(geomSize/real(grid,pReal),1)) &
.and. j /= int(maxloc(geomSize/real(grid,pReal),1))) &
scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal)
enddo
else
@ -260,6 +291,27 @@ subroutine utilities_init()
endif
select case(IO_lc(num%FFTW_plan_mode)) ! setting parameters for the plan creation of FFTW. Basically a translation from fftw3.f
case('estimate','fftw_estimate') ! ordered from slow execution (but fast plan creation) to fast execution
FFTW_planner_flag = 64
case('measure','fftw_measure')
FFTW_planner_flag = 0
case('patient','fftw_patient')
FFTW_planner_flag= 32
case('exhaustive','fftw_exhaustive')
FFTW_planner_flag = 8
case default
call IO_warning(warning_ID=47,ext_msg=trim(IO_lc(num%FFTW_plan_mode)))
FFTW_planner_flag = 32
end select
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) call IO_error(0,ext_msg='Fortran to C') ! check for correct precision in C
call fftw_set_timelimit(num%FFTW_timelimit) ! set timelimit for plan creation
if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6)
!--------------------------------------------------------------------------------------------------
! MPI allocation
gridFFTW = int(grid,C_INTPTR_T)
@ -291,12 +343,12 @@ subroutine utilities_init()
planTensorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
tensorField_real, tensorField_fourier, & ! input data, output data
PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planTensorForth)) call IO_error(810, ext_msg='planTensorForth')
planTensorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
tensorSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
tensorField_fourier,tensorField_real, & ! input data, output data
PETSC_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision
PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
if (.not. C_ASSOCIATED(planTensorBack)) call IO_error(810, ext_msg='planTensorBack')
!--------------------------------------------------------------------------------------------------
@ -304,12 +356,12 @@ subroutine utilities_init()
planVectorForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,&! no. of transforms, default iblock and oblock
vectorField_real, vectorField_fourier, & ! input data, output data
PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planVectorForth)) call IO_error(810, ext_msg='planVectorForth')
planVectorBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
vecSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, & ! no. of transforms, default iblock and oblock
vectorField_fourier,vectorField_real, & ! input data, output data
PETSC_COMM_WORLD, fftw_planner_flag) ! all processors, planer precision
PETSC_COMM_WORLD, FFTW_planner_flag) ! all processors, planer precision
if (.not. C_ASSOCIATED(planVectorBack)) call IO_error(810, ext_msg='planVectorBack')
!--------------------------------------------------------------------------------------------------
@ -317,32 +369,25 @@ subroutine utilities_init()
planScalarForth = fftw_mpi_plan_many_dft_r2c(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
scalarField_real, scalarField_fourier, & ! input data, output data
PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planScalarForth)) call IO_error(810, ext_msg='planScalarForth')
planScalarBack = fftw_mpi_plan_many_dft_c2r(3, [gridFFTW(3),gridFFTW(2),gridFFTW(1)], & ! dimension, logical length in each dimension in reversed order, no. of transforms
scalarSize, FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, &! no. of transforms, default iblock and oblock
scalarField_fourier,scalarField_real, & ! input data, output data
PETSC_COMM_WORLD, fftw_planner_flag) ! use all processors, planer precision
PETSC_COMM_WORLD, FFTW_planner_flag) ! use all processors, planer precision
if (.not. C_ASSOCIATED(planScalarBack)) call IO_error(810, ext_msg='planScalarBack')
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(0_pInt,ext_msg='Fortran to C') ! check for correct precision in C
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
if (debugGeneral) write(6,'(/,a)') ' FFTW initialized'; flush(6)
!--------------------------------------------------------------------------------------------------
! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
do k = grid3Offset+1_pInt, grid3Offset+grid3
k_s(3) = k - 1_pInt
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 > 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
xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s) ! if divergence_correction is set, frequencies are calculated on unit length
do k = grid3Offset+1, grid3Offset+grid3
k_s(3) = k - 1
if(k > grid(3)/2 + 1) 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, grid(2)
k_s(2) = j - 1
if(j > grid(2)/2 + 1) 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, grid1Red
k_s(1) = i - 1 ! symmetry, junst running from 0,1,...,N/2,N/2+1
xi2nd(1:3,i,j,k-grid3Offset) = utilities_getFreqDerivative(k_s)
where(mod(grid,2)==0 .and. [i,j,k] == grid/2+1 .and. &
spectral_derivative_ID == DERIVATIVE_CONTINUOUS_ID) ! for even grids, set the Nyquist Freq component to 0.0
xi1st(1:3,i,j,k-grid3Offset) = cmplx(0.0_pReal,0.0_pReal,pReal)
@ -351,7 +396,7 @@ subroutine utilities_init()
endwhere
enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor
if(num%memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(3,3,3,3,1,1,1), source = cmplx(0.0_pReal,0.0_pReal,pReal))
else ! precalculation of gamma_hat field
allocate (gamma_hat(3,3,3,3,grid1Red,grid(2),grid3), source = cmplx(0.0_pReal,0.0_pReal,pReal))
@ -360,18 +405,17 @@ subroutine utilities_init()
end subroutine utilities_init
!--------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
!> @brief updates reference stiffness and potentially precalculated gamma operator
!> @details Sets the current reference stiffness to the stiffness given as an argument.
!> If the gamma operator is precalculated, it is calculated with this stiffness.
!> In case of an on-the-fly calculation, only the reference stiffness is updated.
!> Also writes out the current reference stiffness for restart.
!--------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------
subroutine utilities_updateGamma(C,saveReference)
use IO, only: &
IO_open_jobFile_binary
use numerics, only: &
memory_efficient, &
worldrank
use mesh, only: &
grid3Offset, &
@ -393,16 +437,14 @@ subroutine utilities_updateGamma(C,saveReference)
logical :: err
C_ref = C
if (saveReference) then
if (worldrank == 0_pInt) then
if (saveReference .and. worldrank == 0) then
write(6,'(/,a)') ' writing reference stiffness to file'
flush(6)
fileUnit = IO_open_jobFile_binary('C_ref','w')
write(fileUnit) C_ref; close(fileUnit)
endif
endif
if(.not. memory_efficient) then
if(.not. num%memory_efficient) then
gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A
do k = grid3Offset+1, grid3Offset+grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
@ -430,11 +472,9 @@ end subroutine utilities_updateGamma
!> @brief forward FFT of data in field_real to field_fourier
!> @details Does an unweighted filtered FFT transform from real to complex
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorForward()
subroutine utilities_FFTtensorForward
implicit none
!--------------------------------------------------------------------------------------------------
! doing the tensor FFT
call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier)
end subroutine utilities_FFTtensorForward
@ -444,7 +484,7 @@ end subroutine utilities_FFTtensorForward
!> @brief backward FFT of data in field_fourier to field_real
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTtensorBackward()
subroutine utilities_FFTtensorBackward
implicit none
call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real)
@ -456,11 +496,9 @@ end subroutine utilities_FFTtensorBackward
!> @brief forward FFT of data in scalarField_real to scalarField_fourier
!> @details Does an unweighted filtered FFT transform from real to complex
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarForward()
subroutine utilities_FFTscalarForward
implicit none
!--------------------------------------------------------------------------------------------------
! doing the scalar FFT
call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier)
end subroutine utilities_FFTscalarForward
@ -470,7 +508,7 @@ end subroutine utilities_FFTscalarForward
!> @brief backward FFT of data in scalarField_fourier to scalarField_real
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTscalarBackward()
subroutine utilities_FFTscalarBackward
implicit none
call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real)
@ -483,11 +521,9 @@ end subroutine utilities_FFTscalarBackward
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
!> @details Does an unweighted filtered FFT transform from real to complex.
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorForward()
subroutine utilities_FFTvectorForward
implicit none
!--------------------------------------------------------------------------------------------------
! doing the vector FFT
call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier)
end subroutine utilities_FFTvectorForward
@ -497,7 +533,7 @@ end subroutine utilities_FFTvectorForward
!> @brief backward FFT of data in field_fourier to field_real
!> @details Does an weighted inverse FFT transform from complex to real
!--------------------------------------------------------------------------------------------------
subroutine utilities_FFTvectorBackward()
subroutine utilities_FFTvectorBackward
implicit none
call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real)
@ -510,8 +546,6 @@ end subroutine utilities_FFTvectorBackward
!> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim
!--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGammaConvolution(fieldAim)
use numerics, only: &
memory_efficient
use math, only: &
math_det33, &
math_invert2
@ -536,7 +570,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation (mechanical equilibrium)
memoryEfficient: if(memory_efficient) then
memoryEfficient: if(num%memory_efficient) then
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any([i,j,k+grid3Offset] /= 1)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
forall(l = 1:3, m = 1:3) &
@ -586,14 +620,14 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
real(pReal), dimension(3,3), intent(in) :: D_ref
real(pReal), intent(in) :: mobility_ref, deltaT
complex(pReal) :: GreenOp_hat
integer(pInt) :: i, j, k
integer :: i, j, k
!--------------------------------------------------------------------------------------------------
! do the actual spectral method calculation
do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red
do k = 1, grid3; do j = 1, grid(2) ;do i = 1, grid1Red
GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal)/ &
(cmplx(mobility_ref,0.0_pReal,pReal) + cmplx(deltaT,0.0_pReal)*&
sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) ! why not use dot_product
sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k))))
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat
enddo; enddo; enddo
@ -612,7 +646,7 @@ real(pReal) function utilities_divergenceRMS()
grid3
implicit none
integer(pInt) :: i, j, k, ierr
integer :: i, j, k, ierr
complex(pReal), dimension(3) :: rescaledGeom
write(6,'(/,a)') ' ... calculating divergence ................................................'
@ -623,8 +657,8 @@ real(pReal) function utilities_divergenceRMS()
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2)
do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
do k = 1, grid3; do j = 1, grid(2)
do i = 2, grid1Red -1 ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)& ! --> sum squared L_2 norm of vector
@ -641,9 +675,9 @@ real(pReal) function utilities_divergenceRMS()
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal)
enddo; enddo
if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
if(grid(1) == 1) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_divergenceRMS')
if(ierr /=0) call IO_error(894, ext_msg='utilities_divergenceRMS')
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
@ -662,7 +696,7 @@ real(pReal) function utilities_curlRMS()
grid3
implicit none
integer(pInt) :: i, j, k, l, ierr
integer :: i, j, k, l, ierr
complex(pReal), dimension(3,3) :: curl_fourier
complex(pReal), dimension(3) :: rescaledGeom
@ -675,9 +709,9 @@ real(pReal) function utilities_curlRMS()
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2);
do i = 2_pInt, grid1Red - 1_pInt
do l = 1_pInt, 3_pInt
do k = 1, grid3; do j = 1, grid(2);
do i = 2, grid1Red - 1
do l = 1, 3
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3))
curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3) &
@ -688,7 +722,7 @@ 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)! Has somewhere a conj. complex counterpart. Therefore count it twice.
enddo
do l = 1_pInt, 3_pInt
do l = 1, 3
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3))
curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3) &
@ -698,7 +732,7 @@ real(pReal) function utilities_curlRMS()
enddo
utilities_curlRMS = utilities_curlRMS &
+ sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal) ! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
do l = 1_pInt, 3_pInt
do l = 1, 3
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) &
@ -711,9 +745,9 @@ real(pReal) function utilities_curlRMS()
enddo; enddo
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_curlRMS')
if(ierr /=0) call IO_error(894, ext_msg='utilities_curlRMS')
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
if(grid(1) == 1) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
end function utilities_curlRMS
@ -738,10 +772,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
real(pReal), intent(in) , dimension(3,3,3,3) :: C !< current average stiffness
real(pReal), intent(in) , dimension(3,3) :: rot_BC !< rotation of load frame
logical, intent(in), dimension(3,3) :: mask_stress !< mask of stress BC
integer(pInt) :: j, k, m, n
integer :: j, k, m, n
logical, dimension(9) :: mask_stressVector
real(pReal), dimension(9,9) :: temp99_Real
integer(pInt) :: size_reduced = 0_pInt
integer :: size_reduced = 0
real(pReal), dimension(:,:), allocatable :: &
s_reduced, & !< reduced compliance matrix (depending on number of stress BC)
c_reduced, & !< reduced stiffness (depending on number of stress BC)
@ -750,8 +784,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
character(len=1024):: formatString
mask_stressVector = reshape(transpose(mask_stress), [9])
size_reduced = int(count(mask_stressVector), pInt)
if(size_reduced > 0_pInt )then
size_reduced = count(mask_stressVector)
if(size_reduced > 0 )then
allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal)
@ -763,37 +797,37 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
transpose(temp99_Real)*1.0e-9_pReal
flush(6)
endif
k = 0_pInt ! calculate reduced stiffness
do n = 1_pInt,9_pInt
k = 0 ! calculate reduced stiffness
do n = 1,9
if(mask_stressVector(n)) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt,9_pInt
k = k + 1
j = 0
do m = 1,9
if(mask_stressVector(m)) then
j = j + 1_pInt
j = j + 1
c_reduced(k,j) = temp99_Real(n,m)
endif; enddo; endif; enddo
call math_invert2(s_reduced, errmatinv, c_reduced) ! invert reduced stiffness
if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
if (errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
if (errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
temp99_Real = 0.0_pReal ! fill up compliance with zeros
k = 0_pInt
do n = 1_pInt,9_pInt
k = 0
do n = 1,9
if(mask_stressVector(n)) then
k = k + 1_pInt
j = 0_pInt
do m = 1_pInt,9_pInt
k = k + 1
j = 0
do m = 1,9
if(mask_stressVector(m)) then
j = j + 1_pInt
j = j + 1
temp99_Real(n,m) = s_reduced(k,j)
endif; enddo; endif; enddo
!--------------------------------------------------------------------------------------------------
! check if inversion was successful
sTimesC = matmul(c_reduced,s_reduced)
do m=1_pInt, size_reduced
do n=1_pInt, size_reduced
do m=1, size_reduced
do n=1, size_reduced
errmatinv = errmatinv &
.or. (m==n .and. abs(sTimesC(m,n)-1.0_pReal) > 1.0e-12_pReal) & ! diagonal elements of S*C should be 1
.or. (m/=n .and. abs(sTimesC(m,n)) > 1.0e-12_pReal) ! off-diagonal elements of S*C should be 0
@ -805,7 +839,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
write(6,trim(formatString),advance='no') ' C * S (load) ', &
transpose(matmul(c_reduced,s_reduced))
write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
endif
else
temp99_real = 0.0_pReal
@ -829,10 +863,10 @@ subroutine utilities_fourierScalarGradient()
grid
implicit none
integer(pInt) :: i, j, k
integer :: i, j, k
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) &
forall(k = 1:grid3, j = 1:grid(2), i = 1:grid1Red) &
vectorField_fourier(1:3,i,j,k) = scalarField_fourier(i,j,k)*xi1st(1:3,i,j,k)
end subroutine utilities_fourierScalarGradient
@ -847,12 +881,12 @@ subroutine utilities_fourierVectorDivergence()
grid
implicit none
integer(pInt) :: i, j, k
integer :: i, j, k
scalarField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
forall(k = 1_pInt:grid3, j = 1_pInt:grid(2), i = 1_pInt:grid1Red) &
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k) + &
sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
forall(k = 1:grid3, j = 1:grid(2), i = 1:grid1Red) &
scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k) &
+ sum(vectorField_fourier(1:3,i,j,k)*conjg(-xi1st(1:3,i,j,k)))
end subroutine utilities_fourierVectorDivergence
@ -866,11 +900,11 @@ subroutine utilities_fourierVectorGradient()
grid
implicit none
integer(pInt) :: i, j, k, m, n
integer :: i, j, k, m, n
tensorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red
do m = 1_pInt, 3_pInt; do n = 1_pInt, 3_pInt
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do m = 1, 3; do n = 1, 3
tensorField_fourier(m,n,i,j,k) = vectorField_fourier(m,i,j,k)*xi1st(n,i,j,k)
enddo; enddo
enddo; enddo; enddo
@ -887,14 +921,13 @@ subroutine utilities_fourierTensorDivergence()
grid
implicit none
integer(pInt) :: i, j, k, m, n
integer :: i, j, k, m, n
vectorField_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid1Red
do m = 1_pInt, 3_pInt; do n = 1_pInt, 3_pInt
vectorField_fourier(m,i,j,k) = &
vectorField_fourier(m,i,j,k) + &
tensorField_fourier(m,n,i,j,k)*conjg(-xi1st(n,i,j,k))
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid1Red
do m = 1, 3; do n = 1, 3
vectorField_fourier(m,i,j,k) = vectorField_fourier(m,i,j,k) &
+ tensorField_fourier(m,n,i,j,k)*conjg(-xi1st(n,i,j,k))
enddo; enddo
enddo; enddo; enddo
@ -934,7 +967,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
integer(pInt) :: &
integer :: &
i,ierr
real(pReal), dimension(3,3,3,3) :: dPdF_max, dPdF_min
real(pReal) :: dPdF_norm_max, dPdF_norm_min
@ -963,7 +996,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_norm_max = 0.0_pReal
dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal)
do i = 1_pInt, product(grid(1:2))*grid3
do i = 1, product(grid(1:2))*grid3
if (dPdF_norm_max < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_max = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_max = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
@ -976,15 +1009,15 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
valueAndRank = [dPdF_norm_max,real(worldrank,pReal)]
call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1, MPI_2DOUBLE_PRECISION, MPI_MAXLOC, PETSC_COMM_WORLD, ierr)
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
if (ierr /= 0) call IO_error(894, ext_msg='MPI_Allreduce max')
call MPI_Bcast(dPdF_max,81,MPI_DOUBLE,int(valueAndRank(2)),PETSC_COMM_WORLD, ierr)
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Bcast max')
if (ierr /= 0) call IO_error(894, ext_msg='MPI_Bcast max')
valueAndRank = [dPdF_norm_min,real(worldrank,pReal)]
call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1, MPI_2DOUBLE_PRECISION, MPI_MINLOC, PETSC_COMM_WORLD, ierr)
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
if (ierr /= 0) call IO_error(894, ext_msg='MPI_Allreduce min')
call MPI_Bcast(dPdF_min,81,MPI_DOUBLE,int(valueAndRank(2)),PETSC_COMM_WORLD, ierr)
if (ierr /= 0_pInt) call IO_error(894_pInt, ext_msg='MPI_Bcast min')
if (ierr /= 0) call IO_error(894, ext_msg='MPI_Bcast min')
C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min)
@ -1074,7 +1107,7 @@ pure function utilities_getFreqDerivative(k_s)
grid
implicit none
integer(pInt), intent(in), dimension(3) :: k_s !< indices of frequency
integer, intent(in), dimension(3) :: k_s !< indices of frequency
complex(pReal), dimension(3) :: utilities_getFreqDerivative
select case (spectral_derivative_ID)
@ -1136,7 +1169,7 @@ subroutine utilities_updateIPcoords(F)
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
integer(pInt) :: i, j, k, m, ierr
integer :: i, j, k, m, ierr
real(pReal), dimension(3) :: step, offset_coords
real(pReal), dimension(3,3) :: Favg
@ -1147,7 +1180,7 @@ subroutine utilities_updateIPcoords(F)
call utilities_FFTtensorForward()
call utilities_fourierTensorDivergence()
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid1Red
if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0,0.0,pReal)))) &
vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ &
sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k))
@ -1157,23 +1190,23 @@ subroutine utilities_updateIPcoords(F)
!--------------------------------------------------------------------------------------------------
! average F
if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
if (grid3Offset == 0) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords')
!--------------------------------------------------------------------------------------------------
! add average to fluctuation and put (0,0,0) on (0,0,0)
step = geomSize/real(grid, pReal)
if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1)
if (grid3Offset == 0) offset_coords = vectorField_real(1:3,1,1,1)
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords')
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
m = 1_pInt
do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)
m = 1
do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) &
+ offset_coords &
+ math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1_pInt,pReal))
m = m+1_pInt
+ math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
m = m+1
enddo; enddo; enddo
end subroutine utilities_updateIPcoords