diff --git a/src/grid_mech_spectral_basic.f90 b/src/grid_mech_spectral_basic.f90 index 295c16b98..f5fdcc21f 100644 --- a/src/grid_mech_spectral_basic.f90 +++ b/src/grid_mech_spectral_basic.f90 @@ -64,6 +64,7 @@ module grid_mech_spectral_basic grid_mech_spectral_basic_solution, & grid_mech_spectral_basic_forward private :: & + converged, & formResidual contains @@ -151,7 +152,7 @@ subroutine grid_mech_spectral_basic_init call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 9, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "_converged" + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -267,7 +268,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ solution%termIll = terminallyIll terminallyIll = .false. - end function grid_mech_spectral_basic_solution diff --git a/src/grid_mech_spectral_polarisation.f90 b/src/grid_mech_spectral_polarisation.f90 index 3fca8fc9f..b7ab8ba12 100644 --- a/src/grid_mech_spectral_polarisation.f90 +++ b/src/grid_mech_spectral_polarisation.f90 @@ -19,10 +19,7 @@ module grid_mech_spectral_polarisation implicit none private - - character (len=*), parameter, public :: & - GRID_MECH_SPECTRAL_POLARISATION_LABEL = 'polarisation' - + !-------------------------------------------------------------------------------------------------- ! derived types type(tSolutionParams), private :: params @@ -74,6 +71,7 @@ module grid_mech_spectral_polarisation grid_mech_spectral_polarisation_solution, & grid_mech_spectral_polarisation_forward private :: & + converged, & formResidual contains @@ -86,10 +84,6 @@ subroutine grid_mech_spectral_polarisation_init IO_intOut, & IO_error, & IO_open_jobFile_binary - use debug, only: & - debug_level, & - debug_spectral, & - debug_spectralRestart use FEsolving, only: & restartInc use numerics, only: & @@ -166,7 +160,7 @@ subroutine grid_mech_spectral_polarisation_init call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr) - call SNESsetConvergenceTest(snes,grid_mech_spectral_polarisation_converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" + call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" CHKERRQ(ierr) call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments @@ -175,12 +169,9 @@ subroutine grid_mech_spectral_polarisation_init call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data F => FandF_tau( 0: 8,:,:,:) F_tau => FandF_tau( 9:17,:,:,:) + restart: if (restartInc > 0) then - if (iand(debug_level(debug_spectral),debug_spectralRestart) /= 0) then - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & - 'reading values of increment ', restartInc, ' from file' - flush(6) - endif + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file' fileUnit = IO_open_jobFile_binary('F_aimDot') read(fileUnit) F_aimDot; close(fileUnit) @@ -218,10 +209,8 @@ subroutine grid_mech_spectral_polarisation_init call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer restartRead: if (restartInc > 0) then - if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0) & - write(6,'(/,a,'//IO_intOut(restartInc)//',a)') & - 'reading more values of increment ', restartInc, ' from file' - flush(6) + write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file' + fileUnit = IO_open_jobFile_binary('C_volAvg') read(fileUnit) C_volAvg; close(fileUnit) fileUnit = IO_open_jobFile_binary('C_volAvgLastInv') @@ -309,77 +298,6 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, end function grid_mech_spectral_polarisation_solution -!-------------------------------------------------------------------------------------------------- -!> @brief convergence check -!-------------------------------------------------------------------------------------------------- -subroutine grid_mech_spectral_polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) - use numerics, only: & - itmax, & - itmin, & - err_div_tolRel, & - err_div_tolAbs, & - err_curl_tolRel, & - err_curl_tolAbs, & - err_stress_tolRel, & - err_stress_tolAbs - use math, only: & - math_mul3333xx33 - use FEsolving, only: & - terminallyIll - - implicit none - SNES :: snes_local - PetscInt :: PETScIter - PetscReal :: & - xnorm, & ! not used - snorm, & ! not used - fnorm ! not used - SNESConvergedReason :: reason - PetscObject :: dummy - PetscErrorCode :: ierr - real(pReal) :: & - curlTol, & - divTol, & - BCTol - -!-------------------------------------------------------------------------------------------------- -! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim-F_av) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc - -!-------------------------------------------------------------------------------------------------- -! error calculation - curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) - divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) - BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) - - converged: if ((totalIter >= itmin .and. & - all([ err_div /divTol, & - err_curl/curlTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then - reason = 1 - elseif (totalIter >= itmax) then converged - reason = -1 - else converged - reason = 0 - endif converged - -!-------------------------------------------------------------------------------------------------- -! report - write(6,'(1/,a)') ' ... reporting .............................................................' - write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & - err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & - err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' - write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & - err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' - write(6,'(/,a)') ' ===========================================================================' - flush(6) - -end subroutine grid_mech_spectral_polarisation_converged - !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !> @details find new boundary conditions and best F estimate for end of current timestep @@ -467,7 +385,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa write(fileUnit) F_tau_lastInc; close (fileUnit) endif - call CPFEM_age() ! age state and kinematics + call CPFEM_age ! age state and kinematics call utilities_updateIPcoords(F) C_volAvgLastInc = C_volAvg @@ -504,7 +422,6 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average math_rotate_backward33(F_aim,rotation_BC)),& [9,grid(1),grid(2),grid3]) @@ -530,6 +447,76 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa end subroutine grid_mech_spectral_polarisation_forward +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check +!-------------------------------------------------------------------------------------------------- +subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) + use numerics, only: & + itmax, & + itmin, & + err_div_tolRel, & + err_div_tolAbs, & + err_curl_tolRel, & + err_curl_tolAbs, & + err_stress_tolRel, & + err_stress_tolAbs + use math, only: & + math_mul3333xx33 + use FEsolving, only: & + terminallyIll + + implicit none + SNES :: snes_local + PetscInt :: PETScIter + PetscReal :: & + xnorm, & ! not used + snorm, & ! not used + fnorm ! not used + SNESConvergedReason :: reason + PetscObject :: dummy + PetscErrorCode :: ierr + real(pReal) :: & + curlTol, & + divTol, & + BCTol + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc + err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim-F_av) + & + params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + +!-------------------------------------------------------------------------------------------------- +! error calculation + curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) + divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) + BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) + + if ((totalIter >= itmin .and. & + all([ err_div /divTol, & + err_curl/curlTol, & + err_BC /BCTol ] < 1.0_pReal)) & + .or. terminallyIll) then + reason = 1 + elseif (totalIter >= itmax) then + reason = -1 + else + reason = 0 + endif + +!-------------------------------------------------------------------------------------------------- +! report + write(6,'(1/,a)') ' ... reporting .............................................................' + write(6,'(/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & + err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & + err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' + write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & + err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' + write(6,'(/,a)') ' ===========================================================================' + flush(6) + +end subroutine converged !-------------------------------------------------------------------------------------------------- !> @brief forms the polarisation residual vector !-------------------------------------------------------------------------------------------------- @@ -570,9 +557,9 @@ subroutine formResidual(in, FandF_tau, & implicit none 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 - 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 PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & @@ -668,7 +655,7 @@ subroutine formResidual(in, FandF_tau, & ! calculating curl tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F - call utilities_FFTtensorForward() + call utilities_FFTtensorForward err_curl = Utilities_curlRMS() end subroutine formResidual