diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index f699ee7de..48994ec15 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -132,7 +132,7 @@ subroutine grid_damage_spectral_init() call DMCreateGlobalVector(damage_grid,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,err_PETSc) ! residual vector of same shape as solution vector CHKERRQ(err_PETSc) - call SNESSetDM(SNES_damage,damage_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da + call SNESSetDM(SNES_damage,damage_grid,err_PETSc); CHKERRQ(err_PETSc) call SNESSetFromOptions(SNES_damage,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments call SNESGetType(SNES_damage,snes_type,err_PETSc); CHKERRQ(err_PETSc) if (trim(snes_type) == 'vinewtonrsls' .or. & @@ -141,12 +141,10 @@ subroutine grid_damage_spectral_init() 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 like contact mechanics, damage etc. + call SNESVISetVariableBounds(SNES_damage,lBound,uBound,err_PETSc) ! variable bounds for variational inequalities call DMRestoreGlobalVector(damage_grid,lBound,err_PETSc); CHKERRQ(err_PETSc) call DMRestoreGlobalVector(damage_grid,uBound,err_PETSc); CHKERRQ(err_PETSc) end if - -!-------------------------------------------------------------------------------------------------- call VecSet(solution_vec,1.0_pReal,err_PETSc); CHKERRQ(err_PETSc) call updateReference() diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index a697ce14a..fb3ec85f6 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 !-------------------------------------------------------------------------------------------------- @@ -159,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) @@ -192,13 +192,13 @@ 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 SNESSetDM(mechanical_snes,mechanical_grid,err_PETSc) + call SNESSetDM(SNES_mechanical,mechanical_grid,err_PETSc) CHKERRQ(err_PETSc) - call SNESSetFromOptions(mechanical_snes,err_PETSc) ! pull it all together with additional cli arguments + call SNESSetFromOptions(SNES_mechanical,err_PETSc) ! pull it all together with additional cli arguments CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- @@ -305,14 +305,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 @@ -533,9 +528,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 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index bf72d1d08..00ff756ec 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) @@ -178,10 +180,10 @@ subroutine grid_mechanical_spectral_basic_init call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 9, i.e. every def grad tensor) 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,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da - call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments + call SNESSetDM(SNES_mechanical,da,err_PETSc); CHKERRQ(err_PETSc) + call SNESsetFromOptions(SNES_mechanical,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- ! init fields @@ -270,13 +272,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 @@ -482,8 +481,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 diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9a258e13e..52c308ae9 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,9 @@ 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) @@ -198,10 +199,10 @@ subroutine grid_mechanical_spectral_polarisation_init call DMcreateGlobalVector(da,solution_vec,err_PETSc); CHKERRQ(err_PETSc) ! global solution vector (grid x 18, i.e. every def grad tensor) 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,da,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da - call SNESsetFromOptions(snes,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments + call SNESSetDM(SNES_mechanical,da,err_PETSc); CHKERRQ(err_PETSc) + call SNESsetFromOptions(SNES_mechanical,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments !-------------------------------------------------------------------------------------------------- ! init fields @@ -285,7 +286,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 +304,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 @@ -542,8 +540,8 @@ subroutine formResidual(in, FandF_tau, & PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_tau, & - residual_F, & - residual_F_tau + r_F, & + r_F_tau PetscInt :: & PETScIter, & nfuncs @@ -555,21 +553,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 => r(1:3,1:3,1,& - X_RANGE, Y_RANGE, Z_RANGE) - residual_F_tau => r(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 @@ -603,13 +603,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) !-------------------------------------------------------------------------------------------------- @@ -620,7 +620,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 @@ -629,11 +629,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 d1d165b83..9c6d10cca 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -134,10 +134,8 @@ subroutine grid_thermal_spectral_init(T_0) 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 SNESSetDM(SNES_thermal,thermal_grid,err_PETSc); CHKERRQ(err_PETSc) ! connect snes to da + call SNESSetDM(SNES_thermal,thermal_grid,err_PETSc); CHKERRQ(err_PETSc) call SNESSetFromOptions(SNES_thermal,err_PETSc); CHKERRQ(err_PETSc) ! pull it all together with additional CLI arguments - -!-------------------------------------------------------------------------------------------------- call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,err_PETSc) CHKERRQ(err_PETSc) T_PETSc = T_current