consistent handling of boundary conditions

A stress boundary condition 'P' indicates the stress at the end of the
load case (same as for 'F') 'Pdot' for given increase of (technical)
strain is not implemented.
Does not change anything for the most common case of zero-stress
boundary conditions, but simplifies the specification of stress ramps

src/grid/DAMASK_grid.f90

@@ -107,8 +107,8 @@ program DAMASK_grid
 !--------------------------------------------------------------------------------------------------
 ! initialize field solver information
   nActiveFields = 1
-  if (any(thermal_type  == THERMAL_conduction_ID  )) nActiveFields = nActiveFields + 1
+  if(any(thermal_type == THERMAL_conduction_ID)) nActiveFields = nActiveFields + 1
-  if (any(damage_type   == DAMAGE_nonlocal_ID     )) nActiveFields = nActiveFields + 1
+  if(any(damage_type  == DAMAGE_nonlocal_ID   )) nActiveFields = nActiveFields + 1
   allocate(solres(nActiveFields))
   allocate(newLoadCase%ID(nActiveFields))
 
@@ -181,7 +181,7 @@ program DAMASK_grid
     if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) &                                            ! sanity check
       call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile))             ! error message for incomplete loadcase
   
-    newLoadCase%stress%myType='stress'
+    newLoadCase%stress%myType='p'
     field = 1
     newLoadCase%ID(field) = FIELD_MECH_ID                                                           ! mechanical active by default
     thermalActive: if (any(thermal_type  == THERMAL_conduction_ID)) then

src/grid/grid_mech_FEM.f90

@@ -66,8 +66,8 @@ module grid_mech_FEM
     F_aim = math_I3, &                                                                              !< current prescribed deformation gradient
     F_aim_lastIter = math_I3, &
     F_aim_lastInc  = math_I3, &                                                                     !< previous average deformation gradient
-    P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
+    P_av = 0.0_pReal, &                                                                             !< average 1st Piola--Kirchhoff stress
-
+    P_aim = 0.0_pReal
   character(len=:), allocatable :: incInfo                                                          !< time and increment information
   real(pReal), dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness
@@ -327,7 +327,6 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
 !--------------------------------------------------------------------------------------------------
 ! set module wide available data
   params%stress_mask = stress_BC%maskFloat
-  params%stress_BC   = stress_BC%values
   params%rotation_BC = rotation_BC
   params%timeinc     = timeinc
   params%timeincOld  = timeinc_old
@@ -374,6 +373,12 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
+  if     (stress_BC%myType=='p') then
+    P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc
+  elseif (stress_BC%myType=='pdot') then !UNTESTED
+    P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc
+  endif
+
   call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
 
   call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
@@ -489,8 +494,6 @@ subroutine formResidual(da_local,x_local, &
   PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
   PetscScalar, dimension(8,3) :: x_elem,  f_elem
   PetscInt             :: i, ii, j, jj, k, kk, ctr, ele
-  real(pReal), dimension(3,3) :: &
-    deltaF_aim
   PetscInt :: &
     PETScIter, &
     nfuncs
@@ -539,10 +542,8 @@ subroutine formResidual(da_local,x_local, &
 
 !--------------------------------------------------------------------------------------------------
 ! stress BC handling
-  F_aim_lastIter = F_aim
+  F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim)                                                 ! S = 0.0 for no bc
-  deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)
+  err_BC = maxval(abs(params%stress_mask * (P_av - P_aim)))                                         ! mask = 0.0 when no stress bc
-  F_aim = F_aim - deltaF_aim
-  err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC)))                              ! mask = 0.0 when no stress bc
 
 !--------------------------------------------------------------------------------------------------
 ! constructing residual

src/grid/grid_mech_spectral_basic.f90

@@ -62,8 +62,8 @@ module grid_mech_spectral_basic
     F_aimDot = 0.0_pReal, &                                                                         !< assumed rate of average deformation gradient
     F_aim = math_I3, &                                                                              !< current prescribed deformation gradient
     F_aim_lastInc = math_I3, &                                                                      !< previous average deformation gradient
-    P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
+    P_av = 0.0_pReal, &                                                                             !< average 1st Piola--Kirchhoff stress
-
+    P_aim = 0.0_pReal
   character(len=:), allocatable :: incInfo                                                          !< time and increment information
   real(pReal), dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness
@@ -129,7 +129,7 @@ subroutine grid_mech_spectral_basic_init
   num%eps_div_atol    = num_grid%get_asFloat  ('eps_div_atol',   defaultVal=1.0e-4_pReal)
   num%eps_div_rtol    = num_grid%get_asFloat  ('eps_div_rtol',   defaultVal=5.0e-4_pReal)
   num%eps_stress_atol = num_grid%get_asFloat  ('eps_stress_atol',defaultVal=1.0e3_pReal)
-  num%eps_stress_rtol = num_grid%get_asFloat  ('eps_stress_rtol',defaultVal=0.01_pReal)
+  num%eps_stress_rtol = num_grid%get_asFloat  ('eps_stress_rtol',defaultVal=0.001_pReal)
   num%itmin           = num_grid%get_asInt    ('itmin',defaultVal=1)
   num%itmax           = num_grid%get_asInt    ('itmax',defaultVal=250)
 
@@ -292,7 +292,6 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 !--------------------------------------------------------------------------------------------------
 ! set module wide available data
   params%stress_mask = stress_BC%maskFloat
-  params%stress_BC   = stress_BC%values
   params%rotation_BC = rotation_BC
   params%timeinc     = timeinc
   params%timeincOld  = timeinc_old
@@ -333,6 +332,12 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
+  if     (stress_BC%myType=='p') then
+    P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc
+  elseif (stress_BC%myType=='pdot') then !UNTESTED
+    P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc
+  endif
+ 
   F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, &                                       ! estimate of F at end of time+timeinc that matches rotated F_aim on average
               rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@@ -484,9 +489,9 @@ subroutine formResidual(in, F, &
 
 !--------------------------------------------------------------------------------------------------
 ! stress BC handling
-  deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC)                                         ! S = 0.0 for no bc
+  deltaF_aim = math_mul3333xx33(S, P_av - P_aim)                                                    ! S = 0.0 for no bc
-  F_aim = F_aim -deltaF_aim
+  F_aim = F_aim - deltaF_aim
-  err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC)))                              ! mask = 0.0 when no stress bc
+  err_BC = maxval(abs(params%stress_mask * (P_av - P_aim)))                                         ! mask = 0.0 when no stress bc
 
 !--------------------------------------------------------------------------------------------------
 ! updated deformation gradient using fix point algorithm of basic scheme
@@ -499,7 +504,7 @@ subroutine formResidual(in, F, &
 
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
-  residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)                                   ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
+  residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)                                  ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
 
 end subroutine formResidual
 

src/grid/grid_mech_spectral_polarisation.f90

@@ -70,8 +70,8 @@ module grid_mech_spectral_polarisation
     F_aim = math_I3, &                                                                              !< current prescribed deformation gradient
     F_aim_lastInc = math_I3, &                                                                      !< previous average deformation gradient
     F_av = 0.0_pReal, &                                                                             !< average incompatible def grad field
-    P_av = 0.0_pReal                                                                                !< average 1st Piola--Kirchhoff stress
+    P_av = 0.0_pReal, &                                                                             !< average 1st Piola--Kirchhoff stress
-
+    P_aim = 0.0_pReal
   character(len=:), allocatable :: incInfo                                                          !< time and increment information
   real(pReal), dimension(3,3,3,3) :: &
     C_volAvg = 0.0_pReal, &                                                                         !< current volume average stiffness
@@ -326,7 +326,6 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 !--------------------------------------------------------------------------------------------------
 ! set module wide available data
   params%stress_mask = stress_BC%maskFloat
-  params%stress_BC   = stress_BC%values
   params%rotation_BC = rotation_BC
   params%timeinc     = timeinc
   params%timeincOld  = timeinc_old
@@ -373,6 +372,12 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
+  if     (stress_BC%myType=='p') then
+    P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc
+  elseif (stress_BC%myType=='pdot') then !UNTESTED
+    P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc
+  endif
+
   F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! estimate of F at end of time+timeinc that matches rotated F_aim on average
                                      rotation_BC%rotate(F_aim,active=.true.)),&
               [9,grid(1),grid(2),grid3])
@@ -587,10 +592,10 @@ subroutine formResidual(in, FandF_tau, &
 
 !--------------------------------------------------------------------------------------------------
 ! stress BC handling
-  F_aim = F_aim - math_mul3333xx33(S, P_av - params%stress_BC)                                     ! S = 0.0 for no bc
+  F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim)                                                 ! S = 0.0 for no bc
   err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
                                                  -params%rotation_BC%rotate(F_av)) + &
-                                 params%stress_mask  * (P_av-params%stress_BC)))                    ! mask = 0.0 for no bc
+                                 params%stress_mask  * (P_av-P_aim)))                               ! mask = 0.0 for no bc
 ! calculate divergence
   tensorField_real = 0.0_pReal
   tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F                                !< stress field in disguise