Merge remote-tracking branch 'origin/development' into None-instead-of-optional
This commit is contained in:
commit
e8d20a8fec
|
@ -27,14 +27,14 @@ module discretization_grid
|
||||||
private
|
private
|
||||||
|
|
||||||
integer, dimension(3), public, protected :: &
|
integer, dimension(3), public, protected :: &
|
||||||
grid !< (global) grid
|
grid !< (global) grid
|
||||||
integer, public, protected :: &
|
integer, public, protected :: &
|
||||||
grid3, & !< (local) grid in 3rd direction
|
grid3, & !< (local) grid in 3rd direction
|
||||||
grid3Offset !< (local) grid offset in 3rd direction
|
grid3Offset !< (local) grid offset in 3rd direction
|
||||||
real(pReal), dimension(3), public, protected :: &
|
real(pReal), dimension(3), public, protected :: &
|
||||||
geomSize !< (global) physical size
|
geomSize !< (global) physical size
|
||||||
real(pReal), public, protected :: &
|
real(pReal), public, protected :: &
|
||||||
size3, & !< (local) size in 3rd direction
|
size3, & !< (local) size in 3rd direction
|
||||||
size3offset !< (local) size offset in 3rd direction
|
size3offset !< (local) size offset in 3rd direction
|
||||||
|
|
||||||
public :: &
|
public :: &
|
||||||
|
|
|
@ -39,9 +39,8 @@ module grid_damage_spectral
|
||||||
type(tSolutionParams) :: params
|
type(tSolutionParams) :: params
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc data
|
! PETSc data
|
||||||
SNES :: damage_snes
|
SNES :: SNES_damage
|
||||||
Vec :: solution_vec
|
Vec :: solution_vec
|
||||||
PetscInt :: xstart, xend, ystart, yend, zstart, zend
|
|
||||||
real(pReal), dimension(:,:,:), allocatable :: &
|
real(pReal), dimension(:,:,:), allocatable :: &
|
||||||
phi_current, & !< field of current damage
|
phi_current, & !< field of current damage
|
||||||
phi_lastInc, & !< field of previous 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)
|
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
|
||||||
CHKERRQ(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
|
! initialize solver specific parts of PETSc
|
||||||
call SNESCreate(PETSC_COMM_WORLD,damage_snes,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESCreate(PETSC_COMM_WORLD,SNES_damage,err_PETSc)
|
||||||
call SNESSetOptionsPrefix(damage_snes,'damage_',err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call SNESSetOptionsPrefix(SNES_damage,'damage_',err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
localK = 0_pPetscInt
|
localK = 0_pPetscInt
|
||||||
localK(worldrank) = int(grid3,pPetscInt)
|
localK(worldrank) = int(grid3,pPetscInt)
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
|
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
|
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
|
||||||
damage_grid,err_PETSc) ! handle, error
|
damage_grid,err_PETSc) ! handle, error
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetDM(damage_snes,damage_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
|
call DMsetFromOptions(damage_grid,err_PETSc)
|
||||||
call DMsetFromOptions(damage_grid,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMsetUp(damage_grid,err_PETSc); CHKERRQ(err_PETSc)
|
call DMsetUp(damage_grid,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)
|
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
|
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetFromOptions(damage_snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
|
call SNESSetDM(SNES_damage,damage_grid,err_PETSc)
|
||||||
call SNESGetType(damage_snes,snes_type,err_PETSc); CHKERRQ(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. &
|
if (trim(snes_type) == 'vinewtonrsls' .or. &
|
||||||
trim(snes_type) == 'vinewtonssls') then
|
trim(snes_type) == 'vinewtonssls') then
|
||||||
call DMGetGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc)
|
call DMGetGlobalVector(damage_grid,lBound,err_PETSc)
|
||||||
call DMGetGlobalVector(damage_grid,uBound,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call VecSet(lBound,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
|
call DMGetGlobalVector(damage_grid,uBound,err_PETSc)
|
||||||
call VecSet(uBound,1.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESVISetVariableBounds(damage_snes,lBound,uBound,err_PETSc) ! variable bounds for variational inequalities like contact mechanics, damage etc.
|
call VecSet(lBound,0.0_pReal,err_PETSc)
|
||||||
call DMRestoreGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMRestoreGlobalVector(damage_grid,uBound,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
|
end if
|
||||||
|
call VecSet(solution_vec,1.0_pReal,err_PETSc)
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
! init fields
|
|
||||||
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc)
|
|
||||||
CHKERRQ(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
|
end subroutine grid_damage_spectral_init
|
||||||
|
|
||||||
|
@ -181,9 +190,9 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
|
||||||
! set module wide availabe data
|
! set module wide availabe data
|
||||||
params%Delta_t = Delta_t
|
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)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESGetConvergedReason(damage_snes,reason,err_PETSc)
|
call SNESGetConvergedReason(SNES_damage,reason,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
if (reason < 1) then
|
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)
|
call homogenization_set_phi(phi_current(i,j,k),ce)
|
||||||
end do; end do; end do
|
end do; end do; end do
|
||||||
|
|
||||||
call VecMin(solution_vec,devNull,phi_min,err_PETSc); CHKERRQ(err_PETSc)
|
call VecMin(solution_vec,devNull,phi_min,err_PETSc)
|
||||||
call VecMax(solution_vec,devNull,phi_max,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call VecMax(solution_vec,devNull,phi_max,err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
if (solution%converged) &
|
if (solution%converged) &
|
||||||
print'(/,1x,a)', '... nonlocal damage 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
|
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
|
logical, intent(in) :: cutBack
|
||||||
integer :: i, j, k, ce
|
integer :: i, j, k, ce
|
||||||
DM :: dm_local
|
DM :: dm_local
|
||||||
PetscScalar, dimension(:,:,:), pointer :: x_scal
|
PetscScalar, dimension(:,:,:), pointer :: phi_PETSc
|
||||||
PetscErrorCode :: err_PETSc
|
PetscErrorCode :: err_PETSc
|
||||||
|
|
||||||
if (cutBack) then
|
if (cutBack) then
|
||||||
|
@ -236,11 +247,12 @@ subroutine grid_damage_spectral_forward(cutBack)
|
||||||
phi_stagInc = phi_lastInc
|
phi_stagInc = phi_lastInc
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! reverting damage field state
|
! reverting damage field state
|
||||||
call SNESGetDM(damage_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESGetDM(SNES_damage,dm_local,err_PETSc)
|
||||||
call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,err_PETSc) !< get the data out of PETSc to work with
|
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
|
call DMDAVecGetArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc) !< get the data out of PETSc to work with
|
||||||
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
phi_PETSc = phi_current
|
||||||
|
call DMDAVecRestoreArrayF90(dm_local,solution_vec,phi_PETSc,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
ce = 0
|
ce = 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)
|
||||||
|
@ -258,7 +270,7 @@ end subroutine grid_damage_spectral_forward
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief forms the spectral damage residual vector
|
!> @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) :: &
|
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
|
||||||
in
|
in
|
||||||
|
@ -267,9 +279,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
|
||||||
x_scal
|
x_scal
|
||||||
PetscScalar, dimension( &
|
PetscScalar, dimension( &
|
||||||
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
||||||
f_scal
|
r
|
||||||
PetscObject :: dummy
|
PetscObject :: dummy
|
||||||
PetscErrorCode :: dummy_err
|
PetscErrorCode :: err_PETSc
|
||||||
integer :: i, j, k, ce
|
integer :: i, j, k, ce
|
||||||
|
|
||||||
|
|
||||||
|
@ -310,7 +322,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! constructing residual
|
! 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
|
end subroutine formResidual
|
||||||
|
|
||||||
|
|
|
@ -50,7 +50,7 @@ module grid_mechanical_FEM
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc data
|
! PETSc data
|
||||||
DM :: mechanical_grid
|
DM :: mechanical_grid
|
||||||
SNES :: mechanical_snes
|
SNES :: SNES_mechanical
|
||||||
Vec :: solution_current, solution_lastInc, solution_rate
|
Vec :: solution_current, solution_lastInc, solution_rate
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -60,7 +60,6 @@ module grid_mechanical_FEM
|
||||||
real(pReal), dimension(3) :: delta
|
real(pReal), dimension(3) :: delta
|
||||||
real(pReal), dimension(3,8) :: BMat
|
real(pReal), dimension(3,8) :: BMat
|
||||||
real(pReal), dimension(8,8) :: HGMat
|
real(pReal), dimension(8,8) :: HGMat
|
||||||
PetscInt :: xstart,ystart,zstart,xend,yend,zend
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! stress, stiffness and compliance average etc.
|
! stress, stiffness and compliance average etc.
|
||||||
|
@ -146,7 +145,7 @@ subroutine grid_mechanical_FEM_init
|
||||||
! set default and user defined options for PETSc
|
! set default and user defined options for PETSc
|
||||||
call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
|
call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
|
||||||
'-mechanical_snes_type newtonls -mechanical_ksp_type fgmres &
|
'-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)
|
err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),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
|
! 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)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',err_PETSc)
|
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
localK = 0_pPetscInt
|
localK = 0_pPetscInt
|
||||||
localK(worldrank) = int(grid3,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
|
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
|
||||||
mechanical_grid,err_PETSc)
|
mechanical_grid,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetDM(mechanical_snes,mechanical_grid,err_PETSc)
|
|
||||||
CHKERRQ(err_PETSc)
|
|
||||||
call DMsetFromOptions(mechanical_grid,err_PETSc)
|
call DMsetFromOptions(mechanical_grid,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMsetUp(mechanical_grid,err_PETSc)
|
call DMsetUp(mechanical_grid,err_PETSc)
|
||||||
|
@ -195,28 +192,28 @@ subroutine grid_mechanical_FEM_init
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,err_PETSc)
|
call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,err_PETSc)
|
||||||
CHKERRQ(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)
|
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)
|
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)
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! init fields
|
! init fields
|
||||||
call VecSet(solution_current,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
|
call VecSet(solution_current,0.0_pReal,err_PETSc)
|
||||||
call VecSet(solution_lastInc,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call VecSet(solution_rate ,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)
|
call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc)
|
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,err_PETSc)
|
||||||
CHKERRQ(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
|
delta = geomSize/real(grid,pReal) ! grid spacing
|
||||||
detJ = product(delta) ! cell volume
|
detJ = product(delta) ! cell volume
|
||||||
|
|
||||||
|
@ -311,14 +308,9 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
|
||||||
! update stiffness (and gamma operator)
|
! update stiffness (and gamma operator)
|
||||||
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
|
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
call SNESsolve(SNES_mechanical,PETSC_NULL_VEC,solution_current,err_PETSc)
|
||||||
! solve BVP
|
|
||||||
call SNESsolve(mechanical_snes,PETSC_NULL_VEC,solution_current,err_PETSc)
|
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
! check convergence
|
|
||||||
call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc)
|
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
solution%converged = reason > 0
|
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)
|
call VecScale(solution_rate,1.0_pReal/Delta_t_old,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
else
|
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
|
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
|
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,')'
|
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
|
||||||
print'(/,1x,a)', '==========================================================================='
|
print'(/,1x,a)', '==========================================================================='
|
||||||
flush(IO_STDOUT)
|
flush(IO_STDOUT)
|
||||||
|
err_PETSc = 0
|
||||||
|
|
||||||
end subroutine converged
|
end subroutine converged
|
||||||
|
|
||||||
|
@ -527,7 +522,7 @@ subroutine formResidual(da_local,x_local, &
|
||||||
|
|
||||||
DM :: da_local
|
DM :: da_local
|
||||||
Vec :: x_local, f_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
|
PetscScalar, dimension(8,3) :: x_elem, f_elem
|
||||||
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
|
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
|
||||||
PetscInt :: &
|
PetscInt :: &
|
||||||
|
@ -538,9 +533,9 @@ subroutine formResidual(da_local,x_local, &
|
||||||
integer(MPI_INTEGER_KIND) :: err_MPI
|
integer(MPI_INTEGER_KIND) :: err_MPI
|
||||||
real(pReal), dimension(3,3,3,3) :: devNull
|
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)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESGetIterationNumber(mechanical_snes,PETScIter,err_PETSc)
|
call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
|
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
|
||||||
|
@ -559,17 +554,18 @@ subroutine formResidual(da_local,x_local, &
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! get deformation gradient
|
! get deformation gradient
|
||||||
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc);CHKERRQ(err_PETSc)
|
call DMDAVecGetArrayF90(da_local,x_local,x_scal,err_PETSc)
|
||||||
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
|
CHKERRQ(err_PETSc)
|
||||||
|
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
|
||||||
ctr = 0
|
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
|
ctr = ctr + 1
|
||||||
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
|
F(1:3,1:3,i,j,k-grid3offset) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
|
||||||
F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
|
|
||||||
enddo; enddo; enddo
|
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
|
! evaluate constitutive response
|
||||||
|
@ -586,47 +582,53 @@ subroutine formResidual(da_local,x_local, &
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! constructing residual
|
! constructing residual
|
||||||
call VecSet(f_local,0.0_pReal,err_PETSc);CHKERRQ(err_PETSc)
|
call VecSet(f_local,0.0_pReal,err_PETSc)
|
||||||
call DMDAVecGetArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMDAVecGetArrayF90(da_local,x_local,x_scal,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
|
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
|
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
|
ctr = ctr + 1
|
||||||
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
|
|
||||||
ele = ele + 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) + &
|
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,ele) + &
|
||||||
homogenization_dPdF(2,2,2,2,ele) + &
|
homogenization_dPdF(2,2,2,2,ele) + &
|
||||||
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
|
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
|
||||||
ctr = 0
|
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
|
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
|
||||||
enddo; enddo; enddo
|
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)
|
||||||
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! applying boundary conditions
|
! applying boundary conditions
|
||||||
call DMDAVecGetArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
|
call DMDAVecGetArrayF90(da_local,f_local,r,err_PETSc)
|
||||||
if (zstart == 0) then
|
CHKERRQ(err_PETSc)
|
||||||
f_scal(0:2,xstart,ystart,zstart) = 0.0
|
if (grid3offset == 0) then
|
||||||
f_scal(0:2,xend+1,ystart,zstart) = 0.0
|
r(0:2,0, 0, 0) = 0.0_pReal
|
||||||
f_scal(0:2,xstart,yend+1,zstart) = 0.0
|
r(0:2,grid(1),0, 0) = 0.0_pReal
|
||||||
f_scal(0:2,xend+1,yend+1,zstart) = 0.0
|
r(0:2,0, grid(2),0) = 0.0_pReal
|
||||||
endif
|
r(0:2,grid(1),grid(2),0) = 0.0_pReal
|
||||||
if (zend + 1 == grid(3)) then
|
end if
|
||||||
f_scal(0:2,xstart,ystart,zend+1) = 0.0
|
if (grid3+grid3offset == grid(3)) then
|
||||||
f_scal(0:2,xend+1,ystart,zend+1) = 0.0
|
r(0:2,0, 0, grid(3)) = 0.0_pReal
|
||||||
f_scal(0:2,xstart,yend+1,zend+1) = 0.0
|
r(0:2,grid(1),0, grid(3)) = 0.0_pReal
|
||||||
f_scal(0:2,xend+1,yend+1,zend+1) = 0.0
|
r(0:2,0, grid(2),grid(3)) = 0.0_pReal
|
||||||
endif
|
r(0:2,grid(1),grid(2),grid(3)) = 0.0_pReal
|
||||||
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,err_PETSc);CHKERRQ(err_PETSc)
|
end if
|
||||||
|
call DMDAVecRestoreArrayF90(da_local,f_local,r,err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
end subroutine formResidual
|
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,pointer,dimension(:,:,:,:) :: x_scal
|
||||||
PetscScalar,dimension(24,24) :: K_ele
|
PetscScalar,dimension(24,24) :: K_ele
|
||||||
PetscScalar,dimension(9,24) :: BMatFull
|
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]
|
PetscInt,dimension(3),parameter :: rows = [0, 1, 2]
|
||||||
PetscScalar :: diag
|
PetscScalar :: diag
|
||||||
PetscObject :: dummy
|
PetscObject :: dummy
|
||||||
|
@ -658,11 +660,12 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,err_PETSc)
|
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call MatZeroEntries(Jac,err_PETSc); CHKERRQ(err_PETSc)
|
call MatZeroEntries(Jac,err_PETSc)
|
||||||
ele = 0
|
CHKERRQ(err_PETSc)
|
||||||
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
|
ce = 0
|
||||||
|
do k = grid3offset+1, grid3offset+grid3; do j = 1, grid(2); do i = 1, grid(1)
|
||||||
ctr = 0
|
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
|
ctr = ctr + 1
|
||||||
col(MatStencil_i,ctr ) = i+ii
|
col(MatStencil_i,ctr ) = i+ii
|
||||||
col(MatStencil_j,ctr ) = j+jj
|
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
|
col(MatStencil_c,ctr+16) = 2
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
row = col
|
row = col
|
||||||
ele = ele + 1
|
ce = ce + 1
|
||||||
K_ele = 0.0
|
K_ele = 0.0
|
||||||
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
|
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
|
||||||
homogenization_dPdF(2,2,2,2,ele) + &
|
homogenization_dPdF(2,2,2,2,ce) + &
|
||||||
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
|
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
|
||||||
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
|
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
|
||||||
homogenization_dPdF(2,2,2,2,ele) + &
|
homogenization_dPdF(2,2,2,2,ce) + &
|
||||||
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
|
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
|
||||||
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ele) + &
|
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,ce) + &
|
||||||
homogenization_dPdF(2,2,2,2,ele) + &
|
homogenization_dPdF(2,2,2,2,ce) + &
|
||||||
homogenization_dPdF(3,3,3,3,ele))/3.0_pReal
|
homogenization_dPdF(3,3,3,3,ce))/3.0_pReal
|
||||||
K_ele = K_ele + &
|
K_ele = K_ele + &
|
||||||
matmul(transpose(BMatFull), &
|
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
|
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)
|
call MatSetValuesStencil(Jac,24_pPETScInt,row,24_pPetscInt,col,K_ele,ADD_VALUES,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
|
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
|
||||||
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
|
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
|
||||||
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(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
|
! applying boundary conditions
|
||||||
diag = (C_volAvg(1,1,1,1)/delta(1)**2 + &
|
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) &
|
||||||
C_volAvg(2,2,2,2)/delta(2)**2 + &
|
* detJ
|
||||||
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)
|
call MatZeroRowsColumns(Jac,size(rows,kind=pPetscInt),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMGetGlobalVector(da_local,coordinates,err_PETSc)
|
call DMGetGlobalVector(da_local,coordinates,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc)
|
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
ele = 0
|
ce = 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)
|
||||||
ele = ele + 1
|
ce = ce + 1
|
||||||
x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele)
|
x_scal(0:2,i-1,j-1,k-1) = discretization_IPcoords(1:3,ce)
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,err_PETSc)
|
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)
|
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)
|
call DMRestoreGlobalVector(da_local,coordinates,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call MatSetNullSpace(Jac,matnull,err_PETSc)
|
call MatSetNullSpace(Jac,matnull,err_PETSc)
|
||||||
|
|
|
@ -50,7 +50,7 @@ module grid_mechanical_spectral_basic
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc data
|
! PETSc data
|
||||||
DM :: da
|
DM :: da
|
||||||
SNES :: snes
|
SNES :: SNES_mechanical
|
||||||
Vec :: solution_vec
|
Vec :: solution_vec
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -158,8 +158,10 @@ subroutine grid_mechanical_spectral_basic_init
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! initialize solver specific parts of PETSc
|
! initialize solver specific parts of PETSc
|
||||||
call SNESCreate(PETSC_COMM_WORLD,snes,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
|
||||||
call SNESSetOptionsPrefix(snes,'mechanical_',err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
localK = 0_pPetscInt
|
localK = 0_pPetscInt
|
||||||
localK(worldrank) = int(grid3,pPetscInt)
|
localK(worldrank) = int(grid3,pPetscInt)
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
|
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
|
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
|
||||||
da,err_PETSc) ! handle, error
|
da,err_PETSc) ! handle, error
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetDM(snes,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
|
call DMsetFromOptions(da,err_PETSc)
|
||||||
call DMsetFromOptions(da,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMsetUp(da,err_PETSc); CHKERRQ(err_PETSc)
|
call DMsetUp(da,err_PETSc)
|
||||||
call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor)
|
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
|
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
|
||||||
CHKERRQ(err_PETSc)
|
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)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! init fields
|
! 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
|
restartRead: if (interface_restartInc > 0) then
|
||||||
print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
|
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
|
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
|
reshape(F,shape(F_lastInc)), & ! target F
|
||||||
0.0_pReal) ! time increment
|
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
|
restartRead2: if (interface_restartInc > 0) then
|
||||||
print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
|
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)
|
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
|
||||||
if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
|
if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
call SNESsolve(SNES_mechanical,PETSC_NULL_VEC,solution_vec,err_PETSc)
|
||||||
! solve BVP
|
CHKERRQ(err_PETSc)
|
||||||
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
! check convergence
|
|
||||||
call SNESGetConvergedReason(snes,reason,err_PETSc); CHKERRQ(err_PETSc)
|
|
||||||
|
|
||||||
solution%converged = reason > 0
|
solution%converged = reason > 0
|
||||||
solution%iterationsNeeded = totalIter
|
solution%iterationsNeeded = totalIter
|
||||||
|
@ -310,7 +316,8 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
|
||||||
PetscScalar, pointer, dimension(:,:,:,:) :: F
|
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
|
if (cutBack) then
|
||||||
C_volAvg = C_volAvgLastInc
|
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
|
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])
|
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
|
! set module wide available data
|
||||||
|
@ -372,9 +380,11 @@ subroutine grid_mechanical_spectral_basic_updateCoords
|
||||||
PetscErrorCode :: err_PETSc
|
PetscErrorCode :: err_PETSc
|
||||||
PetscScalar, dimension(:,:,:,:), pointer :: F
|
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 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
|
end subroutine grid_mechanical_spectral_basic_updateCoords
|
||||||
|
|
||||||
|
@ -388,7 +398,8 @@ subroutine grid_mechanical_spectral_basic_restartWrite
|
||||||
integer(HID_T) :: fileHandle, groupHandle
|
integer(HID_T) :: fileHandle, groupHandle
|
||||||
PetscScalar, dimension(:,:,:,:), pointer :: F
|
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)
|
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
|
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
|
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,')'
|
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
|
||||||
print'(/,1x,a)', '==========================================================================='
|
print'(/,1x,a)', '==========================================================================='
|
||||||
flush(IO_STDOUT)
|
flush(IO_STDOUT)
|
||||||
|
err_PETSc = 0
|
||||||
|
|
||||||
end subroutine converged
|
end subroutine converged
|
||||||
|
|
||||||
|
@ -465,13 +478,13 @@ end subroutine converged
|
||||||
!> @brief forms the residual vector
|
!> @brief forms the residual vector
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine formResidual(in, F, &
|
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)
|
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), &
|
PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
|
||||||
intent(in) :: F !< deformation gradient field
|
intent(in) :: F !< deformation gradient field
|
||||||
PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), &
|
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) :: &
|
real(pReal), dimension(3,3) :: &
|
||||||
deltaF_aim
|
deltaF_aim
|
||||||
PetscInt :: &
|
PetscInt :: &
|
||||||
|
@ -481,8 +494,10 @@ subroutine formResidual(in, F, &
|
||||||
PetscErrorCode :: err_PETSc
|
PetscErrorCode :: err_PETSc
|
||||||
integer(MPI_INTEGER_KIND) :: err_MPI
|
integer(MPI_INTEGER_KIND) :: err_MPI
|
||||||
|
|
||||||
call SNESGetNumberFunctionEvals(snes,nfuncs,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
|
||||||
call SNESGetIterationNumber(snes,PETScIter,err_PETSc); CHKERRQ(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
|
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
|
||||||
|
|
||||||
|
@ -500,7 +515,7 @@ subroutine formResidual(in, F, &
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! evaluate constitutive response
|
! 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, &
|
P_av,C_volAvg,C_minMaxAvg, &
|
||||||
F,params%Delta_t,params%rotation_BC)
|
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)
|
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
|
! 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) = r ! 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(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier
|
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
|
! 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
|
end subroutine formResidual
|
||||||
|
|
||||||
|
|
|
@ -55,7 +55,7 @@ module grid_mechanical_spectral_polarisation
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc data
|
! PETSc data
|
||||||
DM :: da
|
DM :: da
|
||||||
SNES :: snes
|
SNES :: SNES_mechanical
|
||||||
Vec :: solution_vec
|
Vec :: solution_vec
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -178,8 +178,10 @@ subroutine grid_mechanical_spectral_polarisation_init
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! initialize solver specific parts of PETSc
|
! initialize solver specific parts of PETSc
|
||||||
call SNESCreate(PETSC_COMM_WORLD,snes,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESCreate(PETSC_COMM_WORLD,SNES_mechanical,err_PETSc)
|
||||||
call SNESSetOptionsPrefix(snes,'mechanical_',err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call SNESSetOptionsPrefix(SNES_mechanical,'mechanical_',err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
localK = 0_pPetscInt
|
localK = 0_pPetscInt
|
||||||
localK(worldrank) = int(grid3,pPetscInt)
|
localK(worldrank) = int(grid3,pPetscInt)
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
|
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
|
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
|
||||||
da,err_PETSc) ! handle, error
|
da,err_PETSc) ! handle, error
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetDM(snes,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
|
call DMsetFromOptions(da,err_PETSc)
|
||||||
call DMsetFromOptions(da,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMsetUp(da,err_PETSc); CHKERRQ(err_PETSc)
|
call DMsetUp(da,err_PETSc)
|
||||||
call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor)
|
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
|
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
|
||||||
CHKERRQ(err_PETSc)
|
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)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! init fields
|
! 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 => FandF_tau(0: 8,:,:,:)
|
||||||
F_tau => FandF_tau(9:17,:,:,:)
|
F_tau => FandF_tau(9:17,:,:,:)
|
||||||
|
|
||||||
|
@ -285,7 +293,7 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
|
||||||
! input data for solution
|
! input data for solution
|
||||||
character(len=*), intent(in) :: &
|
character(len=*), intent(in) :: &
|
||||||
incInfoIn
|
incInfoIn
|
||||||
type(tSolutionState) :: &
|
type(tSolutionState) :: &
|
||||||
solution
|
solution
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc Data
|
! PETSc Data
|
||||||
|
@ -303,13 +311,10 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
|
||||||
S_scale = math_invSym3333(C_minMaxAvg)
|
S_scale = math_invSym3333(C_minMaxAvg)
|
||||||
end if
|
end if
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
call SNESSolve(SNES_mechanical,PETSC_NULL_VEC,solution_vec,err_PETSc)
|
||||||
! solve BVP
|
CHKERRQ(err_PETSc)
|
||||||
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESGetConvergedReason(SNES_mechanical,reason,err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
!--------------------------------------------------------------------------------------------------
|
|
||||||
! check convergence
|
|
||||||
call SNESGetConvergedReason(snes,reason,err_PETSc); CHKERRQ(err_PETSc)
|
|
||||||
|
|
||||||
solution%converged = reason > 0
|
solution%converged = reason > 0
|
||||||
solution%iterationsNeeded = totalIter
|
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
|
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 => FandF_tau(0: 8,:,:,:)
|
||||||
F_tau => FandF_tau(9:17,:,:,:)
|
F_tau => FandF_tau(9:17,:,:,:)
|
||||||
|
|
||||||
|
@ -446,7 +452,8 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
|
||||||
integer(HID_T) :: fileHandle, groupHandle
|
integer(HID_T) :: fileHandle, groupHandle
|
||||||
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
|
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 => FandF_tau(0: 8,:,:,:)
|
||||||
F_tau => FandF_tau(9:17,:,:,:)
|
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,')'
|
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
|
||||||
print'(/,1x,a)', '==========================================================================='
|
print'(/,1x,a)', '==========================================================================='
|
||||||
flush(IO_STDOUT)
|
flush(IO_STDOUT)
|
||||||
|
err_PETSc = 0
|
||||||
|
|
||||||
end subroutine converged
|
end subroutine converged
|
||||||
|
|
||||||
|
@ -531,18 +539,18 @@ end subroutine converged
|
||||||
!> @brief forms the residual vector
|
!> @brief forms the residual vector
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine formResidual(in, FandF_tau, &
|
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)
|
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), &
|
PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), &
|
||||||
target, intent(in) :: FandF_tau
|
target, intent(in) :: FandF_tau
|
||||||
PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),&
|
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(:,:,:,:,:) :: &
|
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
|
||||||
F, &
|
F, &
|
||||||
F_tau, &
|
F_tau, &
|
||||||
residual_F, &
|
r_F, &
|
||||||
residual_F_tau
|
r_F_tau
|
||||||
PetscInt :: &
|
PetscInt :: &
|
||||||
PETScIter, &
|
PETScIter, &
|
||||||
nfuncs
|
nfuncs
|
||||||
|
@ -554,21 +562,23 @@ subroutine formResidual(in, FandF_tau, &
|
||||||
|
|
||||||
!---------------------------------------------------------------------------------------------------
|
!---------------------------------------------------------------------------------------------------
|
||||||
|
|
||||||
F => FandF_tau(1:3,1:3,1,&
|
F => FandF_tau(1:3,1:3,1,&
|
||||||
XG_RANGE,YG_RANGE,ZG_RANGE)
|
XG_RANGE,YG_RANGE,ZG_RANGE)
|
||||||
F_tau => FandF_tau(1:3,1:3,2,&
|
F_tau => FandF_tau(1:3,1:3,2,&
|
||||||
XG_RANGE,YG_RANGE,ZG_RANGE)
|
XG_RANGE,YG_RANGE,ZG_RANGE)
|
||||||
residual_F => residuum(1:3,1:3,1,&
|
r_F => r(1:3,1:3,1,&
|
||||||
X_RANGE, Y_RANGE, Z_RANGE)
|
X_RANGE, Y_RANGE, Z_RANGE)
|
||||||
residual_F_tau => residuum(1:3,1:3,2,&
|
r_F_tau => r(1:3,1:3,2,&
|
||||||
X_RANGE, Y_RANGE, Z_RANGE)
|
X_RANGE, Y_RANGE, Z_RANGE)
|
||||||
|
|
||||||
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
|
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)
|
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'
|
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
|
||||||
|
|
||||||
call SNESGetNumberFunctionEvals(snes,nfuncs,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc)
|
||||||
call SNESGetIterationNumber(snes,PETScIter,err_PETSc); CHKERRQ(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
|
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
|
||||||
|
|
||||||
|
@ -602,13 +612,13 @@ subroutine formResidual(in, FandF_tau, &
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! constructing residual
|
! 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
|
! 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, &
|
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)
|
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)))
|
params%stress_mask)))
|
||||||
! calculate divergence
|
! calculate divergence
|
||||||
tensorField_real = 0.0_pReal
|
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
|
call utilities_FFTtensorForward
|
||||||
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
|
err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress
|
||||||
|
|
||||||
|
@ -628,11 +638,11 @@ subroutine formResidual(in, FandF_tau, &
|
||||||
e = 0
|
e = 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)
|
||||||
e = e + 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), &
|
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))) &
|
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
|
end do; end do; end do
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
|
|
@ -38,9 +38,8 @@ module grid_thermal_spectral
|
||||||
type(tSolutionParams) :: params
|
type(tSolutionParams) :: params
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! PETSc data
|
! PETSc data
|
||||||
SNES :: thermal_snes
|
SNES :: SNES_thermal
|
||||||
Vec :: solution_vec
|
Vec :: solution_vec
|
||||||
PetscInt :: xstart, xend, ystart, yend, zstart, zend
|
|
||||||
real(pReal), dimension(:,:,:), allocatable :: &
|
real(pReal), dimension(:,:,:), allocatable :: &
|
||||||
T_current, & !< field of current temperature
|
T_current, & !< field of current temperature
|
||||||
T_lastInc, & !< field of previous 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)
|
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),err_PETSc)
|
||||||
CHKERRQ(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
|
! initialize solver specific parts of PETSc
|
||||||
call SNESCreate(PETSC_COMM_WORLD,thermal_snes,err_PETSc); CHKERRQ(err_PETSc)
|
call SNESCreate(PETSC_COMM_WORLD,SNES_thermal,err_PETSc)
|
||||||
call SNESSetOptionsPrefix(thermal_snes,'thermal_',err_PETSc);CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call SNESSetOptionsPrefix(SNES_thermal,'thermal_',err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
localK = 0_pPetscInt
|
localK = 0_pPetscInt
|
||||||
localK(worldrank) = int(grid3,pPetscInt)
|
localK(worldrank) = int(grid3,pPetscInt)
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
|
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
|
[int(grid(1),pPetscInt)],[int(grid(2),pPetscInt)],localK, & ! local grid
|
||||||
thermal_grid,err_PETSc) ! handle, error
|
thermal_grid,err_PETSc) ! handle, error
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetDM(thermal_snes,thermal_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da
|
call DMsetFromOptions(thermal_grid,err_PETSc)
|
||||||
call DMsetFromOptions(thermal_grid,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMsetUp(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)
|
call DMCreateGlobalVector(thermal_grid,solution_vec,err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
|
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESSetFromOptions(thermal_snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments
|
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
|
||||||
! init fields
|
|
||||||
call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,err_PETSc)
|
|
||||||
CHKERRQ(err_PETSc)
|
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)
|
call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
|
||||||
CHKERRQ(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)
|
call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
call updateReference
|
call updateReference()
|
||||||
|
|
||||||
end subroutine grid_thermal_spectral_init
|
end subroutine grid_thermal_spectral_init
|
||||||
|
|
||||||
|
@ -179,9 +175,9 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
|
||||||
! set module wide availabe data
|
! set module wide availabe data
|
||||||
params%Delta_t = Delta_t
|
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)
|
CHKERRQ(err_PETSc)
|
||||||
call SNESGetConvergedReason(thermal_snes,reason,err_PETSc)
|
call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
|
||||||
if (reason < 1) then
|
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)
|
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
|
end do; end do; end do
|
||||||
|
|
||||||
call VecMin(solution_vec,devNull,T_min,err_PETSc); CHKERRQ(err_PETSc)
|
call VecMin(solution_vec,devNull,T_min,err_PETSc)
|
||||||
call VecMax(solution_vec,devNull,T_max,err_PETSc); CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
|
call VecMax(solution_vec,devNull,T_max,err_PETSc)
|
||||||
|
CHKERRQ(err_PETSc)
|
||||||
if (solution%converged) &
|
if (solution%converged) &
|
||||||
print'(/,1x,a)', '... thermal conduction 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
|
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
|
logical, intent(in) :: cutBack
|
||||||
integer :: i, j, k, ce
|
integer :: i, j, k, ce
|
||||||
DM :: dm_local
|
DM :: dm_local
|
||||||
PetscScalar, dimension(:,:,:), pointer :: x_scal
|
PetscScalar, dimension(:,:,:), pointer :: T_PETSc
|
||||||
PetscErrorCode :: err_PETSc
|
PetscErrorCode :: err_PETSc
|
||||||
|
|
||||||
if (cutBack) then
|
if (cutBack) then
|
||||||
|
@ -235,12 +233,12 @@ subroutine grid_thermal_spectral_forward(cutBack)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! reverting thermal field state
|
! reverting thermal field state
|
||||||
call SNESGetDM(thermal_snes,dm_local,err_PETSc)
|
call SNESGetDM(SNES_thermal,dm_local,err_PETSc)
|
||||||
CHKERRQ(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)
|
CHKERRQ(err_PETSc)
|
||||||
x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
|
T_PETSc = T_current
|
||||||
call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,err_PETSc)
|
call DMDAVecRestoreArrayF90(dm_local,solution_vec,T_PETSc,err_PETSc)
|
||||||
CHKERRQ(err_PETSc)
|
CHKERRQ(err_PETSc)
|
||||||
ce = 0
|
ce = 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)
|
||||||
|
@ -258,7 +256,7 @@ end subroutine grid_thermal_spectral_forward
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief forms the spectral thermal residual vector
|
!> @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) :: &
|
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
|
||||||
in
|
in
|
||||||
|
@ -267,9 +265,9 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
|
||||||
x_scal
|
x_scal
|
||||||
PetscScalar, dimension( &
|
PetscScalar, dimension( &
|
||||||
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: &
|
||||||
f_scal
|
r
|
||||||
PetscObject :: dummy
|
PetscObject :: dummy
|
||||||
PetscErrorCode :: dummy_err
|
PetscErrorCode :: err_PETSc
|
||||||
integer :: i, j, k, ce
|
integer :: i, j, k, ce
|
||||||
|
|
||||||
T_current = x_scal
|
T_current = x_scal
|
||||||
|
@ -304,7 +302,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,dummy_err)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! constructing residual
|
! 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
|
end subroutine formResidual
|
||||||
|
|
||||||
|
|
|
@ -53,6 +53,7 @@ contains
|
||||||
subroutine parallelization_init
|
subroutine parallelization_init
|
||||||
|
|
||||||
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
|
integer(MPI_INTEGER_KIND) :: err_MPI, typeSize
|
||||||
|
character(len=4) :: rank_str
|
||||||
!$ integer :: got_env, threadLevel
|
!$ integer :: got_env, threadLevel
|
||||||
!$ integer(pI32) :: OMP_NUM_THREADS
|
!$ integer(pI32) :: OMP_NUM_THREADS
|
||||||
!$ character(len=6) NumThreadsString
|
!$ character(len=6) NumThreadsString
|
||||||
|
@ -112,6 +113,7 @@ subroutine parallelization_init
|
||||||
|
|
||||||
if (worldrank /= 0) then
|
if (worldrank /= 0) then
|
||||||
close(OUTPUT_UNIT) ! disable output
|
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
|
open(OUTPUT_UNIT,file='/dev/null',status='replace') ! close() alone will leave some temp files in cwd
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue