diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 4c5f9dd6d..f8836e401 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -13,8 +13,7 @@ module DAMASK_spectral_solverAL use math, only: & math_I3 use DAMASK_spectral_utilities, only: & - tSolutionState, & - tSolutionParams + tSolutionState implicit none private @@ -25,6 +24,14 @@ module DAMASK_spectral_solverAL !-------------------------------------------------------------------------------------------------- ! derived types + type tSolutionParams !< @todo use here the type definition for a full loadcase including mask + real(pReal), dimension(3,3) :: P_BC, rotation_BC + real(pReal) :: timeinc + real(pReal) :: timeincOld + real(pReal) :: temperature + real(pReal) :: density + end type tSolutionParams + type(tSolutionParams), private :: params real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal @@ -89,7 +96,9 @@ module DAMASK_spectral_solverAL SNESSetConvergenceTest, & SNESSetFromOptions, & SNESCreate, & - MPI_Abort + MPI_Abort, & + MPI_Bcast, & + MPI_Allreduce contains @@ -115,9 +124,12 @@ subroutine AL_init(temperature) use DAMASK_interface, only: & getSolverJobName use DAMASK_spectral_Utilities, only: & + Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords + Utilities_updateIPcoords, & + grid1Red, & + wgt use mesh, only: & gridLocal, & gridGlobal @@ -138,6 +150,7 @@ subroutine AL_init(temperature) integer(pInt) :: proc character(len=1024) :: rankStr + call Utilities_init() if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(a)') ' $Id$' @@ -156,7 +169,6 @@ subroutine AL_init(temperature) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) @@ -170,7 +182,6 @@ subroutine AL_init(temperature) gridLocal (1),gridLocal (2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,ierr) CHKERRQ(ierr) @@ -255,7 +266,7 @@ end subroutine AL_init !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, & - rotation_BC) + rotation_BC,density) use numerics, only: & update_gamma use math, only: & @@ -276,7 +287,8 @@ 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 - temperature_bc + temperature_bc, & + density logical, intent(in) :: & guess type(tBoundaryCondition), intent(in) :: & @@ -312,6 +324,7 @@ type(tSolutionState) function & params%timeinc = timeinc params%timeincOld = timeinc_old params%temperature = temperature_bc + params%density = density !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -350,14 +363,16 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33, & math_mul3333xx33, & math_invSym3333, & - math_mul33x33 + math_mul33x33, & + PI use DAMASK_spectral_Utilities, only: & wgt, & - tensorField_realMPI, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & + field_realMPI, & + field_fourierMPI, & + Utilities_FFTforward, & + Utilities_fourierConvolution, & Utilities_inverseLaplace, & - utilities_FFTtensorBackward, & + Utilities_FFTbackward, & Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS @@ -429,9 +444,9 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! - tensorField_realMPI = 0.0_pReal + field_realMPI = 0.0_pReal do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) - tensorField_realMPI(1:3,1:3,i,j,k) = & + field_realMPI(1:3,1:3,i,j,k) = & polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3)) @@ -440,13 +455,13 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call utilities_FFTtensorForward() - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call utilities_FFTtensorBackward() + call Utilities_FFTforward() + call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_lambda = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + residual_F_lambda = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -458,11 +473,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculate divergence - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F - call utilities_FFTtensorForward() + field_realMPI = 0.0_pReal + field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F + call Utilities_FFTforward() err_div = Utilities_divergenceRMS() - call utilities_FFTtensorBackward() + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -478,11 +493,11 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F - call utilities_FFTtensorForward() + field_realMPI = 0.0_pReal + field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + call Utilities_FFTforward() err_curl = Utilities_curlRMS() - call utilities_FFTtensorBackward() + call Utilities_FFTbackward() end subroutine AL_formResidual @@ -712,6 +727,7 @@ subroutine AL_destroy() call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) + call Utilities_destroy() end subroutine AL_destroy diff --git a/code/DAMASK_spectral_solverPolarisation.f90 b/code/DAMASK_spectral_solverPolarisation.f90 index 874dea863..fc65e7053 100644 --- a/code/DAMASK_spectral_solverPolarisation.f90 +++ b/code/DAMASK_spectral_solverPolarisation.f90 @@ -13,8 +13,7 @@ module DAMASK_spectral_solverPolarisation use math, only: & math_I3 use DAMASK_spectral_utilities, only: & - tSolutionState, & - tSolutionParams + tSolutionState implicit none private @@ -25,6 +24,14 @@ module DAMASK_spectral_solverPolarisation !-------------------------------------------------------------------------------------------------- ! derived types + type tSolutionParams !< @todo use here the type definition for a full loadcase including mask + real(pReal), dimension(3,3) :: P_BC, rotation_BC + real(pReal) :: timeinc + real(pReal) :: timeincOld + real(pReal) :: temperature + real(pReal) :: density + end type tSolutionParams + type(tSolutionParams), private :: params real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal @@ -117,9 +124,12 @@ subroutine Polarisation_init(temperature) use DAMASK_interface, only: & getSolverJobName use DAMASK_spectral_Utilities, only: & + Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - Utilities_updateIPcoords + Utilities_updateIPcoords, & + grid1Red, & + wgt use mesh, only: & gridLocal, & gridGlobal @@ -140,6 +150,7 @@ subroutine Polarisation_init(temperature) integer(pInt) :: proc character(len=1024) :: rankStr + call Utilities_init() if (worldrank == 0_pInt) then write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' write(6,'(a)') ' $Id$' @@ -158,7 +169,6 @@ subroutine Polarisation_init(temperature) !-------------------------------------------------------------------------------------------------- ! PETSc Init call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3) do proc = 1, worldsize call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) @@ -172,10 +182,10 @@ subroutine Polarisation_init(temperature) gridLocal (1),gridLocal (2),localK, & ! local grid da,ierr) ! handle, error CHKERRQ(ierr) - call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) call DMDASNESSetFunctionLocal(da,INSERT_VALUES,Polarisation_formResidual,dummy,ierr) CHKERRQ(ierr) + call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) call SNESSetConvergenceTest(snes,Polarisation_converged,dummy,PETSC_NULL_FUNCTION,ierr) CHKERRQ(ierr) call SNESSetFromOptions(snes,ierr); CHKERRQ(ierr) @@ -255,7 +265,7 @@ end subroutine Polarisation_init !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & Polarisation_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc, & - rotation_BC) + rotation_BC,density) use numerics, only: & update_gamma use math, only: & @@ -267,6 +277,8 @@ type(tSolutionState) function & use FEsolving, only: & restartWrite, & terminallyIll + use numerics, only: & + worldrank implicit none @@ -276,7 +288,8 @@ 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 - temperature_bc + temperature_bc, & + density logical, intent(in) :: & guess type(tBoundaryCondition), intent(in) :: & @@ -311,6 +324,7 @@ type(tSolutionState) function & params%timeinc = timeinc params%timeincOld = timeinc_old params%temperature = temperature_bc + params%density = density !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -349,14 +363,16 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) math_transpose33, & math_mul3333xx33, & math_invSym3333, & - math_mul33x33 + math_mul33x33, & + PI use DAMASK_spectral_Utilities, only: & wgt, & - tensorField_realMPI, & - utilities_FFTtensorForward, & - utilities_fourierGammaConvolution, & + field_realMPI, & + field_fourierMPI, & + Utilities_FFTforward, & + Utilities_fourierConvolution, & Utilities_inverseLaplace, & - utilities_FFTtensorBackward, & + Utilities_FFTbackward, & Utilities_constitutiveResponse, & Utilities_divergenceRMS, & Utilities_curlRMS @@ -428,9 +444,9 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! - tensorField_realMPI = 0.0_pReal + field_realMPI = 0.0_pReal do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1) - tensorField_realMPI(1:3,1:3,i,j,k) = & + field_realMPI(1:3,1:3,i,j,k) = & polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarAlpha*math_mul33x33(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)) @@ -438,13 +454,13 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! doing convolution in Fourier space - call utilities_FFTtensorForward() - call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) - call utilities_FFTtensorBackward() + call Utilities_FFTforward() + call Utilities_fourierConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual - residual_F_tau = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) + residual_F_tau = polarBeta*F - field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response @@ -456,11 +472,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculate divergence - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F - call utilities_FFTtensorForward() + field_realMPI = 0.0_pReal + field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F + call Utilities_FFTforward() err_div = Utilities_divergenceRMS() - call utilities_FFTtensorBackward() + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -476,11 +492,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr) !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_realMPI = 0.0_pReal - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F - call utilities_FFTtensorForward() + field_realMPI = 0.0_pReal + field_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F + call Utilities_FFTforward() err_curl = Utilities_curlRMS() - call utilities_FFTtensorBackward() + call Utilities_FFTbackward() end subroutine Polarisation_formResidual @@ -711,6 +727,7 @@ subroutine Polarisation_destroy() call VecDestroy(solution_vec,ierr); CHKERRQ(ierr) call SNESDestroy(snes,ierr); CHKERRQ(ierr) call DMDestroy(da,ierr); CHKERRQ(ierr) + call Utilities_destroy() end subroutine Polarisation_destroy