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

Merge branch '146_grid-indexing+allocation' into 'development' 

avoid duplicated variables + name adjustments

See merge request damask/DAMASK!505
This commit is contained in:
Franz Roters 2022-01-26 14:52:43 +00:00
parent 4f926b9686 68ea33c591
commit e7a999ffb1
7 changed files with 298 additions and 253 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/grid/discretization_grid.f90
View File

@ -27,14 +27,14 @@ module discretization_grid
private
integer, dimension(3), public, protected :: &
grid !< (global) grid
grid !< (global) grid
integer, public, protected :: &
grid3, & !< (local) grid in 3rd direction
grid3, & !< (local) grid in 3rd direction
grid3Offset !< (local) grid offset in 3rd direction
real(pReal), dimension(3), public, protected :: &
geomSize !< (global) physical size
geomSize !< (global) physical size
real(pReal), public, protected :: &
size3, & !< (local) size in 3rd direction
size3, & !< (local) size in 3rd direction
size3offset !< (local) size offset in 3rd direction
public :: &

99
src/grid/grid_damage_spectral.f90
View File

@ -39,9 +39,8 @@ module grid_damage_spectral
type(tSolutionParams) :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES :: damage_snes
SNES :: SNES_damage
Vec :: solution_vec
PetscInt :: xstart, xend, ystart, yend, zstart, zend
real(pReal), dimension(:,:,:), allocatable :: &
phi_current, & !< field of current damage
phi_lastInc, & !< field of previous damage
@ -105,10 +104,18 @@ subroutine grid_damage_spectral_init()
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! init fields
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,err_PETSc); CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(damage_snes,'damage_',err_PETSc);CHKERRQ(err_PETSc)
call SNESCreate(PETSC_COMM_WORLD,SNES_damage,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -122,39 +129,41 @@ subroutine grid_damage_spectral_init()
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
damage_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call SNESSetDM(damage_snes,damage_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
call DMsetFromOptions(damage_grid,err_PETSc); CHKERRQ(err_PETSc)
call DMsetUp(damage_grid,err_PETSc); CHKERRQ(err_PETSc)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMsetFromOptions(damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
call SNESSetFromOptions(damage_snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
call SNESGetType(damage_snes,snes_type,err_PETSc); CHKERRQ(err_PETSc)
call SNESSetDM(SNES_damage,damage_grid,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetFromOptions(SNES_damage,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc)
call SNESGetType(SNES_damage,snes_type,err_PETSc)
CHKERRQ(err_PETSc)
if (trim(snes_type) == 'vinewtonrsls' .or. &
trim(snes_type) == 'vinewtonssls') then
call DMGetGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc)
call DMGetGlobalVector(damage_grid,uBound,err_PETSc); CHKERRQ(err_PETSc)
call VecSet(lBound,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
call VecSet(uBound,1.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
call SNESVISetVariableBounds(damage_snes,lBound,uBound,err_PETSc) ! variable bounds for variational inequalities like contact mechanics, damage etc.
call DMRestoreGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc)
call DMRestoreGlobalVector(damage_grid,uBound,err_PETSc); CHKERRQ(err_PETSc)
call DMGetGlobalVector(damage_grid,lBound,err_PETSc)
CHKERRQ(err_PETSc)
call DMGetGlobalVector(damage_grid,uBound,err_PETSc)
CHKERRQ(err_PETSc)
call VecSet(lBound,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call VecSet(uBound,1.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call SNESVISetVariableBounds(SNES_damage,lBound,uBound,err_PETSc) ! variable bounds for variational inequalities
CHKERRQ(err_PETSc)
call DMRestoreGlobalVector(damage_grid,lBound,err_PETSc)
CHKERRQ(err_PETSc)
call DMRestoreGlobalVector(damage_grid,uBound,err_PETSc)
CHKERRQ(err_PETSc)
end if
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc)
call VecSet(solution_vec,1.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
call VecSet(solution_vec,1.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
call updateReference
call updateReference()
end subroutine grid_damage_spectral_init
@ -181,9 +190,9 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
! set module wide availabe data
params%Delta_t = Delta_t
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,err_PETSc)
call SNESSolve(SNES_damage,PETSC_NULL_VEC,solution_vec,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetConvergedReason(damage_snes,reason,err_PETSc)
call SNESGetConvergedReason(SNES_damage,reason,err_PETSc)
CHKERRQ(err_PETSc)
if (reason < 1) then
@ -209,8 +218,10 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
call homogenization_set_phi(phi_current(i,j,k),ce)
end do; end do; end do
call VecMin(solution_vec,devNull,phi_min,err_PETSc); CHKERRQ(err_PETSc)
call VecMax(solution_vec,devNull,phi_max,err_PETSc); CHKERRQ(err_PETSc)
call VecMin(solution_vec,devNull,phi_min,err_PETSc)
CHKERRQ(err_PETSc)
call VecMax(solution_vec,devNull,phi_max,err_PETSc)
CHKERRQ(err_PETSc)
if (solution%converged) &
print'(/,1x,a)', '... nonlocal damage converged .....................................'
print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm
@ -228,7 +239,7 @@ subroutine grid_damage_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, ce
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscScalar, dimension(:,:,:), pointer :: phi_PETSc
PetscErrorCode :: err_PETSc
if (cutBack) then
@ -236,11 +247,12 @@ subroutine grid_damage_spectral_forward(cutBack)
phi_stagInc = phi_lastInc
!--------------------------------------------------------------------------------------------------
! reverting damage field state
call SNESGetDM(damage_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,err_PETSc) !< get the data out of PETSc to work with
call SNESGetDM(SNES_damage,dm_local,err_PETSc)
CHKERRQ(err_PETSc)
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc) !< get the data out of PETSc to work with
CHKERRQ(err_PETSc)
phi_PETSc = phi_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
@ -258,7 +270,7 @@ end subroutine grid_damage_spectral_forward
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral damage residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
@ -267,9 +279,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
x_scal
PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
r
PetscObject :: dummy
PetscErrorCode :: dummy_err
PetscErrorCode :: err_PETSc
integer :: i, j, k, ce
@ -310,7 +322,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
r = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
err_PETSc = 0
end subroutine formResidual

186
src/grid/grid_mech_FEM.f90
View File

@ -50,7 +50,7 @@ module grid_mechanical_FEM
!--------------------------------------------------------------------------------------------------
! PETSc data
DM :: mechanical_grid
SNES :: mechanical_snes
SNES :: SNES_mechanical
Vec :: solution_current, solution_lastInc, solution_rate
!--------------------------------------------------------------------------------------------------
@ -60,7 +60,6 @@ module grid_mechanical_FEM
real(pReal), dimension(3) :: delta
real(pReal), dimension(3,8) :: BMat
real(pReal), dimension(8,8) :: HGMat
PetscInt :: xstart,ystart,zstart,xend,yend,zend
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
@ -146,7 +145,7 @@ subroutine grid_mechanical_FEM_init
! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
'-mechanical_snes_type newtonls -mechanical_ksp_type fgmres &
&-mechanical_ksp_max_it 25 -mechanical_mg_levels_ksp_type chebyshev', &
&-mechanical_ksp_max_it 25', &
err_PETSc)
CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
@ -160,9 +159,9 @@ subroutine grid_mechanical_FEM_init
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,err_PETSc)
call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',err_PETSc)
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
@ -177,8 +176,6 @@ subroutine grid_mechanical_FEM_init
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetDM(mechanical_snes,mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetFromOptions(mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(mechanical_grid,err_PETSc)
@ -195,28 +192,28 @@ subroutine grid_mechanical_FEM_init
CHKERRQ(err_PETSc)
call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetConvergenceTest(mechanical_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "_converged"
call SNESSetConvergenceTest(SNES_mechanical,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "_converged"
CHKERRQ(err_PETSc)
call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1_pPetscInt), err_PETSc) ! ignore linear solve failures
call SNESSetMaxLinearSolveFailures(SNES_mechanical, huge(1_pPetscInt), err_PETSc) ! ignore linear solve failures
CHKERRQ(err_PETSc)
call SNESSetFromOptions(mechanical_snes,err_PETSc) ! pull it all together with additional cli arguments
call SNESSetDM(SNES_mechanical,mechanical_grid,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional cli arguments
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! init fields
call VecSet(solution_current,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
call VecSet(solution_lastInc,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
call VecSet(solution_rate ,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
call VecSet(solution_current,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call VecSet(solution_lastInc,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call VecSet(solution_rate ,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAGetCorners(mechanical_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc) ! local grid extent
CHKERRQ(err_PETSc)
xend = xstart+xend-1
yend = ystart+yend-1
zend = zstart+zend-1
delta = geomSize/real(grid,pReal) ! grid spacing
detJ = product(delta) ! cell volume
@ -311,14 +308,9 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
! update stiffness (and gamma operator)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(mechanical_snes,PETSC_NULL_VEC,solution_current,err_PETSc)
call SNESsolve(SNES_mechanical,PETSC_NULL_VEC,solution_current,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc)
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
CHKERRQ(err_PETSc)
solution%converged = reason > 0
@ -386,9 +378,11 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
call VecScale(solution_rate,1.0_pReal/Delta_t_old,err_PETSc)
CHKERRQ(err_PETSc)
else
call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
call VecSet(solution_rate,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
endif
call VecCopy(solution_current,solution_lastInc,err_PETSc); CHKERRQ(err_PETSc)
call VecCopy(solution_current,solution_lastInc,err_PETSc)
CHKERRQ(err_PETSc)
F_lastInc = F
@ -515,6 +509,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,e
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT)
err_PETSc = 0
end subroutine converged
@ -527,7 +522,7 @@ subroutine formResidual(da_local,x_local, &
DM :: da_local
Vec :: x_local, f_local
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, r
PetscScalar, dimension(8,3) :: x_elem, f_elem
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
PetscInt :: &
@ -538,9 +533,9 @@ subroutine formResidual(da_local,x_local, &
integer(MPI_INTEGER_KIND) :: err_MPI
real(pReal), dimension(3,3,3,3) :: devNull
call SNESGetNumberFunctionEvals(mechanical_snes,nfuncs,err_PETSc)
call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetIterationNumber(mechanical_snes,PETScIter,err_PETSc)
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
CHKERRQ(err_PETSc)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -559,17 +554,18 @@ subroutine formResidual(da_local,x_local, &
!--------------------------------------------------------------------------------------------------
! get deformation gradient
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc)
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
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) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -586,47 +582,53 @@ subroutine formResidual(da_local,x_local, &
!--------------------------------------------------------------------------------------------------
! constructing residual
call VecSet(f_local,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc)
call VecSet(f_local,0.0_pReal,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
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 + &
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,i,j,k-grid3offset)))*detJ + &
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
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)
r(0:2,i+ii,j+jj,k+kk) = r(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,err_PETSc);CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
call DMDAVecGetArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
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,err_PETSc);CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
if (grid3offset == 0) then
r(0:2,0, 0, 0) = 0.0_pReal
r(0:2,grid(1),0, 0) = 0.0_pReal
r(0:2,0, grid(2),0) = 0.0_pReal
r(0:2,grid(1),grid(2),0) = 0.0_pReal
end if
if (grid3+grid3offset == grid(3)) then
r(0:2,0, 0, grid(3)) = 0.0_pReal
r(0:2,grid(1),0, grid(3)) = 0.0_pReal
r(0:2,0, grid(2),grid(3)) = 0.0_pReal
r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal
end if
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine formResidual
@ -643,7 +645,7 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
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 :: i, ii, j, jj, k, kk, ctr, ce
PetscInt,dimension(3),parameter :: rows = [0, 1, 2]
PetscScalar :: diag
PetscObject :: dummy
@ -658,11 +660,12 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
CHKERRQ(err_PETSc)
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,err_PETSc)
CHKERRQ(err_PETSc)
call MatZeroEntries(Jac,err_PETSc); CHKERRQ(err_PETSc)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
call MatZeroEntries(Jac,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
do kk = -1, 0; do jj = -1, 0; do ii = -1, 0
ctr = ctr + 1
col(MatStencil_i,ctr ) = i+ii
col(MatStencil_j,ctr ) = j+jj
@ -678,49 +681,52 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
col(MatStencil_c,ctr+16) = 2
enddo; enddo; enddo
row = col
ele = ele + 1
ce = ce + 1
K_ele = 0.0
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
homogenization_dPdF(2,2,2,2,ele) + &
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
homogenization_dPdF(2,2,2,2,ce) + &
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
K_ele = K_ele + &
matmul(transpose(BMatFull), &
matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,ele), &
matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,ce), &
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
call MatSetValuesStencil(Jac,24_pPETScInt,row,24_pPetscInt,col,K_ele,ADD_VALUES,err_PETSc)
CHKERRQ(err_PETSc)
enddo; enddo; enddo
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
CHKERRQ(err_PETSc)
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
CHKERRQ(err_PETSc)
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc)
CHKERRQ(err_PETSc)
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! applying boundary conditions
diag = (C_volAvg(1,1,1,1)/delta(1)**2 + &
C_volAvg(2,2,2,2)/delta(2)**2 + &
C_volAvg(3,3,3,3)/delta(3)**2)*detJ
diag = (C_volAvg(1,1,1,1)/delta(1)**2 + C_volAvg(2,2,2,2)/delta(2)**2 + C_volAvg(3,3,3,3)/delta(3)**2) &
* detJ
call MatZeroRowsColumns(Jac,size(rows,kind=pPetscInt),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,err_PETSc)
CHKERRQ(err_PETSc)
call DMGetGlobalVector(da_local,coordinates,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc)
CHKERRQ(err_PETSc)
ele = 0
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
ele = ele + 1
x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele)
ce = 0
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
ce = ce + 1
x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce)
enddo; enddo; enddo
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,err_PETSc)
CHKERRQ(err_PETSc) ! initialize to undeformed coordinates (ToDo: use ip coordinates)
CHKERRQ(err_PETSc) ! initialize to undeformed coordinates (ToDo: use ip coordinates)
call MatNullSpaceCreateRigidBody(coordinates,matnull,err_PETSc)
CHKERRQ(err_PETSc) ! get rigid body deformation modes
CHKERRQ(err_PETSc) ! get rigid body deformation modes
call DMRestoreGlobalVector(da_local,coordinates,err_PETSc)
CHKERRQ(err_PETSc)
call MatSetNullSpace(Jac,matnull,err_PETSc)

77
src/grid/grid_mech_spectral_basic.f90
View File

@ -50,7 +50,7 @@ module grid_mechanical_spectral_basic
!--------------------------------------------------------------------------------------------------
! PETSc data
DM :: da
SNES :: snes
SNES :: SNES_mechanical
Vec :: solution_vec
!--------------------------------------------------------------------------------------------------
@ -158,8 +158,10 @@ subroutine grid_mechanical_spectral_basic_init
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,err_PETSc); CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(snes,'mechanical_',err_PETSc);CHKERRQ(err_PETSc)
call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -173,19 +175,25 @@ subroutine grid_mechanical_spectral_basic_init
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call SNESSetDM(snes,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
call DMsetFromOptions(da,err_PETSc); CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc); CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor)
call DMsetFromOptions(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged"
call SNESsetConvergenceTest(SNES_mechanical,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged"
CHKERRQ(err_PETSc)
call SNESSetDM(SNES_mechanical,da,err_PETSc)
CHKERRQ(err_PETSc)
call SNESsetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc)
call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) ! places pointer on PETSc data
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc) ! places pointer on PETSc data
CHKERRQ(err_PETSc)
restartRead: if (interface_restartInc > 0) then
print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
@ -218,7 +226,8 @@ subroutine grid_mechanical_spectral_basic_init
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc) ! deassociate pointer
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc) ! deassociate pointer
CHKERRQ(err_PETSc)
restartRead2: if (interface_restartInc > 0) then
print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
@ -270,13 +279,10 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,err_PETSc); CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,err_PETSc); CHKERRQ(err_PETSc)
call SNESsolve(SNES_mechanical,PETSC_NULL_VEC,solution_vec,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
CHKERRQ(err_PETSc)
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
@ -310,7 +316,8 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
PetscScalar, pointer, dimension(:,:,:,:) :: F
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
if (cutBack) then
C_volAvg = C_volAvgLastInc
@ -353,7 +360,8 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! set module wide available data
@ -372,9 +380,11 @@ subroutine grid_mechanical_spectral_basic_updateCoords
PetscErrorCode :: err_PETSc
PetscScalar, dimension(:,:,:,:), pointer :: F
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_updateCoords(F)
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine grid_mechanical_spectral_basic_updateCoords
@ -388,7 +398,8 @@ subroutine grid_mechanical_spectral_basic_restartWrite
integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:,:), pointer :: F
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
@ -415,7 +426,8 @@ subroutine grid_mechanical_spectral_basic_restartWrite
if (num%update_gamma) call utilities_saveReferenceStiffness
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecRestoreArrayF90(da,solution_vec,F,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine grid_mechanical_spectral_basic_restartWrite
@ -457,6 +469,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT)
err_PETSc = 0
end subroutine converged
@ -465,13 +478,13 @@ end subroutine converged
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, F, &
residuum, dummy, err_PETSc)
r, dummy, err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
intent(in) :: F !< deformation gradient field
PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
intent(out) :: residuum !< residuum field
intent(out) :: r !< residuum field
real(pReal), dimension(3,3) :: &
deltaF_aim
PetscInt :: &
@ -481,8 +494,10 @@ subroutine formResidual(in, F, &
PetscErrorCode :: err_PETSc
integer(MPI_INTEGER_KIND) :: err_MPI
call SNESGetNumberFunctionEvals(snes,nfuncs,err_PETSc); CHKERRQ(err_PETSc)
call SNESGetIterationNumber(snes,PETScIter,err_PETSc); CHKERRQ(err_PETSc)
call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
CHKERRQ(err_PETSc)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -500,7 +515,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
call utilities_constitutiveResponse(r, & ! residuum gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
@ -515,7 +530,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
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) = r ! store fPK field for subsequent FFT forward transform
call utilities_FFTtensorForward ! FFT forward of global "tensorField_real"
err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
@ -523,7 +538,7 @@ subroutine 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
r = 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 formResidual

92
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -55,7 +55,7 @@ module grid_mechanical_spectral_polarisation
!--------------------------------------------------------------------------------------------------
! PETSc data
DM :: da
SNES :: snes
SNES :: SNES_mechanical
Vec :: solution_vec
!--------------------------------------------------------------------------------------------------
@ -178,8 +178,10 @@ subroutine grid_mechanical_spectral_polarisation_init
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,err_PETSc); CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(snes,'mechanical_',err_PETSc);CHKERRQ(err_PETSc)
call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -193,19 +195,25 @@ subroutine grid_mechanical_spectral_polarisation_init
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
da,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call SNESSetDM(snes,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
call DMsetFromOptions(da,err_PETSc); CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc); CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMsetFromOptions(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(da,err_PETSc)
CHKERRQ(err_PETSc)
call DMcreateGlobalVector(da,solution_vec,err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged"
call SNESsetConvergenceTest(SNES_mechanical,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,err_PETSc) ! specify custom convergence check function "converged"
CHKERRQ(err_PETSc)
call SNESSetDM(SNES_mechanical,da,err_PETSc)
CHKERRQ(err_PETSc)
call SNESsetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc)
call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc); CHKERRQ(err_PETSc) ! places pointer on PETSc data
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc) ! places pointer on PETSc data
CHKERRQ(err_PETSc)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
@ -285,7 +293,7 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
! input data for solution
character(len=*), intent(in) :: &
incInfoIn
type(tSolutionState) :: &
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
! PETSc Data
@ -303,13 +311,10 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
S_scale = math_invSym3333(C_minMaxAvg)
end if
!--------------------------------------------------------------------------------------------------
! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,err_PETSc); CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! check convergence
call SNESGetConvergedReason(snes,reason,err_PETSc); CHKERRQ(err_PETSc)
call SNESSolve(SNES_mechanical,PETSC_NULL_VEC,solution_vec,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
CHKERRQ(err_PETSc)
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
@ -345,7 +350,8 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
real(pReal), dimension(3,3) :: F_lambda33
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
@ -446,7 +452,8 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc); CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,err_PETSc)
CHKERRQ(err_PETSc)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
@ -523,6 +530,7 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
print'(/,1x,a)', '==========================================================================='
flush(IO_STDOUT)
err_PETSc = 0
end subroutine converged
@ -531,18 +539,18 @@ end subroutine converged
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, &
residuum, dummy,err_PETSc)
r, dummy,err_PETSc)
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
target, intent(out) :: r !< residuum field
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, &
F_tau, &
residual_F, &
residual_F_tau
r_F, &
r_F_tau
PetscInt :: &
PETScIter, &
nfuncs
@ -554,21 +562,23 @@ subroutine formResidual(in, FandF_tau, &
!---------------------------------------------------------------------------------------------------
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 => 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)
r_F => r(1:3,1:3,1,&
X_RANGE, Y_RANGE, Z_RANGE)
r_F_tau => r(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_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
call SNESGetNumberFunctionEvals(snes,nfuncs,err_PETSc); CHKERRQ(err_PETSc)
call SNESGetIterationNumber(snes,PETScIter,err_PETSc); CHKERRQ(err_PETSc)
call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
CHKERRQ(err_PETSc)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -602,13 +612,13 @@ subroutine formResidual(in, FandF_tau, &
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
r_F_tau = num%beta*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)
call utilities_constitutiveResponse(r_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/num%beta,params%Delta_t,params%rotation_BC)
F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
!--------------------------------------------------------------------------------------------------
@ -619,7 +629,7 @@ subroutine formResidual(in, FandF_tau, &
params%stress_mask)))
! 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
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = r_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
@ -628,11 +638,11 @@ subroutine formResidual(in, FandF_tau, &
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) = &
r_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
r_F(1:3,1:3,i,j,k) - matmul(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)
+ r_F_tau(1:3,1:3,i,j,k)
end do; end do; end do
!--------------------------------------------------------------------------------------------------

87
src/grid/grid_thermal_spectral.f90
View File

@ -38,9 +38,8 @@ module grid_thermal_spectral
type(tSolutionParams) :: params
!--------------------------------------------------------------------------------------------------
! PETSc data
SNES :: thermal_snes
Vec :: solution_vec
PetscInt :: xstart, xend, ystart, yend, zstart, zend
SNES :: SNES_thermal
Vec :: solution_vec
real(pReal), dimension(:,:,:), allocatable :: &
T_current, & !< field of current temperature
T_lastInc, & !< field of previous temperature
@ -100,10 +99,24 @@ subroutine grid_thermal_spectral_init(T_0)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc)
!--------------------------------------------------------------------------------------------------
! init fields
allocate(T_current(grid(1),grid(2),grid3), source=T_0)
allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0)
allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
call homogenization_thermal_setField(T_0,0.0_pReal,ce)
end do; end do; end do
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,thermal_snes,err_PETSc); CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(thermal_snes,'thermal_',err_PETSc);CHKERRQ(err_PETSc)
call SNESCreate(PETSC_COMM_WORLD,SNES_thermal,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc)
CHKERRQ(err_PETSc)
localK = 0_pPetscInt
localK(worldrank) = int(grid3,pPetscInt)
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
@ -117,42 +130,25 @@ subroutine grid_thermal_spectral_init(T_0)
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
thermal_grid,err_PETSc) ! handle, error
CHKERRQ(err_PETSc)
call SNESSetDM(thermal_snes,thermal_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
call DMsetFromOptions(thermal_grid,err_PETSc); CHKERRQ(err_PETSc)
call DMsetUp(thermal_grid,err_PETSc); CHKERRQ(err_PETSc)
call DMsetFromOptions(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMsetUp(thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
CHKERRQ(err_PETSc)
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
CHKERRQ(err_PETSc)
call SNESSetFromOptions(thermal_snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc)
call SNESSetDM(SNES_thermal,thermal_grid,err_PETSc)
CHKERRQ(err_PETSc)
call SNESSetFromOptions(SNES_thermal,err_PETSc) ! pull it all together with additional CLI arguments
CHKERRQ(err_PETSc)
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1
T_current(i,j,k) = T_0
T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k)
call homogenization_thermal_setField(T_0,0.0_pReal,ce)
end do; end do; end do
call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current
T_PETSc = T_current
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
call updateReference
call updateReference()
end subroutine grid_thermal_spectral_init
@ -179,9 +175,9 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
! set module wide availabe data
params%Delta_t = Delta_t
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,err_PETSc)
call SNESSolve(SNES_thermal,PETSC_NULL_VEC,solution_vec,err_PETSc)
CHKERRQ(err_PETSc)
call SNESGetConvergedReason(thermal_snes,reason,err_PETSc)
call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc)
CHKERRQ(err_PETSc)
if (reason < 1) then
@ -207,8 +203,10 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
end do; end do; end do
call VecMin(solution_vec,devNull,T_min,err_PETSc); CHKERRQ(err_PETSc)
call VecMax(solution_vec,devNull,T_max,err_PETSc); CHKERRQ(err_PETSc)
call VecMin(solution_vec,devNull,T_min,err_PETSc)
CHKERRQ(err_PETSc)
call VecMax(solution_vec,devNull,T_max,err_PETSc)
CHKERRQ(err_PETSc)
if (solution%converged) &
print'(/,1x,a)', '... thermal conduction converged ..................................'
print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
@ -226,7 +224,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
logical, intent(in) :: cutBack
integer :: i, j, k, ce
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscScalar, dimension(:,:,:), pointer :: T_PETSc
PetscErrorCode :: err_PETSc
if (cutBack) then
@ -235,12 +233,12 @@ subroutine grid_thermal_spectral_forward(cutBack)
!--------------------------------------------------------------------------------------------------
! reverting thermal field state
call SNESGetDM(thermal_snes,dm_local,err_PETSc)
call SNESGetDM(SNES_thermal,dm_local,err_PETSc)
CHKERRQ(err_PETSc)
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,err_PETSc) !< get the data out of PETSc to work with
call DMDAVecGetArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc) !< get the data out of PETSc to work with
CHKERRQ(err_PETSc)
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,err_PETSc)
T_PETSc = T_current
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc)
CHKERRQ(err_PETSc)
ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
@ -258,7 +256,7 @@ end subroutine grid_thermal_spectral_forward
!--------------------------------------------------------------------------------------------------
!> @brief forms the spectral thermal residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
subroutine formResidual(in,x_scal,r,dummy,err_PETSc)
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
@ -267,9 +265,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
x_scal
PetscScalar, dimension( &
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
f_scal
r
PetscObject :: dummy
PetscErrorCode :: dummy_err
PetscErrorCode :: err_PETSc
integer :: i, j, k, ce
T_current = x_scal
@ -304,7 +302,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
r = T_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
err_PETSc = 0
end subroutine formResidual

2
src/parallelization.f90
View File

@ -53,6 +53,7 @@ contains
subroutine parallelization_init
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
character(len=4) :: rank_str
!$ integer :: got_env, threadLevel
!$ integer(pI32) :: OMP_NUM_THREADS
!$ character(len=6) NumThreadsString
@ -112,6 +113,7 @@ subroutine parallelization_init
if (worldrank /= 0) then
close(OUTPUT_UNIT) ! disable output
write(rank_str,'(i4.4)') worldrank ! use for MPI debug filenames
open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd
endif