Merge remote-tracking branch 'origin/development' into misc-improvements

PRIVATE

@@ -1 +1 @@
-Subproject commit 6b8ba6d844b70695233b02eae85f4300315f339d
+Subproject commit e1a1048e1f593683b4b432d41455bd236008c3ad

VERSION

@@ -1 +1 @@
-v3.0.0-alpha2-26-gaad123f41
+v3.0.0-alpha2-35-g1ebd10745

python/damask/grid_filters.py

@@ -234,7 +234,7 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
     origin[_np.where(cells==1)] = 0.0
 
     if cells.prod() != len(coordinates0):
-        raise ValueError('Data count {len(coordinates0)} does not match cells {cells}.')
+        raise ValueError(f'Data count {len(coordinates0)} does not match cells {cells}.')
 
     start = origin + delta*.5
     end   = origin - delta*.5 + size
@@ -387,7 +387,7 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
     origin    = mincorner
 
     if (cells+1).prod() != len(coordinates0):
-        raise ValueError('Data count {len(coordinates0)} does not match cells {cells}.')
+        raise ValueError(f'Data count {len(coordinates0)} does not match cells {cells}.')
 
     atol = _np.max(size)*5e-2
     if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],cells[0]+1),atol=atol) and \

src/grid/DAMASK_grid.f90

@@ -42,8 +42,8 @@ program DAMASK_grid
 
 !--------------------------------------------------------------------------------------------------
 ! variables related to information from load case and geom file
-  real(pReal), dimension(9) :: temp_valueVector = 0.0_pReal                                         !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
+  real(pReal), dimension(9) :: temp_valueVector                                                     !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
-  logical,     dimension(9) :: temp_maskVector  = .false.                                           !< temporarily from loadcase file when reading in tensors
+  logical,     dimension(9) :: temp_maskVector                                                      !< temporarily from loadcase file when reading in tensors
 
 !--------------------------------------------------------------------------------------------------
 ! loop variables, convergence etc.
@@ -143,8 +143,6 @@ program DAMASK_grid
       mech_restartWrite => grid_mech_spectral_basic_restartWrite
 
     case ('Polarisation')
-     if(debug_grid%contains('basic')) &
-        call IO_warning(42, ext_msg='debug Divergence')
       mech_init         => grid_mech_spectral_polarisation_init
       mech_forward      => grid_mech_spectral_polarisation_forward
       mech_solution     => grid_mech_spectral_polarisation_solution
@@ -152,8 +150,6 @@ program DAMASK_grid
       mech_restartWrite => grid_mech_spectral_polarisation_restartWrite
 
     case ('FEM')
-     if(debug_grid%contains('basic')) &
-        call IO_warning(42, ext_msg='debug Divergence')
       mech_init         => grid_mech_FEM_init
       mech_forward      => grid_mech_FEM_forward
       mech_solution     => grid_mech_FEM_solution
@@ -190,33 +186,35 @@ program DAMASK_grid
     load_step => load_steps%get(l)
 
     step_mech => load_step%get('mechanics')
-    loadCases(l)%stress%myType='P'
+    loadCases(l)%stress%myType=''
     readMech: do m = 1, step_mech%length
       select case (step_mech%getKey(m))
-        case('dot_F','L','F')                                                                       ! assign values for the deformation BC matrix
+        case ('L','dot_F','F')                                                                      ! assign values for the deformation BC matrix
           loadCases(l)%deformation%myType = step_mech%getKey(m)
-          temp_valueVector = 0.0_pReal
-
           step_deformation => step_mech%get(m)
-          do j = 1, 9
+
-            temp_maskVector(j) = step_deformation%get_asString(j) /= 'x'                            ! true if not a 'x'
-            if (temp_maskVector(j)) temp_valueVector(j) = step_deformation%get_asFloat(j)           ! read value where applicable
-          enddo
-          loadCases(l)%deformation%mask   = transpose(reshape(temp_maskVector,[ 3,3]))              ! mask in 3x3 notation
-          loadCases(l)%deformation%values = math_9to33(temp_valueVector)                            ! values in 3x3 notation
-        case('P')
           temp_valueVector = 0.0_pReal
-          step_stress => step_mech%get(m)
           do j = 1, 9
-            temp_maskVector(j) = step_stress%get_asString(j) /= 'x'                                 ! true if not a 'x'
+            temp_maskVector(j) = step_deformation%get_asString(j) /= 'x'
-            if (temp_maskVector(j)) temp_valueVector(j) = step_stress%get_asFloat(j)                ! read value where applicable
+            if (temp_maskVector(j)) temp_valueVector(j) = step_deformation%get_asFloat(j)
           enddo
-          loadCases(l)%stress%mask   = transpose(reshape(temp_maskVector,[ 3,3]))
+          loadCases(l)%deformation%mask   = transpose(reshape(temp_maskVector,[3,3]))
+          loadCases(l)%deformation%values = math_9to33(temp_valueVector)
+        case ('dot_P','P')
+          loadCases(l)%stress%myType = step_mech%getKey(m)
+          step_stress => step_mech%get(m)
+
+          temp_valueVector = 0.0_pReal
+          do j = 1, 9
+            temp_maskVector(j) = step_stress%get_asString(j) /= 'x'
+            if (temp_maskVector(j)) temp_valueVector(j) = step_stress%get_asFloat(j)
+          enddo
+          loadCases(l)%stress%mask   = transpose(reshape(temp_maskVector,[3,3]))
           loadCases(l)%stress%values = math_9to33(temp_valueVector)
       end select
       call loadCases(l)%rot%fromAxisAngle(step_mech%get_asFloats('R',defaultVal = real([0.0,0.0,1.0,0.0],pReal)),degrees=.true.)
     enddo readMech
-    if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/F/dot_F missing')
+    if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing')
 
     step_discretization => load_step%get('discretization')
     if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
@@ -239,50 +237,60 @@ program DAMASK_grid
           if (any(loadCases(l)%deformation%mask(j,1:3) .eqv. .true.) .and. &
               any(loadCases(l)%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832                ! each row should be either fully or not at all defined
         enddo
-        print*, ' L:'
+      endif
-      else if (loadCases(l)%deformation%myType == 'F') then
+      if (loadCases(l)%deformation%myType == 'F') then
         print*, ' F:'
-      else if (loadCases(l)%deformation%myType == 'dot_F') then
+      else
-        print*, ' dot_F:'
+        print*, ' '//loadCases(l)%deformation%myType//' / 1/s:'
       endif
       do i = 1, 3; do j = 1, 3
-        if(loadCases(l)%deformation%mask(i,j)) then
+        if (loadCases(l)%deformation%mask(i,j)) then
           write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
         else
           write(IO_STDOUT,'(2x,12a)',advance='no') '     x      '
         endif
         enddo; write(IO_STDOUT,'(/)',advance='no')
       enddo
-      if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831          ! exclusive or masking only
+      if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831
       if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) &
         errorID = 838                                                                               ! no rotation is allowed by stress BC
-      print*, ' P / MPa:'
+
+      if (loadCases(l)%stress%myType == 'P')     print*, ' P / MPa:'
+      if (loadCases(l)%stress%myType == 'dot_P') print*, ' dot_P / MPa/s:'
+
+      if (loadCases(l)%stress%myType /= '') then
         do i = 1, 3; do j = 1, 3
-        if(loadCases(l)%stress%mask(i,j)) then
+          if (loadCases(l)%stress%mask(i,j)) then
             write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
           else
             write(IO_STDOUT,'(2x,12a)',advance='no') '     x      '
           endif
           enddo; write(IO_STDOUT,'(/)',advance='no')
         enddo
+      endif
       if (any(dNeq(loadCases(l)%rot%asMatrix(), math_I3))) &
         write(IO_STDOUT,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'R:',&
                  transpose(loadCases(l)%rot%asMatrix())
 
-      if (loadCases(l)%t < 0.0_pReal) errorID = 834
-      print'(a,f0.3)', '  t: ', loadCases(l)%t
-      if (loadCases(l)%N < 1)    errorID = 835
-      print'(a,i0)',   '  N: ', loadCases(l)%N
-      if (loadCases(l)%f_out < 1)  errorID = 836
-      print'(a,i0)',   '  f_out: ', loadCases(l)%f_out
       if (loadCases(l)%r <= 0.0)       errorID = 833
-      print'(a,f0.3)', '  r: ', loadCases(l)%r
+      if (loadCases(l)%t < 0.0_pReal)  errorID = 834
-
+      if (loadCases(l)%N < 1)          errorID = 835
+      if (loadCases(l)%f_out < 1)      errorID = 836
       if (loadCases(l)%f_restart < 1)  errorID = 839
+
+      if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then
+        print'(a)', '  r: 1 (constant step widths)'
+      else
+        print'(a,f0.3)', '  r: ', loadCases(l)%r
+      endif
+      print'(a,f0.3)', '  t: ', loadCases(l)%t
+      print'(a,i0)',   '  N: ', loadCases(l)%N
+      print'(a,i0)',   '  f_out: ', loadCases(l)%f_out
       if (loadCases(l)%f_restart < huge(0)) &
         print'(a,i0)', '  f_restart: ', loadCases(l)%f_restart
 
       if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string)                 ! exit with error message
+
     endif reportAndCheck
   enddo
 
@@ -309,8 +317,6 @@ program DAMASK_grid
     writeHeader: if (interface_restartInc < 1) then
       open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
       write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                ! statistics file
-      if (debug_grid%contains('basic')) print'(/,a)', ' header of statistics file written out'
-      flush(IO_STDOUT)
     else writeHeader
       open(newunit=statUnit,file=trim(getSolverJobName())//&
                                   '.sta',form='FORMATTED', position='APPEND', status='OLD')
@@ -319,6 +325,7 @@ program DAMASK_grid
 
   writeUndeformed: if (interface_restartInc < 1) then
     print'(/,a)', ' ... writing initial configuration to file ........................'
+    flush(IO_STDOUT)
     call CPFEM_results(0,0.0_pReal)
   endif writeUndeformed
 

src/grid/grid_mech_FEM.f90

@@ -93,10 +93,8 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_FEM_init
 
-  real(pReal) :: HGCoeff = 0.0e-2_pReal
+  real(pReal), parameter :: HGCoeff = 0.0e-2_pReal
-  real(pReal), dimension(3,3) :: &
+  real(pReal), parameter, dimension(4,8) :: &
-    temp33_Real = 0.0_pReal
-  real(pReal), dimension(4,8) :: &
     HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, &
                        1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, &
                       -1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, &
@@ -121,18 +119,19 @@ subroutine grid_mech_FEM_init
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
-  debug_grid => config_debug%get('grid', defaultVal=emptyList)
+  debug_grid => config_debug%get('grid',defaultVal=emptyList)
   debugRotation = debug_grid%contains('rotation')
 
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
+
   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=1.0e-3_pReal)
-  num%itmin           = num_grid%get_asInt  ('itmin',          defaultVal=1)
+  num%itmin           = num_grid%get_asInt  ('itmin',defaultVal=1)
-  num%itmax           = num_grid%get_asInt  ('itmax',          defaultVal=250)
+  num%itmax           = num_grid%get_asInt  ('itmax',defaultVal=250)
 
   if (num%eps_div_atol <= 0.0_pReal)             call IO_error(301,ext_msg='eps_div_atol')
   if (num%eps_div_rtol < 0.0_pReal)              call IO_error(301,ext_msg='eps_div_rtol')
@@ -225,6 +224,7 @@ subroutine grid_mech_FEM_init
     fileHandle  = HDF5_openFile(fileName)
     groupHandle = HDF5_openGroup(fileHandle,'solver')
 
+    call HDF5_read(groupHandle,P_aim,        'P_aim')
     call HDF5_read(groupHandle,F_aim,        'F_aim')
     call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
     call HDF5_read(groupHandle,F_aimDot,     'F_aimDot')
@@ -240,7 +240,7 @@ subroutine grid_mech_FEM_init
 
   homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])                          ! set starting condition for materialpoint_stressAndItsTangent
   call utilities_updateCoords(F)
-  call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, &                     ! stress field, stress avg, global average of stiffness and (min+max)/2
+  call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, &                            ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       F, &                                                          ! target F
                                       0.0_pReal)                                                    ! time increment
   call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
@@ -295,6 +295,7 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
   solution%iterationsNeeded = totalIter
   solution%termIll = terminallyIll
   terminallyIll = .false.
+  P_aim = merge(P_aim,P_av,params%stress_mask)
 
 end function grid_mech_FEM_solution
 
@@ -302,18 +303,17 @@ end function grid_mech_FEM_solution
 !--------------------------------------------------------------------------------------------------
 !> @brief forwarding routine
 !> @details find new boundary conditions and best F estimate for end of current timestep
-!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
+subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
                                  deformation_BC,stress_BC,rotation_BC)
 
   logical,                  intent(in) :: &
     cutBack, &
     guess
   real(pReal),              intent(in) :: &
-    timeinc_old, &
+    Delta_t_old, &
-    timeinc, &
+    Delta_t, &
-    loadCaseTime                                                                                    !< remaining time of current load case
+    t_remaining                                                                                     !< remaining time of current load case
   type(tBoundaryCondition), intent(in) :: &
     stress_BC, &
     deformation_BC
@@ -323,13 +323,6 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
   PetscScalar, pointer, dimension(:,:,:,:) :: &
     u_current,u_lastInc
 
-!--------------------------------------------------------------------------------------------------
-! set module wide available data
-  params%stress_mask = stress_BC%mask
-  params%rotation_BC = rotation_BC
-  params%timeinc     = timeinc
-  params%timeincOld  = timeinc_old
-
 
   call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
   call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
@@ -339,7 +332,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
   else
     C_volAvgLastInc    = C_volAvg
 
-    F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
+    F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)  ! estimate deformation rate for prescribed stress components
     F_aim_lastInc = F_aim
 
     !-----------------------------------------------------------------------------------------------
@@ -347,18 +340,18 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
     if     (deformation_BC%myType=='L') then                                                        ! calculate F_aimDot from given L and current F
       F_aimDot = F_aimDot &
                + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
-    elseif(deformation_BC%myType=='dot_F') then                                                     ! F_aimDot is prescribed
+    elseif (deformation_BC%myType=='dot_F') then                                                    ! F_aimDot is prescribed
       F_aimDot = F_aimDot &
                + merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
     elseif (deformation_BC%myType=='F') then                                                        ! aim at end of load case is prescribed
       F_aimDot = F_aimDot &
-               + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
+               + merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask)
     endif
 
     if (guess) then
       call VecWAXPY(solution_rate,-1.0_pReal,solution_lastInc,solution_current,ierr)
       CHKERRQ(ierr)
-      call VecScale(solution_rate,1.0_pReal/timeinc_old,ierr); CHKERRQ(ierr)
+      call VecScale(solution_rate,1.0_pReal/Delta_t_old,ierr); CHKERRQ(ierr)
     else
       call VecSet(solution_rate,0.0_pReal,ierr); CHKERRQ(ierr)
     endif
@@ -371,23 +364,28 @@ 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
+  F_aim = F_aim_lastInc + F_aimDot * Delta_t
-  if     (stress_BC%myType=='P') then
+  if (stress_BC%myType=='P')     P_aim = P_aim &
-    P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
+                                       + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t
-  elseif (stress_BC%myType=='dot_P') then !UNTESTED
+  if (stress_BC%myType=='dot_P') P_aim = P_aim &
-    P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
+                                       + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t
-  endif
 
-  call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
+  call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr)
 
-  call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
+  call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
-  call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr)
+  call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
+
+!--------------------------------------------------------------------------------------------------
+! set module wide available data
+  params%stress_mask = stress_BC%mask
+  params%rotation_BC = rotation_BC
+  params%timeinc     = Delta_t
 
 end subroutine grid_mech_FEM_forward
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Age
+!> @brief Update coordinates
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_FEM_updateCoords
 
@@ -415,6 +413,7 @@ subroutine grid_mech_FEM_restartWrite
   fileHandle  = HDF5_openFile(fileName,'w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
 
+  call HDF5_write(groupHandle,P_aim,        'P_aim')
   call HDF5_write(groupHandle,F_aim,        'F_aim')
   call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
   call HDF5_write(groupHandle,F_aimDot,     'F_aimDot')
@@ -458,10 +457,10 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
   divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
   BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
 
-  if ((totalIter >= num%itmin .and. &
+  if (terminallyIll .or. &
+      (totalIter >= num%itmin .and. &
        all([ err_div/divTol, &
-                                  err_BC /BCTol       ] < 1.0_pReal)) &
+             err_BC /BCTol   ] < 1.0_pReal))) then
-              .or.    terminallyIll) then
     reason = 1
   elseif (totalIter >= num%itmax) then
     reason = -1
@@ -665,16 +664,16 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
           C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ
   call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
   CHKERRQ(ierr)
-  call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
+  call DMGetGlobalVector(da_local,coordinates,ierr); CHKERRQ(ierr)
-  call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr)
+  call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr); CHKERRQ(ierr)
   ele = 0
   do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
     ele = ele + 1
     x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele)
   enddo; enddo; enddo
-  call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr)                       ! initialize to undeformed coordinates (ToDo: use ip coordinates)
+  call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr); CHKERRQ(ierr)                      ! initialize to undeformed coordinates (ToDo: use ip coordinates)
-  call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr)                          ! get rigid body deformation modes
+  call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr); CHKERRQ(ierr)                         ! get rigid body deformation modes
-  call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
+  call DMRestoreGlobalVector(da_local,coordinates,ierr); CHKERRQ(ierr)
   call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
   call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
   call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)

src/grid/grid_mech_spectral_basic.f90

@@ -94,8 +94,6 @@ contains
 subroutine grid_mech_spectral_basic_init
 
   real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
-  real(pReal), dimension(3,3) :: &
-    temp33_Real = 0.0_pReal
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
     F                                                                                               ! pointer to solution data
@@ -118,7 +116,7 @@ subroutine grid_mech_spectral_basic_init
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
-  debug_grid => config_debug%get('grid', defaultVal=emptyList)
+  debug_grid => config_debug%get('grid',defaultVal=emptyList)
   debugRotation = debug_grid%contains('rotation')
 
 !-------------------------------------------------------------------------------------------------
@@ -126,10 +124,10 @@ subroutine grid_mech_spectral_basic_init
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
 
   num%update_gamma    = num_grid%get_asBool ('update_gamma',   defaultVal=.false.)
-  num%eps_div_atol    = num_grid%get_asFloat ('eps_div_atol',   defaultVal=1.0e-4_pReal)
+  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_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_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
-  num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=1.0e-3_pReal)
+  num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
   num%itmin           = num_grid%get_asInt  ('itmin',defaultVal=1)
   num%itmax           = num_grid%get_asInt  ('itmax',defaultVal=250)
 
@@ -149,8 +147,8 @@ subroutine grid_mech_spectral_basic_init
 
 !--------------------------------------------------------------------------------------------------
 ! allocate global fields
-  allocate (F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
+  allocate(F_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
-  allocate (Fdot     (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
+  allocate(Fdot     (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
 
 !--------------------------------------------------------------------------------------------------
 ! initialize solver specific parts of PETSc
@@ -189,6 +187,7 @@ subroutine grid_mech_spectral_basic_init
     fileHandle  = HDF5_openFile(fileName)
     groupHandle = HDF5_openGroup(fileHandle,'solver')
 
+    call HDF5_read(groupHandle,P_aim,        'P_aim')
     call HDF5_read(groupHandle,F_aim,        'F_aim')
     call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
     call HDF5_read(groupHandle,F_aimDot,     'F_aimDot')
@@ -202,7 +201,7 @@ subroutine grid_mech_spectral_basic_init
 
   homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])                          ! set starting condition for materialpoint_stressAndItsTangent
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
-  call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, &                         ! stress field, stress avg, global average of stiffness and (min+max)/2
+  call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
                                       0.0_pReal)                                                    ! time increment
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! deassociate pointer
@@ -262,6 +261,7 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
   solution%iterationsNeeded = totalIter
   solution%termIll = terminallyIll
   terminallyIll = .false.
+  P_aim = merge(P_aim,P_av,params%stress_mask)
 
 end function grid_mech_spectral_basic_solution
 
@@ -269,18 +269,17 @@ end function grid_mech_spectral_basic_solution
 !--------------------------------------------------------------------------------------------------
 !> @brief forwarding routine
 !> @details find new boundary conditions and best F estimate for end of current timestep
-!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
+subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
                                             deformation_BC,stress_BC,rotation_BC)
 
   logical,                  intent(in) :: &
     cutBack, &
     guess
   real(pReal),              intent(in) :: &
-    timeinc_old, &
+    Delta_t_old, &
-    timeinc, &
+    Delta_t, &
-    loadCaseTime                                                                                     !< remaining time of current load case
+    t_remaining                                                                                     !< remaining time of current load case
   type(tBoundaryCondition), intent(in) :: &
     stress_BC, &
     deformation_BC
@@ -289,12 +288,6 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: F
 
-!--------------------------------------------------------------------------------------------------
-! set module wide available data
-  params%stress_mask = stress_BC%mask
-  params%rotation_BC = rotation_BC
-  params%timeinc     = timeinc
-  params%timeincOld  = timeinc_old
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
@@ -305,7 +298,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
     C_volAvgLastInc    = C_volAvg
     C_minMaxAvgLastInc = C_minMaxAvg
 
-    F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
+    F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)  ! estimate deformation rate for prescribed stress components
     F_aim_lastInc = F_aim
 
     !-----------------------------------------------------------------------------------------------
@@ -313,40 +306,45 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
     if     (deformation_BC%myType=='L') then                                                        ! calculate F_aimDot from given L and current F
       F_aimDot = F_aimDot &
                + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
-    elseif(deformation_BC%myType=='dot_F') then                                                     ! F_aimDot is prescribed
+    elseif (deformation_BC%myType=='dot_F') then                                                    ! F_aimDot is prescribed
       F_aimDot = F_aimDot &
                + merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
     elseif (deformation_BC%myType=='F') then                                                        ! aim at end of load case is prescribed
       F_aimDot = F_aimDot &
-               + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
+               + merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask)
     endif
 
     Fdot = utilities_calculateRate(guess, &
-                                   F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
+                                   F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
                                    rotation_BC%rotate(F_aimDot,active=.true.))
     F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
 
-    homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
+    homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
   endif
 
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
-  F_aim = F_aim_lastInc + F_aimDot * timeinc
+  F_aim = F_aim_lastInc + F_aimDot * Delta_t
-  if     (stress_BC%myType=='P') then
+  if (stress_BC%myType=='P')     P_aim = P_aim &
-    P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
+                                       + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t
-  elseif (stress_BC%myType=='dot_P') then !UNTESTED
+  if (stress_BC%myType=='dot_P') P_aim = P_aim &
-    P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
+                                       + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t
-  endif
 
-  F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! estimate of F at end of time+timeinc that matches rotated F_aim on average
+  F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, &                                      ! estimate of F at end of time+Delta_t 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)
 
+!--------------------------------------------------------------------------------------------------
+! set module wide available data
+  params%stress_mask = stress_BC%mask
+  params%rotation_BC = rotation_BC
+  params%timeinc     = Delta_t
+
 end subroutine grid_mech_spectral_basic_forward
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Age
+!> @brief Update coordinates
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_spectral_basic_updateCoords
 
@@ -378,6 +376,7 @@ subroutine grid_mech_spectral_basic_restartWrite
   fileHandle  = HDF5_openFile(fileName,'w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
 
+  call HDF5_write(groupHandle,P_aim,        'P_aim')
   call HDF5_write(groupHandle,F_aim,        'F_aim')
   call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
   call HDF5_write(groupHandle,F_aimDot,     'F_aimDot')

src/grid/grid_mech_spectral_polarisation.f90

@@ -105,8 +105,6 @@ contains
 subroutine grid_mech_spectral_polarisation_init
 
   real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
-  real(pReal), dimension(3,3) :: &
-    temp33_Real = 0.0_pReal
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
     FandF_tau, &                                                                                    ! overall pointer to solution data
@@ -211,6 +209,7 @@ subroutine grid_mech_spectral_polarisation_init
     fileHandle  = HDF5_openFile(fileName)
     groupHandle = HDF5_openGroup(fileHandle,'solver')
 
+    call HDF5_read(groupHandle,P_aim,        'P_aim')
     call HDF5_read(groupHandle,F_aim,        'F_aim')
     call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
     call HDF5_read(groupHandle,F_aimDot,     'F_aimDot')
@@ -228,7 +227,7 @@ subroutine grid_mech_spectral_polarisation_init
 
   homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3])                          ! set starting condition for materialpoint_stressAndItsTangent
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
-  call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, &                         ! stress field, stress avg, global average of stiffness and (min+max)/2
+  call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
                                       0.0_pReal)                                                    ! time increment
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                        ! deassociate pointer
@@ -294,6 +293,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
   solution%iterationsNeeded = totalIter
   solution%termIll = terminallyIll
   terminallyIll = .false.
+  P_aim = merge(P_aim,P_av,params%stress_mask)
 
 end function grid_mech_spectral_polarisation_solution
 
@@ -301,18 +301,17 @@ end function grid_mech_spectral_polarisation_solution
 !--------------------------------------------------------------------------------------------------
 !> @brief forwarding routine
 !> @details find new boundary conditions and best F estimate for end of current timestep
-!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
+subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
                                                    deformation_BC,stress_BC,rotation_BC)
 
   logical,                  intent(in) :: &
     cutBack, &
     guess
   real(pReal),              intent(in) :: &
-    timeinc_old, &
+    Delta_t_old, &
-    timeinc, &
+    Delta_t, &
-    loadCaseTime                                                                                    !< remaining time of current load case
+    t_remaining                                                                                     !< remaining time of current load case
   type(tBoundaryCondition), intent(in) :: &
     stress_BC, &
     deformation_BC
@@ -323,12 +322,6 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
   integer :: i, j, k
   real(pReal), dimension(3,3) :: F_lambda33
 
-!--------------------------------------------------------------------------------------------------
-! set module wide available data
-  params%stress_mask = stress_BC%mask
-  params%rotation_BC = rotation_BC
-  params%timeinc     = timeinc
-  params%timeincOld  = timeinc_old
 
   call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
   F     => FandF_tau(0: 8,:,:,:)
@@ -341,7 +334,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
     C_volAvgLastInc    = C_volAvg
     C_minMaxAvgLastInc = C_minMaxAvg
 
-    F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
+    F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)  ! estimate deformation rate for prescribed stress components
     F_aim_lastInc = F_aim
 
     !-----------------------------------------------------------------------------------------------
@@ -349,19 +342,19 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
     if     (deformation_BC%myType=='L') then                                                        ! calculate F_aimDot from given L and current F
       F_aimDot = F_aimDot &
                + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
-    elseif(deformation_BC%myType=='dot_F') then                                                     ! F_aimDot is prescribed
+    elseif (deformation_BC%myType=='dot_F') then                                                    ! F_aimDot is prescribed
       F_aimDot = F_aimDot &
                + merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
     elseif (deformation_BC%myType=='F') then                                                        ! aim at end of load case is prescribed
       F_aimDot = F_aimDot &
-               + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
+               + merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask)
     endif
 
     Fdot     = utilities_calculateRate(guess, &
-                                       F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
+                                       F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),Delta_t_old, &
                                        rotation_BC%rotate(F_aimDot,active=.true.))
     F_tauDot = utilities_calculateRate(guess, &
-                                       F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
+                                       F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), Delta_t_old, &
                                        rotation_BC%rotate(F_aimDot,active=.true.))
     F_lastInc     = reshape(F,    [3,3,grid(1),grid(2),grid3])
     F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
@@ -371,38 +364,41 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 
 !--------------------------------------------------------------------------------------------------
 ! update average and local deformation gradients
-  F_aim = F_aim_lastInc + F_aimDot * timeinc
+  F_aim = F_aim_lastInc + F_aimDot * Delta_t
-  if     (stress_BC%myType=='P') then
+  if(stress_BC%myType=='P')     P_aim = P_aim &
-    P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
+                                      + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t
-  elseif (stress_BC%myType=='dot_P') then !UNTESTED
+  if(stress_BC%myType=='dot_P') P_aim = P_aim &
-    P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
+                                      + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t
-  endif
 
-  F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! estimate of F at end of time+timeinc that matches rotated F_aim on average
+  F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, &                                      ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
                                      rotation_BC%rotate(F_aim,active=.true.)),&
               [9,grid(1),grid(2),grid3])
   if (guess) then
-     F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
+     F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), &
                      [9,grid(1),grid(2),grid3])                                                     ! does not have any average value as boundary condition
    else
     do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
        F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
-       F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
+       F_lambda33 = math_I3 &
-                                   math_mul3333xx33(C_scale,&
+                  + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, &
-                                                    matmul(transpose(F_lambda33),&
+                    math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-math_I3)))
-                                                           F_lambda33)-math_I3))*0.5_pReal) &
-                  + math_I3
        F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
     enddo; enddo; enddo
   endif
 
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
 
+!--------------------------------------------------------------------------------------------------
+! set module wide available data
+  params%stress_mask = stress_BC%mask
+  params%rotation_BC = rotation_BC
+  params%timeinc     = Delta_t
+
 end subroutine grid_mech_spectral_polarisation_forward
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Age
+!> @brief Update coordinates
 !--------------------------------------------------------------------------------------------------
 subroutine grid_mech_spectral_polarisation_updateCoords
 
@@ -436,6 +432,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite
   fileHandle  = HDF5_openFile(fileName,'w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
 
+  call HDF5_write(groupHandle,F_aim,        'P_aim')
   call HDF5_write(groupHandle,F_aim,        'F_aim')
   call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
   call HDF5_write(groupHandle,F_aimDot,     'F_aimDot')
@@ -480,11 +477,11 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
   divTol     = max(maxval(abs(P_av))         *num%eps_div_rtol   ,num%eps_div_atol)
   BCTol      = max(maxval(abs(P_av))         *num%eps_stress_rtol,num%eps_stress_atol)
 
-  if ((totalIter >= num%itmin .and. &
+  if (terminallyIll .or. &
+      (totalIter >= num%itmin .and. &
        all([ err_div /divTol, &
              err_curl/curlTol, &
-                                  err_BC  /BCTol       ] < 1.0_pReal)) &
+             err_BC  /BCTol   ] < 1.0_pReal))) then
-              .or.    terminallyIll) then
     reason = 1
   elseif (totalIter >= num%itmax) then
     reason = -1

src/grid/spectral_utilities.f90

@@ -93,7 +93,7 @@ module spectral_utilities
     real(pReal), dimension(3,3) :: stress_BC
     logical, dimension(3,3)     :: stress_mask
     type(rotation)              :: rotation_BC
-    real(pReal) :: timeinc, timeincOld
+    real(pReal) :: timeinc
   end type tSolutionParams
 
   type :: tNumerics

src/mesh/mesh_mech_FEM.f90

@@ -32,7 +32,6 @@ module mesh_mech_FEM
   type tSolutionParams
     type(tFieldBC)  :: fieldBC
     real(pReal)     :: timeinc
-    real(pReal)     :: timeincOld
   end type tSolutionParams
 
   type(tSolutionParams)  :: params
@@ -283,7 +282,6 @@ type(tSolutionState) function FEM_mech_solution( &
 !--------------------------------------------------------------------------------------------------
 ! set module wide availabe data
   params%timeinc = timeinc
-  params%timeincOld = timeinc_old
   params%fieldBC = fieldBC
 
   call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)                             ! solve mech_snes based on solution guess (result in solution)