diff --git a/VERSION b/VERSION index 6cc584f93..9c296336f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.0-290-g12ac5e3 +v2.0.0-292-g2ebc5ec diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index e999bf5dc..53c9c4423 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -139,6 +139,7 @@ program DAMASK_spectral integer(MPI_OFFSET_KIND) :: fileOffset integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize integer(pInt), parameter :: maxByteOut = 2147483647-4096 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742 + integer(pInt), parameter :: maxRealOut = maxByteOut/pReal integer(pLongInt), dimension(2) :: outputIndex PetscErrorCode :: ierr external :: & @@ -444,8 +445,8 @@ program DAMASK_spectral if (.not. appendToOutFile) then ! if not restarting, write 0th increment do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output - outputIndex=int([(i-1_pInt)*((maxByteOut/pReal)/materialpoint_sizeResults)+1_pInt, & - min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) + outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & + min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),& [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), & (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,& @@ -471,7 +472,7 @@ program DAMASK_spectral ! forwarding time timeIncOld = timeinc if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale - timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used + timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used else if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale @@ -653,8 +654,8 @@ program DAMASK_spectral call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr) if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_seek') do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1 ! slice the output of my process in chunks not exceeding the limit for one output - outputIndex=int([(i-1_pInt)*((maxByteOut/pReal)/materialpoint_sizeResults)+1_pInt, & - min(i*((maxByteOut/pReal)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) + outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, & + min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt) call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),& [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), & (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,& diff --git a/code/hydrogenflux_cahnhilliard.f90 b/code/hydrogenflux_cahnhilliard.f90 index 898e7ed8d..35168da2d 100644 --- a/code/hydrogenflux_cahnhilliard.f90 +++ b/code/hydrogenflux_cahnhilliard.f90 @@ -223,8 +223,7 @@ function hydrogenflux_cahnhilliard_getMobility33(ip,el) enddo hydrogenflux_cahnhilliard_getMobility33 = & - hydrogenflux_cahnhilliard_getMobility33/ & - homogenization_Ngrains(mesh_element(3,el)) + hydrogenflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function hydrogenflux_cahnhilliard_getMobility33 @@ -258,8 +257,7 @@ function hydrogenflux_cahnhilliard_getDiffusion33(ip,el) enddo hydrogenflux_cahnhilliard_getDiffusion33 = & - hydrogenflux_cahnhilliard_getDiffusion33/ & - homogenization_Ngrains(mesh_element(3,el)) + hydrogenflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function hydrogenflux_cahnhilliard_getDiffusion33 @@ -295,8 +293,7 @@ function hydrogenflux_cahnhilliard_getFormationEnergy(ip,el) enddo hydrogenflux_cahnhilliard_getFormationEnergy = & - hydrogenflux_cahnhilliard_getFormationEnergy/ & - homogenization_Ngrains(mesh_element(3,el)) + hydrogenflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal) end function hydrogenflux_cahnhilliard_getFormationEnergy @@ -331,10 +328,9 @@ function hydrogenflux_cahnhilliard_getEntropicCoeff(ip,el) lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el)) enddo - hydrogenflux_cahnhilliard_getEntropicCoeff = & - hydrogenflux_cahnhilliard_getEntropicCoeff* & + hydrogenflux_cahnhilliard_getEntropicCoeff = hydrogenflux_cahnhilliard_getEntropicCoeff* & temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ & - homogenization_Ngrains(material_homog(ip,el)) + real(homogenization_Ngrains(material_homog(ip,el)),pReal) end function hydrogenflux_cahnhilliard_getEntropicCoeff @@ -393,8 +389,8 @@ subroutine hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_d enddo enddo - KPot = KPot/homogenization_Ngrains(material_homog(ip,el)) - dKPot_dCh = dKPot_dCh/homogenization_Ngrains(material_homog(ip,el)) + KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal) + dKPot_dCh = dKPot_dCh/real(homogenization_Ngrains(material_homog(ip,el)),pReal) end subroutine hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent 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) diff --git a/code/spectral_utilities.f90 b/code/spectral_utilities.f90 index 1bf322843..b9ed5245d 100644 --- a/code/spectral_utilities.f90 +++ b/code/spectral_utilities.f90 @@ -422,7 +422,7 @@ subroutine utilities_updateGamma(C,saveReference) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k-grid3Offset))*xi1st(m,i,j,k-grid3Offset) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_complex(l,m) = sum(C_ref(l,1:3,m,1:3)*xiDyad_cmplx) + temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex) matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex) if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then @@ -558,7 +558,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & xiDyad_cmplx(l,m) = conjg(-xi1st(l,i,j,k))*xi1st(m,i,j,k) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_complex(l,m) = sum(C_ref(l,1:3,m,1:3)*xiDyad_cmplx) + temp33_complex(l,m) = sum(cmplx(C_ref(l,1:3,m,1:3),0.0_pReal)*xiDyad_cmplx) matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex) matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex) if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then @@ -610,8 +610,8 @@ subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) ! do the actual spectral method calculation do k = 1_pInt, grid3; do j = 1_pInt, grid(2) ;do i = 1_pInt, grid1Red GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal)/ & - (cmplx(mobility_ref,0.0_pReal,pReal) + & - deltaT*sum(conjg(xi1st(1:3,i,j,k))*matmul(D_ref,xi1st(1:3,i,j,k)))) ! why not use dot_product + (cmplx(mobility_ref,0.0_pReal,pReal) + cmplx(deltaT,0.0_pReal)*& + sum(conjg(xi1st(1:3,i,j,k))* matmul(cmplx(D_ref,0.0_pReal),xi1st(1:3,i,j,k)))) ! why not use dot_product scalarField_fourier(i,j,k) = scalarField_fourier(i,j,k)*GreenOp_hat enddo; enddo; enddo