diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 738206f9a..8de2de6ef 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -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 :: & diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 8955f9d66..96b72dbcc 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -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 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 010f11cea..6a3f1323d 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -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) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 74400f160..c60c542d0 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -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 diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 765a1a5a7..810047403 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 8862df725..f102cf2c1 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -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 diff --git a/src/parallelization.f90 b/src/parallelization.f90 index d9ca15981..29deaf724 100644 --- a/src/parallelization.f90 +++ b/src/parallelization.f90 @@ -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