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

no pInt + same indentation everywhere

This commit is contained in:
Martin Diehl 2019-03-12 11:36:18 +01:00
parent 818338ca93
commit fdc8a848a5
3 changed files with 634 additions and 641 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
src/grid_damage_spectral.f90
View File

@ -1,4 +1,5 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH !> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Spectral solver for nonlocal damage !> @brief Spectral solver for nonlocal damage
@ -16,7 +17,6 @@ module grid_damage_spectral
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
@ -51,8 +51,6 @@ contains
! ToDo: Restart not implemented ! ToDo: Restart not implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_init subroutine grid_damage_spectral_init
use IO, only: &
IO_intOut
use spectral_utilities, only: & use spectral_utilities, only: &
wgt wgt
use mesh, only: & use mesh, only: &
@ -94,22 +92,21 @@ subroutine grid_damage_spectral_init
localK(worldrank+1) = grid3 localK(worldrank+1) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3D(PETSC_COMM_WORLD, & call DMDACreate3D(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & !< global grid grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, & 1, 1, worldsize, &
1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap) 1, 0, & ! #dof (damage phase field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & !< local grid [grid(1)],[grid(2)],localK, & ! local grid
damage_grid,ierr) !< handle, error damage_grid,ierr) ! handle, error
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr) call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr)
call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor) call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,& call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
if (trim(snes_type) == 'vinewtonrsls' .or. & if (trim(snes_type) == 'vinewtonrsls' .or. &
trim(snes_type) == 'vinewtonssls') then trim(snes_type) == 'vinewtonssls') then
@ -117,7 +114,7 @@ subroutine grid_damage_spectral_init
call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr)
call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr)
call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) ! variable bounds for variational inequalities like contact mechanics, damage etc.
call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
endif endif
@ -129,10 +126,10 @@ subroutine grid_damage_spectral_init
xend = xstart + xend - 1 xend = xstart + xend - 1
yend = ystart + yend - 1 yend = ystart + yend - 1
zend = zstart + zend - 1 zend = zstart + zend - 1
call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! damage reference diffusion update ! damage reference diffusion update
@ -151,6 +148,7 @@ subroutine grid_damage_spectral_init
end subroutine grid_damage_spectral_init end subroutine grid_damage_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations !> @brief solution for the spectral damage scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -170,13 +168,12 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case loadCaseTime !< remaining time of current load case
integer :: i, j, k, cell integer :: i, j, k, cell
type(tSolutionState) :: solution
PetscInt ::position PetscInt ::position
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
PetscErrorCode :: ierr PetscErrorCode :: ierr
type(tSolutionState) :: &
solution
SNESConvergedReason :: reason SNESConvergedReason :: reason
solution%converged =.false. solution%converged =.false.
@ -201,7 +198,7 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
damage_stagInc = damage_current damage_stagInc = damage_current
solution%stagConverged = stagNorm < min(err_damage_tolAbs,err_damage_tolRel*solnNorm) solution%stagConverged = stagNorm < min(err_damage_tolAbs, err_damage_tolRel*solnNorm)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating damage state ! updating damage state
@ -226,7 +223,7 @@ end function grid_damage_spectral_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief spectral damage forwarding routine !> @brief spectral damage forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_damage_spectral_forward() subroutine grid_damage_spectral_forward
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
@ -324,18 +321,18 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
! evaluate polarization field ! evaluate polarization field
scalarField_real = 0.0_pReal scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current
call utilities_FFTscalarForward() call utilities_FFTscalarForward
call utilities_fourierScalarGradient() !< calculate gradient of damage field call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward() call utilities_FFTvectorBackward
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward() call utilities_FFTvectorForward
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
@ -350,9 +347,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator ! convolution of damage field with green operator
call utilities_FFTscalarForward() call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > damage_lastInc) &
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_lastInc
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < residualStiffness) &

48
src/grid_mech_spectral_basic.f90
View File

@ -10,7 +10,6 @@ module grid_mech_spectral_basic
use PETScdmda use PETScdmda
use PETScsnes use PETScsnes
use prec, only: & use prec, only: &
pInt, &
pReal pReal
use math, only: & use math, only: &
math_I3 math_I3
@ -59,8 +58,8 @@ module grid_mech_spectral_basic
err_BC, & !< deviation from stress BC err_BC, & !< deviation from stress BC
err_div !< RMS of div of P err_div !< RMS of div of P
integer(pInt), private :: & integer, private :: &
totalIter = 0_pInt !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_spectral_basic_init, & grid_mech_spectral_basic_init, &
@ -130,7 +129,7 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
allocate (F_lastInc (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)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -153,9 +152,9 @@ subroutine grid_mech_spectral_basic_init
call DMsetFromOptions(da,ierr); CHKERRQ(ierr) call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(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 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,grid_mech_spectral_basic_formResidual,PETSC_NULL_SNES,ierr)! residual vector of same shape as solution vector
CHKERRQ(ierr) 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,grid_mech_spectral_basic_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "_converged"
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
@ -163,7 +162,7 @@ subroutine grid_mech_spectral_basic_init
! init fields ! 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) ! get the data out of PETSc to work with
restart: if (restartInc > 0_pInt) then restart: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading values of increment ', restartInc, ' from file' 'reading values of increment ', restartInc, ' from file'
@ -182,11 +181,11 @@ subroutine grid_mech_spectral_basic_init
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F 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) 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') 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 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) 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') if(ierr /=0) call IO_error(894, 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_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 = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restart endif restart
@ -200,7 +199,7 @@ subroutine grid_mech_spectral_basic_init
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
! QUESTION: why not writing back right after reading (l.189)? ! QUESTION: why not writing back right after reading (l.189)?
restartRead: if (restartInc > 0_pInt) then restartRead: if (restartInc > 0) then
if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) &
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & write(6,'(/,a,'//IO_intOut(restartInc)//',a)') &
'reading more values of increment ', restartInc, ' from file' 'reading more values of increment ', restartInc, ' from file'
@ -287,11 +286,8 @@ end function grid_mech_spectral_basic_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector !> @brief forms the basic residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_formResidual(in, & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work) subroutine grid_mech_spectral_basic_formResidual(in, F, &
F, & ! defgrad field on grid residuum, dummy, ierr)
residuum, & ! residuum field on grid
dummy, &
ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin itmin
@ -318,11 +314,11 @@ subroutine grid_mech_spectral_basic_formResidual(in, &
terminallyIll terminallyIll
implicit none implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, & PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
dimension(3,3, XG_RANGE,YG_RANGE,ZG_RANGE), intent(in) :: F intent(in) :: F !< deformation gradient field
PetscScalar, & PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
dimension(3,3, X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: residuum intent(out) :: residuum !< residuum field
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
deltaF_aim deltaF_aim
PetscInt :: & PetscInt :: &
@ -334,11 +330,11 @@ subroutine grid_mech_spectral_basic_formResidual(in, &
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1_pInt totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
@ -366,10 +362,10 @@ subroutine grid_mech_spectral_basic_formResidual(in, &
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward() ! FFT forward of global "tensorField_real" call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC)) ! convolution of Gamma and tensorField_fourier, with arg
call utilities_FFTtensorBackward() ! FFT backward of global tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
@ -507,7 +503,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
write(fileUnit) F_lastInc; close (fileUnit) write(fileUnit) F_lastInc; close (fileUnit)
endif endif
call CPFEM_age() ! age state and kinematics call CPFEM_age ! age state and kinematics
call utilities_updateIPcoords(F) call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg

30
src/grid_thermal_spectral.f90
View File

@ -1,4 +1,5 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH !> @author Shaokang Zhang, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Spectral solver for thermal conduction !> @brief Spectral solver for thermal conduction
@ -16,7 +17,6 @@ module grid_thermal_spectral
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
@ -70,7 +70,7 @@ subroutine grid_thermal_spectral_init
petsc_options petsc_options
implicit none implicit none
integer, dimension(worldsize) :: localK PetscInt, dimension(worldsize) :: localK
integer :: i, j, k, cell integer :: i, j, k, cell
DM :: thermal_grid DM :: thermal_grid
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
@ -100,9 +100,9 @@ subroutine grid_thermal_spectral_init
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, & 1, 1, worldsize, &
1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap) 1, 0, & ! #dof (thermal phase field), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & !< local grid [grid(1)],[grid(2)],localK, & ! local grid
thermal_grid,ierr) !< handle, error thermal_grid,ierr) ! handle, error
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr) call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr)
@ -152,6 +152,7 @@ subroutine grid_thermal_spectral_init
end subroutine grid_thermal_spectral_init end subroutine grid_thermal_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral thermal scheme with internal iterations !> @brief solution for the spectral thermal scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -167,7 +168,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result
thermal_conduction_putTemperatureAndItsRate thermal_conduction_putTemperatureAndItsRate
implicit none implicit none
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment timeinc_old, & !< increment in time of last increment
@ -229,7 +229,7 @@ end function grid_thermal_spectral_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_thermal_spectral_forward() subroutine grid_thermal_spectral_forward
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
@ -332,18 +332,18 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
! evaluate polarization field ! evaluate polarization field
scalarField_real = 0.0_pReal scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current
call utilities_FFTscalarForward() call utilities_FFTscalarForward
call utilities_fourierScalarGradient() !< calculate gradient of damage field call utilities_fourierScalarGradient !< calculate gradient of damage field
call utilities_FFTvectorBackward() call utilities_FFTvectorBackward
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward() call utilities_FFTvectorForward
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
@ -358,9 +358,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator ! convolution of damage field with green operator
call utilities_FFTscalarForward() call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc) call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual