diff --git a/code/spectral_damage.f90 b/code/spectral_damage.f90 index 27d513175..6260fe43a 100644 --- a/code/spectral_damage.f90 +++ b/code/spectral_damage.f90 @@ -42,7 +42,6 @@ module spectral_damage integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment real(pReal), dimension(3,3), private :: D_ref real(pReal), private :: mobility_ref - character(len=1024), private :: incInfo public :: & spectral_damage_init, & @@ -50,21 +49,7 @@ module spectral_damage spectral_damage_forward, & spectral_damage_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASNESSetFunctionLocal, & PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & MPI_Abort, & MPI_Bcast, & MPI_Allreduce @@ -90,15 +75,30 @@ subroutine spectral_damage_init() damage_nonlocal_getMobility implicit none + integer(pInt), dimension(:), allocatable :: localK + integer(pInt) :: proc + integer(pInt) :: i, j, k, cell DM :: damage_grid Vec :: uBound, lBound PetscErrorCode :: ierr PetscObject :: dummy - integer(pInt), dimension(:), allocatable :: localK - integer(pInt) :: proc - integer(pInt) :: i, j, k, cell character(len=100) :: snes_type + external :: & + SNESCreate, & + SNESSetOptionsPrefix, & + DMDACreate3D, & + SNESSetDM, & + DMDAGetCorners, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + SNESSetFromOptions, & + SNESGetType, & + VecSet, & + DMGetGlobalVector, & + DMRestoreGlobalVector, & + SNESVISetVariableBounds + mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -194,12 +194,18 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old integer(pInt) :: i, j, k, cell PetscInt ::position PetscReal :: minDamage, maxDamage, stagNorm, solnNorm - + !-------------------------------------------------------------------------------------------------- ! PETSc Data PetscErrorCode :: ierr SNESConvergedReason :: reason - + + external :: & + VecMin, & + VecMax, & + SNESSolve, & + SNESGetConvergedReason + spectral_damage_solution%converged =.false. !-------------------------------------------------------------------------------------------------- @@ -353,10 +359,13 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime) timeinc, & loadCaseTime !< remaining time of current load case logical, intent(in) :: guess - PetscErrorCode :: ierr integer(pInt) :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal + PetscErrorCode :: ierr + + external :: & + SNESGetDM if (cutBack) then damage_current = damage_lastInc @@ -400,6 +409,10 @@ subroutine spectral_damage_destroy() implicit none PetscErrorCode :: ierr + external :: & + VecDestroy, & + SNESDestroy + call VecDestroy(solution,ierr); CHKERRQ(ierr) call SNESDestroy(damage_snes,ierr); CHKERRQ(ierr) diff --git a/code/spectral_mech_AL.f90 b/code/spectral_mech_AL.f90 index b6a8c9353..0077ba2c3 100644 --- a/code/spectral_mech_AL.f90 +++ b/code/spectral_mech_AL.f90 @@ -22,7 +22,7 @@ module spectral_mech_AL DAMASK_spectral_solverAL_label = 'al' !-------------------------------------------------------------------------------------------------- -! derived types +! derived types type(tSolutionParams), private :: params real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal @@ -31,7 +31,7 @@ module spectral_mech_AL DM, private :: da SNES, private :: snes Vec, private :: solution_vec - + !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: & @@ -72,21 +72,7 @@ module spectral_mech_AL AL_forward, & AL_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASNESSetFunctionLocal, & PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & MPI_Abort, & MPI_Bcast, & MPI_Allreduce @@ -136,11 +122,22 @@ subroutine AL_init integer(pInt) :: proc character(len=1024) :: rankStr - if (worldrank == 0_pInt) then + external :: & + SNESCreate, & + SNESSetOptionsPrefix, & + DMDACreate3D, & + SNESSetDM, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions + + mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif + endif mainProcess !-------------------------------------------------------------------------------------------------- ! allocate global fields @@ -150,7 +147,7 @@ subroutine AL_init allocate (F_lambdaDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- -! PETSc Init +! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 @@ -185,10 +182,10 @@ subroutine AL_init 'reading values of increment ', restartInc - 1_pInt, ' from file' flush(6) write(rankStr,'(a1,i0)')'_',worldrank - call IO_read_realFile(777,'F'//trim(rankStr), trim(getSolverJobName()),size(F)) + call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) read (777,rec=1) F close (777) - call IO_read_realFile(777,'F_lastInc'//trim(rankStr), trim(getSolverJobName()),size(F_lastInc)) + call IO_read_realFile(777,'F_lastInc'//trim(rankStr),trim(getSolverJobName()),size(F_lastInc)) read (777,rec=1) F_lastInc close (777) call IO_read_realFile(777,'F_lambda'//trim(rankStr),trim(getSolverJobName()),size(F_lambda)) @@ -214,15 +211,14 @@ subroutine AL_init F_lambda_lastInc = F_lastInc endif restart - call Utilities_updateIPcoords(reshape(F,shape(F_lastInc))) call Utilities_constitutiveResponse(F_lastInc, reshape(F,shape(F_lastInc)), & 0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3) nullify(F) nullify(F_lambda) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc - - readRestart: if (restartInc > 1_pInt) then + + restartRead: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' @@ -236,7 +232,7 @@ subroutine AL_init call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) - endif readRestart + endif restartRead call Utilities_updateGamma(C_minMaxAvg,.True.) C_scale = C_minMaxAvg @@ -263,7 +259,7 @@ type(tSolutionState) function & use FEsolving, only: & restartWrite, & terminallyIll - + implicit none !-------------------------------------------------------------------------------------------------- @@ -286,6 +282,10 @@ type(tSolutionState) function & PetscErrorCode :: ierr SNESConvergedReason :: reason + external :: & + SNESSolve, & + SNESGetConvergedReason + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- @@ -298,7 +298,7 @@ type(tSolutionState) function & endif !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data mask_stress = P_BC%maskFloat params%P_BC = P_BC%values params%rotation_BC = rotation_BC @@ -387,6 +387,10 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) integer(pInt) :: & i, j, k, e + external :: & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber + F => x_scal(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) F_lambda => x_scal(1:3,1:3,2,& @@ -414,7 +418,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & - math_transpose33(F_aim) + math_transpose33(F_aim) flush(6) endif endif newIteration @@ -507,7 +511,7 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr fnorm SNESConvergedReason :: reason PetscObject :: dummy - PetscErrorCode ::ierr + PetscErrorCode :: ierr real(pReal) :: & curlTol, & divTol, & @@ -704,6 +708,11 @@ subroutine AL_destroy() implicit none PetscErrorCode :: ierr + external :: & + VecDestroy, & + SNESDestroy, & + DMDestroy + call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) diff --git a/code/spectral_mech_Basic.f90 b/code/spectral_mech_Basic.f90 index 358a095d1..445fbdf10 100644 --- a/code/spectral_mech_Basic.f90 +++ b/code/spectral_mech_Basic.f90 @@ -48,7 +48,7 @@ module spectral_mech_basic C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness - S = 0.0_pReal !< current compliance (filled up with zeros) + S = 0.0_pReal !< current compliance (filled up with zeros) real(pReal), private :: err_stress, err_div logical, private :: ForwardData integer(pInt), private :: & @@ -61,21 +61,7 @@ module spectral_mech_basic BasicPETSc_forward, & basicPETSc_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASNESSetFunctionLocal, & PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & MPI_Abort, & MPI_Bcast, & MPI_Allreduce @@ -105,7 +91,7 @@ subroutine basicPETSc_init use spectral_utilities, only: & Utilities_constitutiveResponse, & Utilities_updateGamma, & - utilities_updateIPcoords, & + Utilities_updateIPcoords, & wgt use mesh, only: & grid, & @@ -115,15 +101,28 @@ subroutine basicPETSc_init implicit none real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P - PetscScalar, dimension(:,:,:,:), pointer :: F - PetscErrorCode :: ierr - PetscObject :: dummy real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal + + PetscErrorCode :: ierr + PetscObject :: dummy + PetscScalar, pointer, dimension(:,:,:,:) :: F + integer(pInt), dimension(:), allocatable :: localK integer(pInt) :: proc character(len=1024) :: rankStr - + + external :: & + SNESCreate, & + SNESSetOptionsPrefix, & + DMDACreate3D, & + SNESSetDM, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions + mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -147,9 +146,9 @@ subroutine basicPETSc_init DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid - 1, 1, worldsize, & + 1 , 1, worldsize, & 9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - grid (1),grid (2),localK, & ! local grid + grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) @@ -195,10 +194,9 @@ subroutine basicPETSc_init temp33_Real, & .false., & math_I3) - call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc - restartRead: if (restartInc > 1_pInt) then + restartRead: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' @@ -243,19 +241,24 @@ type(tSolutionState) function & timeinc, & !< increment in time for current solution timeinc_old, & !< increment in time of last increment loadCaseTime !< remaining time of current load case + logical, intent(in) :: & + guess type(tBoundaryCondition), intent(in) :: & P_BC, & F_BC character(len=*), intent(in) :: & incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC - logical, intent(in) :: & - guess !-------------------------------------------------------------------------------------------------- ! PETSc Data PetscErrorCode :: ierr SNESConvergedReason :: reason + + external :: & + SNESSolve, & + SNESGetConvergedReason + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- @@ -263,9 +266,9 @@ type(tSolutionState) function & S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite) - + !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data mask_stress = P_BC%maskFloat params%P_BC = P_BC%values params%rotation_BC = rotation_BC @@ -292,7 +295,7 @@ end function BasicPETSc_solution !-------------------------------------------------------------------------------------------------- -!> @brief forms the AL residual vector +!> @brief forms the basic residual vector !-------------------------------------------------------------------------------------------------- subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) use numerics, only: & @@ -312,10 +315,11 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) debug_spectral, & debug_spectralRotation use spectral_utilities, only: & + wgt, & tensorField_real, & utilities_FFTtensorForward, & - utilities_FFTtensorBackward, & utilities_fourierGammaConvolution, & + utilities_FFTtensorBackward, & Utilities_constitutiveResponse, & Utilities_divergenceRMS use IO, only: & @@ -338,11 +342,15 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr + external :: & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment - newIteration: if (totalIter <= PETScIter) then + newIteration: if(totalIter <= PETScIter) then !-------------------------------------------------------------------------------------------------- ! report begin of new iteration totalIter = totalIter + 1_pInt @@ -351,7 +359,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr) ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim (lab) =', & - math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) + math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC)) write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', & math_transpose33(F_aim) flush(6) @@ -401,7 +409,7 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du worldrank use FEsolving, only: & terminallyIll - + implicit none SNES :: snes_local PetscInt :: PETScIter @@ -415,10 +423,10 @@ subroutine BasicPETSc_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,du real(pReal) :: & divTol, & stressTol - + divTol = max(maxval(abs(P_av))*err_div_tolRel,err_div_tolAbs) stressTol = max(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) - + converged: if ((totalIter >= itmin .and. & all([ err_div/divTol, & err_stress/stressTol ] < 1.0_pReal)) & @@ -451,21 +459,21 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r use math, only: & math_mul33x33 ,& math_rotate_backward33 + use numerics, only: & + worldrank use mesh, only: & grid, & grid3 use spectral_utilities, only: & Utilities_calculateRate, & Utilities_forwardField, & - utilities_updateIPcoords, & + Utilities_updateIPcoords, & tBoundaryCondition, & cutBack use IO, only: & IO_write_JobRealFile use FEsolving, only: & restartWrite - use numerics, only: & - worldrank implicit none real(pReal), intent(in) :: & @@ -478,8 +486,9 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r real(pReal), dimension(3,3), intent(in) :: rotation_BC logical, intent(in) :: & guess + PetscErrorCode :: ierr PetscScalar, pointer :: F(:,:,:,:) - PetscErrorCode :: ierr + character(len=1024) :: rankStr call DMDAVecGetArrayF90(da,solution_vec,F,ierr) @@ -508,7 +517,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r write (777,rec=1) C_volAvgLastInc close(777) endif - endif + endif call utilities_updateIPcoords(F) @@ -538,6 +547,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3])) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) endif + F_aim = F_aim + f_aimDot * timeinc !-------------------------------------------------------------------------------------------------- @@ -558,6 +568,11 @@ subroutine BasicPETSc_destroy() implicit none PetscErrorCode :: ierr + external :: & + VecDestroy, & + SNESDestroy, & + DMDestroy + call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) diff --git a/code/spectral_mech_Polarisation.f90 b/code/spectral_mech_Polarisation.f90 index d7f1599e5..fc9fae7d7 100644 --- a/code/spectral_mech_Polarisation.f90 +++ b/code/spectral_mech_Polarisation.f90 @@ -22,7 +22,7 @@ module spectral_mech_Polarisation DAMASK_spectral_solverPolarisation_label = 'polarisation' !-------------------------------------------------------------------------------------------------- -! derived types +! derived types type(tSolutionParams), private :: params real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal @@ -31,7 +31,7 @@ module spectral_mech_Polarisation DM, private :: da SNES, private :: snes Vec, private :: solution_vec - + !-------------------------------------------------------------------------------------------------- ! common pointwise data real(pReal), private, dimension(:,:,:,:,:), allocatable :: & @@ -57,7 +57,7 @@ module spectral_mech_Polarisation S = 0.0_pReal, & !< current compliance (filled up with zeros) C_scale = 0.0_pReal, & S_scale = 0.0_pReal - + real(pReal), private :: & err_BC, & !< deviation from stress BC err_curl, & !< RMS of curl of F @@ -72,21 +72,7 @@ module spectral_mech_Polarisation Polarisation_forward, & Polarisation_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASNESSetFunctionLocal, & PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & MPI_Abort, & MPI_Bcast, & MPI_Allreduce @@ -136,11 +122,22 @@ subroutine Polarisation_init integer(pInt) :: proc character(len=1024) :: rankStr - if (worldrank == 0_pInt) then + external :: & + SNESCreate, & + SNESSetOptionsPrefix, & + DMDACreate3D, & + SNESSetDM, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + SNESGetConvergedReason, & + SNESSetConvergenceTest, & + SNESSetFromOptions + + mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - endif + endif mainProcess !-------------------------------------------------------------------------------------------------- ! allocate global fields @@ -150,7 +147,7 @@ subroutine Polarisation_init allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- -! PETSc Init +! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 @@ -163,7 +160,7 @@ subroutine Polarisation_init grid(1),grid(2),grid(3), & ! global grid 1 , 1, worldsize, & 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) - grid (1),grid (2),localK, & ! local grid + grid(1),grid(2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) @@ -182,7 +179,7 @@ subroutine Polarisation_init restart: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & - 'reading values of increment', restartInc - 1_pInt, 'from file' + 'reading values of increment ', restartInc - 1_pInt, ' from file' flush(6) write(rankStr,'(a1,i0)')'_',worldrank call IO_read_realFile(777,'F'//trim(rankStr),trim(getSolverJobName()),size(F)) @@ -221,7 +218,7 @@ subroutine Polarisation_init nullify(F_tau) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc - readRestart: if (restartInc > 1_pInt) then + restartRead: if (restartInc > 1_pInt) then if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) & write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') & 'reading more values of increment', restartInc - 1_pInt, 'from file' @@ -235,7 +232,7 @@ subroutine Polarisation_init call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg)) read (777,rec=1) C_minMaxAvg close (777) - endif readRestart + endif restartRead call Utilities_updateGamma(C_minMaxAvg,.True.) C_scale = C_minMaxAvg @@ -262,7 +259,7 @@ type(tSolutionState) function & use FEsolving, only: & restartWrite, & terminallyIll - + implicit none !-------------------------------------------------------------------------------------------------- @@ -285,6 +282,10 @@ type(tSolutionState) function & PetscErrorCode :: ierr SNESConvergedReason :: reason + external :: & + SNESSolve, & + SNESGetConvergedReason + incInfo = incInfoIn !-------------------------------------------------------------------------------------------------- @@ -385,7 +386,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) PetscErrorCode :: ierr integer(pInt) :: & i, j, k, e - + + external :: & + SNESGetNumberFunctionEvals, & + SNESGetIterationNumber + F => x_scal(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) F_tau => x_scal(1:3,1:3,2,& @@ -505,7 +510,7 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, fnorm SNESConvergedReason :: reason PetscObject :: dummy - PetscErrorCode ::ierr + PetscErrorCode :: ierr real(pReal) :: & curlTol, & divTol, & @@ -631,7 +636,8 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC write (777,rec=1) C_volAvgLastInc close(777) endif - endif + endif + call utilities_updateIPcoords(F) if (cutBack) then @@ -701,6 +707,11 @@ subroutine Polarisation_destroy() implicit none PetscErrorCode :: ierr + external :: & + VecDestroy, & + SNESDestroy, & + DMDestroy + call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) diff --git a/code/spectral_thermal.f90 b/code/spectral_thermal.f90 index 985bdf0c6..ab980e091 100644 --- a/code/spectral_thermal.f90 +++ b/code/spectral_thermal.f90 @@ -42,7 +42,6 @@ module spectral_thermal integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment real(pReal), dimension(3,3), private :: D_ref real(pReal), private :: mobility_ref - character(len=1024), private :: incInfo public :: & spectral_thermal_init, & @@ -50,21 +49,7 @@ module spectral_thermal spectral_thermal_forward, & spectral_thermal_destroy external :: & - VecDestroy, & - DMDestroy, & - DMDACreate3D, & - DMCreateGlobalVector, & - DMDASNESSetFunctionLocal, & PETScFinalize, & - SNESDestroy, & - SNESGetNumberFunctionEvals, & - SNESGetIterationNumber, & - SNESSolve, & - SNESSetDM, & - SNESGetConvergedReason, & - SNESSetConvergenceTest, & - SNESSetFromOptions, & - SNESCreate, & MPI_Abort, & MPI_Bcast, & MPI_Allreduce @@ -99,10 +84,20 @@ subroutine spectral_thermal_init integer(pInt) :: proc integer(pInt) :: i, j, k, cell DM :: thermal_grid - PetscScalar, pointer :: x_scal(:,:,:) + PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr PetscObject :: dummy + external :: & + SNESCreate, & + SNESSetOptionsPrefix, & + DMDACreate3D, & + SNESSetDM, & + DMDAGetCorners, & + DMCreateGlobalVector, & + DMDASNESSetFunctionLocal, & + SNESSetFromOptions + mainProcess: if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -154,6 +149,8 @@ subroutine spectral_thermal_init x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) +!-------------------------------------------------------------------------------------------------- +! thermal reference diffusion update cell = 0_pInt D_ref = 0.0_pReal mobility_ref = 0.0_pReal @@ -171,7 +168,7 @@ subroutine spectral_thermal_init end subroutine spectral_thermal_init !-------------------------------------------------------------------------------------------------- -!> @brief solution for the Basic PETSC scheme with internal iterations +!> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_old,loadCaseTime) use numerics, only: & @@ -196,12 +193,18 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol integer(pInt) :: i, j, k, cell PetscInt :: position PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm - + !-------------------------------------------------------------------------------------------------- ! PETSc Data PetscErrorCode :: ierr SNESConvergedReason :: reason - + + external :: & + VecMin, & + VecMax, & + SNESSolve, & + SNESGetConvergedReason + spectral_thermal_solution%converged =.false. !-------------------------------------------------------------------------------------------------- @@ -355,8 +358,11 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime) logical, intent(in) :: guess integer(pInt) :: i, j, k, cell DM :: dm_local - PetscScalar, pointer :: x_scal(:,:,:) + PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr + + external :: & + SNESGetDM if (cutBack) then temperature_current = temperature_lastInc @@ -405,6 +411,10 @@ subroutine spectral_thermal_destroy() implicit none PetscErrorCode :: ierr + external :: & + VecDestroy, & + SNESDestroy + call VecDestroy(solution,ierr); CHKERRQ(ierr) call SNESDestroy(thermal_snes,ierr); CHKERRQ(ierr)