From 61981617d7a03f845a1a078fbbbc1fa8d241b93d Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Wed, 18 Dec 2013 09:35:05 +0000 Subject: [PATCH] separated forwarding and solution subroutines for better control at the load step looping level --- code/DAMASK_spectral_driver.f90 | 31 ++++ code/DAMASK_spectral_solverAL.f90 | 174 ++++++++++++-------- code/DAMASK_spectral_solverBasicPETSc.f90 | 130 +++++++++------ code/DAMASK_spectral_solverPolarisation.f90 | 172 +++++++++++-------- 4 files changed, 316 insertions(+), 191 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index bfa9accd5..f11d16b06 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -433,6 +433,37 @@ program DAMASK_spectral_Driver stepFraction = stepFraction + 1_pInt remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc !-------------------------------------------------------------------------------------------------- +! forward solution + select case(myspectralsolver) + case (DAMASK_spectral_SolverBasic_label) +#ifdef PETSc + case (DAMASK_spectral_SolverBasicPETSC_label) + call BasicPETSC_forward (& + guess,timeinc,timeIncOld,remainingLoadCaseTime, & + P_BC = loadCases(currentLoadCase)%P, & + F_BC = loadCases(currentLoadCase)%deformation, & + rotation_BC = loadCases(currentLoadCase)%rotation) + + case (DAMASK_spectral_SolverAL_label) + call AL_forward (& + guess,timeinc,timeIncOld,remainingLoadCaseTime, & + P_BC = loadCases(currentLoadCase)%P, & + F_BC = loadCases(currentLoadCase)%deformation, & + rotation_BC = loadCases(currentLoadCase)%rotation) + + case (DAMASK_spectral_SolverPolarisation_label) + call Polarisation_forward (& + guess,timeinc,timeIncOld,remainingLoadCaseTime, & + P_BC = loadCases(currentLoadCase)%P, & + F_BC = loadCases(currentLoadCase)%deformation, & + rotation_BC = loadCases(currentLoadCase)%rotation) +#endif + end select + + end select + enddo + +!-------------------------------------------------------------------------------------------------- ! report begin of new increment write(6,'(/,a)') ' ###########################################################################' write(6,'(1x,a,es12.5'//& diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 9b0cb5b71..ad541925f 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -100,6 +100,7 @@ module DAMASK_spectral_solverAL public :: & AL_init, & AL_solution, & + AL_forward, & AL_destroy external :: & VecDestroy, & @@ -277,25 +278,15 @@ type(tSolutionState) function & use numerics, only: & update_gamma use math, only: & - math_mul33x33 ,& - math_mul3333xx33, & - math_rotate_backward33, & - math_invSym3333, & - math_transpose33 -use mesh, only: & - mesh_ipCoordinates, & - mesh_deformedCoordsFFT + math_invSym3333 use IO, only: & IO_write_jobRealFile use DAMASK_spectral_Utilities, only: & grid, & geomSize, & tBoundaryCondition, & - Utilities_forwardField, & - Utilities_calculateRate, & Utilities_maskedCompliance, & - Utilities_updateGamma, & - cutBack + Utilities_updateGamma use FEsolving, only: & restartWrite, & terminallyIll @@ -327,16 +318,13 @@ use mesh, only: & PetscErrorCode :: ierr SNESConvergedReason :: reason - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: F_lambda33 - incInfo = incInfoIn - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart' flush(6) @@ -365,61 +353,10 @@ use mesh, only: & write (777,rec=1) C_volAvgLastInc close(777) endif + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + AL_solution%converged =.false. - if (cutBack) then - F_aim = F_aim_lastInc - F_lambda= reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid(3)]) - F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) - C_volAvg = C_volAvgLastInc - else - C_volAvgLastInc = C_volAvg -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime - endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old - F_aim_lastInc = F_aim - -!-------------------------------------------------------------------------------------------------- -! update coordinates and rate and forward last inc - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& - F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) - F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])) - F_lastInc2 = F_lastInc - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) - F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]) - endif - F_aim = F_aim + f_aimDot * timeinc - -!-------------------------------------------------------------------------------------------------- -! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim - math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)]) - F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & - [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition - if (.not. guess) then ! large strain forwarding - do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) - F_lambda33 = reshape(F_lambda(:,i,j,k),[3,3]) - F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & - math_mul3333xx33(C_scale,& - math_mul33x33(math_transpose33(F_lambda33),& - F_lambda33) -math_I3))*0.5_pReal)& - + math_I3 - F_lambda(:,i,j,k) = reshape(F_lambda33,[9]) - enddo; enddo; enddo - endif - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) - CHKERRQ(ierr) - !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) @@ -696,6 +633,101 @@ subroutine AL_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr end subroutine AL_converged +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) + use math, only: & + math_mul33x33, & + math_mul3333xx33, & + math_transpose33, & + math_rotate_backward33 + use DAMASK_spectral_Utilities, only: & + grid, & + geomSize, & + Utilities_calculateRate, & + Utilities_forwardField, & + tBoundaryCondition, & + cutBack + use mesh, only: & + mesh_ipCoordinates,& + mesh_deformedCoordsFFT + + implicit none +#include + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) + if (cutBack) then + F_aim = F_aim_lastInc + F_lambda= reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid(3)]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) + C_volAvg = C_volAvgLastInc + else +!-------------------------------------------------------------------------------------------------- +! calculate rate for aim + C_volAvgLastInc = C_volAvg + if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F + f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = F_BC%maskFloat * F_BC%values + elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + endif + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) + F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)])) + F_lastInc2 = F_lastInc + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) + F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid(3)]) + endif + + F_aim = F_aim + f_aimDot * timeinc + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)), & + [9,grid(1),grid(2),grid(3)]) + F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & + [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition + if (.not. guess) then ! large strain forwarding + do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + F_lambda33 = reshape(F_lambda(:,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & + math_mul3333xx33(C_scale,& + math_mul33x33(math_transpose33(F_lambda33),& + F_lambda33) -math_I3))*0.5_pReal)& + + math_I3 + F_lambda(:,i,j,k) = reshape(F_lambda33,[9]) + enddo; enddo; enddo + endif + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + +end subroutine AL_forward + !-------------------------------------------------------------------------------------------------- !> @brief destroy routine !-------------------------------------------------------------------------------------------------- diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 035b917d0..ca4d39e01 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -86,7 +86,8 @@ module DAMASK_spectral_SolverBasicPETSc public :: & basicPETSc_init, & - basicPETSc_solution ,& + basicPETSc_solution, & + BasicPETSc_forward, & basicPETSc_destroy external :: & VecDestroy, & @@ -249,23 +250,13 @@ type(tSolutionState) function basicPETSc_solution( & use numerics, only: & update_gamma, & itmax - use math, only: & - math_mul33x33 ,& - math_rotate_backward33 - use mesh, only: & - mesh_ipCoordinates,& - mesh_deformedCoordsFFT use IO, only: & IO_write_JobRealFile use DAMASK_spectral_Utilities, only: & grid, & - geomSize, & tBoundaryCondition, & - Utilities_forwardField, & - Utilities_calculateRate, & Utilities_maskedCompliance, & - Utilities_updateGamma, & - cutBack + Utilities_updateGamma use FEsolving, only: & restartWrite, & terminallyIll @@ -282,14 +273,14 @@ type(tSolutionState) function basicPETSc_solution( & loadCaseTime, & !< remaining time of current load case temperature_bc, & density - 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 @@ -298,9 +289,9 @@ type(tSolutionState) function basicPETSc_solution( & SNESConvergedReason :: reason incInfo = incInfoIn - call DMDAVecGetArrayF90(da,solution_vec,F,ierr) !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver + call DMDAVecGetArrayF90(da,solution_vec,F,ierr) if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart' flush(6) @@ -323,42 +314,9 @@ type(tSolutionState) function basicPETSc_solution( & write (777,rec=1) C_volAvgLastInc close(777) endif - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& - F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - if (cutBack) then - F_aim = F_aim_lastInc - F = reshape(F_lastInc,[9,grid(1),grid(2),grid(3)]) - C_volAvg = C_volAvgLastInc - else - C_volAvgLastInc = C_volAvg - - -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime - endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old - F_aim_lastInc = F_aim - -!-------------------------------------------------------------------------------------------------- -! update coordinates and rate and forward last inc - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& - F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,params%rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) - F_lastInc2 = F_lastInc - F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid(3)]) - endif - F_aim = F_aim + f_aimDot * timeinc - - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, & - rotation_BC)),[9,grid(1),grid(2),grid(3)]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + + BasicPETSc_solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) @@ -607,6 +565,78 @@ subroutine BasicPETSc_convergedKSP(ksp_local,PETScIter,fnorm,reason,dummy,ierr) end subroutine BasicPETSc_convergedKSP +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) + use math, only: & + math_mul33x33 ,& + math_rotate_backward33 + use DAMASK_spectral_Utilities, only: & + grid, & + geomSize, & + Utilities_calculateRate, & + Utilities_forwardField, & + tBoundaryCondition, & + cutBack + use mesh, only: & + mesh_ipCoordinates,& + mesh_deformedCoordsFFT + + implicit none +#include + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + PetscScalar, pointer :: F(:,:,:,:) + PetscErrorCode :: ierr + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call DMDAVecGetArrayF90(da,solution_vec,F,ierr) + if (cutBack) then + F_aim = F_aim_lastInc + F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) + C_volAvg = C_volAvgLastInc + else + C_volAvgLastInc = C_volAvg +!-------------------------------------------------------------------------------------------------- +! calculate rate for aim + if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F + f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = F_BC%maskFloat * F_BC%values + elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + endif + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) + F_lastInc2 = F_lastInc + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) + endif + F_aim = F_aim + f_aimDot * timeinc + +!-------------------------------------------------------------------------------------------------- +! update local deformation gradient + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)]) + call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) + +end subroutine BasicPETSc_forward !-------------------------------------------------------------------------------------------------- !> @brief destroy routine diff --git a/code/DAMASK_spectral_solverPolarisation.f90 b/code/DAMASK_spectral_solverPolarisation.f90 index 9feec95e5..98e150902 100644 --- a/code/DAMASK_spectral_solverPolarisation.f90 +++ b/code/DAMASK_spectral_solverPolarisation.f90 @@ -100,6 +100,7 @@ module DAMASK_spectral_solverPolarisation public :: & Polarisation_init, & Polarisation_solution, & + Polarisation_forward, & Polarisation_destroy external :: & VecDestroy, & @@ -277,25 +278,15 @@ type(tSolutionState) function & use numerics, only: & update_gamma use math, only: & - math_mul33x33 ,& - math_mul3333xx33, & - math_rotate_backward33, & - math_invSym3333, & - math_transpose33 -use mesh, only: & - mesh_ipCoordinates, & - mesh_deformedCoordsFFT + math_invSym3333 use IO, only: & IO_write_jobRealFile use DAMASK_spectral_Utilities, only: & grid, & geomSize, & tBoundaryCondition, & - Utilities_forwardField, & - Utilities_calculateRate, & Utilities_maskedCompliance, & - Utilities_updateGamma, & - cutBack + Utilities_updateGamma use FEsolving, only: & restartWrite, & terminallyIll @@ -327,16 +318,13 @@ use mesh, only: & PetscErrorCode :: ierr SNESConvergedReason :: reason - integer(pInt) :: i, j, k - real(pReal), dimension(3,3) :: F_lambda33 - incInfo = incInfoIn - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) - F => xx_psc(0:8,:,:,:) - F_tau => xx_psc(9:17,:,:,:) !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_tau => xx_psc(9:17,:,:,:) if (restartWrite) then write(6,'(/,a)') ' writing converged results for restart' flush(6) @@ -365,60 +353,9 @@ use mesh, only: & write (777,rec=1) C_volAvgLastInc close(777) endif + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) Polarisation_solution%converged =.false. - if (cutBack) then - F_aim = F_aim_lastInc - F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)]) - F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) - C_volAvg = C_volAvgLastInc - else - C_volAvgLastInc = C_volAvg -!-------------------------------------------------------------------------------------------------- -! calculate rate for aim - if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F - f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed - f_aimDot = F_BC%maskFloat * F_BC%values - elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed - f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime - endif - if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old - F_aim_lastInc = F_aim - -!-------------------------------------------------------------------------------------------------- -! update coordinates and rate and forward last inc - mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& - F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) - Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & - timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) - F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), & - timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])) - F_lastInc2 = F_lastInc - F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) - F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]) - endif - F_aim = F_aim + f_aimDot * timeinc - -!-------------------------------------------------------------------------------------------------- -! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim - math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)]) - F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition - if (.not. guess) then ! large strain forwarding - do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) - F_lambda33 = reshape(F_tau(:,i,j,k)-F(:,i,j,k),[3,3]) - F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & - math_mul3333xx33(C_scale,& - math_mul33x33(math_transpose33(F_lambda33),& - F_lambda33) -math_I3))*0.5_pReal)& - + math_I3 - F_tau(:,i,j,k) = reshape(F_lambda33,[9])+F(:,i,j,k) - enddo; enddo; enddo - endif - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) - CHKERRQ(ierr) - !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C_volAvg) @@ -694,6 +631,101 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason, end subroutine Polarisation_converged +!-------------------------------------------------------------------------------------------------- +!> @brief forwarding routine +!-------------------------------------------------------------------------------------------------- +subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_BC) + use math, only: & + math_mul33x33, & + math_mul3333xx33, & + math_transpose33, & + math_rotate_backward33 + use DAMASK_spectral_Utilities, only: & + grid, & + geomSize, & + Utilities_calculateRate, & + Utilities_forwardField, & + tBoundaryCondition, & + cutBack + use mesh, only: & + mesh_ipCoordinates,& + mesh_deformedCoordsFFT + + implicit none +#include + real(pReal), intent(in) :: & + timeinc_old, & + timeinc, & + loadCaseTime !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & + P_BC, & + F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + logical, intent(in) :: & + guess + PetscErrorCode :: ierr + PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_tau + integer(pInt) :: i, j, k + real(pReal), dimension(3,3) :: F_lambda33 + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) + F => xx_psc(0:8,:,:,:) + F_tau => xx_psc(9:17,:,:,:) + if (cutBack) then + F_aim = F_aim_lastInc + F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)]) + F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) + C_volAvg = C_volAvgLastInc + else + C_volAvgLastInc = C_volAvg +!-------------------------------------------------------------------------------------------------- +! calculate rate for aim + if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F + f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + f_aimDot = F_BC%maskFloat * F_BC%values + elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + endif + if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(& + F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)]) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)])) + F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), & + timeinc_old,guess,F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid(3)])) + F_lastInc2 = F_lastInc + F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid(3)]) + F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid(3)]) + endif + F_aim = F_aim + f_aimDot * timeinc + +!-------------------------------------------------------------------------------------------------- +! update local deformation gradient + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid(3)]) + F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), [9,grid(1),grid(2),grid(3)]) ! does not have any average value as boundary condition + if (.not. guess) then ! large strain forwarding + do k = 1_pInt, grid(3); do j = 1_pInt, grid(2); do i = 1_pInt, grid(1) + F_lambda33 = reshape(F_tau(:,i,j,k)-F(:,i,j,k),[3,3]) + F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & + math_mul3333xx33(C_scale,& + math_mul33x33(math_transpose33(F_lambda33),& + F_lambda33) -math_I3))*0.5_pReal)& + + math_I3 + F_tau(:,i,j,k) = reshape(F_lambda33,[9])+F(:,i,j,k) + enddo; enddo; enddo + endif + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) + +end subroutine Polarisation_forward + !-------------------------------------------------------------------------------------------------- !> @brief destroy routine !--------------------------------------------------------------------------------------------------