From ce7a0571fd75a5b3f897a27d13a1514e24d32c5c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Jan 2013 10:12:03 +0000 Subject: [PATCH] fixed bug in forwarding fields for AL solver --- code/DAMASK_spectral_solverAL.f90 | 11 +++++----- code/DAMASK_spectral_solverBasic.f90 | 2 +- code/DAMASK_spectral_solverBasicPETSc.f90 | 4 ++-- code/DAMASK_spectral_utilities.f90 | 25 ++++++++++++++--------- 4 files changed, 24 insertions(+), 18 deletions(-) diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index cc2ed44d8..818de7798 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -46,7 +46,8 @@ module DAMASK_spectral_solverAL F_lambda_lastInc, & !< field of previous incompatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient F_lambdaDot !< field of assumed rate of incopatible deformation gradient - real(pReal), private :: temperature !< temperature, no spatial quantity at the moment + real(pReal), private :: & + temperature !< temperature, no spatial quantity at the moment !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. @@ -325,10 +326,10 @@ type(tSolutionState) function & !-------------------------------------------------------------------------------------------------- ! update local deformation gradient - F = reshape(Utilities_forwardField(timeinc,math_rotate_backward33(F_aim,rotation_BC),& - F_lastInc,Fdot),[9,res(1),res(2),res(3)]) - F_lambda = reshape(Utilities_forwardField(timeinc,math_rotate_backward33(F_aim,rotation_BC),& - F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)]) + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim + math_rotate_backward33(F_aim,rotation_BC)),[9,res(1),res(2),res(3)]) + F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), & ! does not have any average value as boundary condition + [9,res(1),res(2),res(3)]) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr) CHKERRQ(ierr) diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index 8a286fb42..a8768a921 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -261,7 +261,7 @@ type(tSolutionState) function & F_lastInc = F endif F_aim = F_aim + f_aimDot * timeinc - F = Utilities_forwardField(timeinc,math_rotate_backward33(F_aim,rotation_BC),F_lastInc,Fdot) + F = Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim,rotation_BC)) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 57ea8c595..d90b1e8e7 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -292,8 +292,8 @@ type(tSolutionState) function & endif F_aim = F_aim + f_aimDot * timeinc - F = reshape(Utilities_forwardField(timeinc,math_rotate_backward33(F_aim,rotation_BC),F_lastInc,& - Fdot),[9,res(1),res(2),res(3)]) + F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot,math_rotate_backward33(F_aim, & + rotation_BC)),[9,res(1),res(2),res(3)]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) CHKERRQ(ierr) diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index a86d4d120..19b726b62 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -787,26 +787,31 @@ end function utilities_calculateRate !-------------------------------------------------------------------------------------------------- -!> @brief forwards a field with a pointwise given rate, ensures that the average matches the aim +!> @brief forwards a field with a pointwise given rate, if aim is given, +!> ensures that the average matches the aim !-------------------------------------------------------------------------------------------------- -pure function utilities_forwardField(timeinc,aim,field_lastInc,rate) +pure function utilities_forwardField(timeinc,field_lastInc,rate,aim) use mesh, only: & res, & wgt implicit none - real(pReal), intent(in) :: timeinc !< timeinc of current step - real(pReal), intent(in), dimension(3,3) :: aim !< average field value aim - real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & - field_lastInc,& !< initial field + real(pReal), intent(in) :: & + timeinc !< timeinc of current step + real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: & + field_lastInc, & !< initial field rate !< rate by which to forward - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_forwardField - real(pReal), dimension(3,3) :: fieldDiff !< - aim + real(pReal), intent(in), optional, dimension(3,3) :: & + aim !< average field value aim + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_forwardField + real(pReal), dimension(3,3) :: fieldDiff !< - aim utilities_forwardField = field_lastInc + rate*timeinc - fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim - utilities_forwardField = utilities_forwardField - & + if (present(aim)) then !< correct to match average + fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim + utilities_forwardField = utilities_forwardField - & spread(spread(spread(fieldDiff,3,res(1)),4,res(2)),5,res(3)) + endif end function utilities_forwardField