From 2dd520b4a289a437f37f396fcf1f529688743948 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Dec 2020 08:06:50 +0100 Subject: [PATCH] P_aim should be independent from P_av P_av is not defined after restart or cutback. Restart with change of load case is probably still an issue --- PRIVATE | 2 +- python/damask/grid_filters.py | 4 +- src/grid/DAMASK_grid.f90 | 111 +++++++++-------- src/grid/grid_mech_FEM.f90 | 123 +++++++++---------- src/grid/grid_mech_spectral_basic.f90 | 97 ++++++++------- src/grid/grid_mech_spectral_polarisation.f90 | 107 ++++++++-------- src/grid/spectral_utilities.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 2 - 8 files changed, 224 insertions(+), 224 deletions(-) diff --git a/PRIVATE b/PRIVATE index 08f8aea46..de65e1df5 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 08f8aea465a1b5e476b584bcae7927d113919b1d +Subproject commit de65e1df5a76362de93667e9820dbf330b56f96d diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 9e06e075a..ee929b3bf 100644 --- a/python/damask/grid_filters.py +++ b/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 \ diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index a8271cffc..514443dbb 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/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) - logical, dimension(9) :: temp_maskVector = .false. !< temporarily from loadcase file when reading in tensors + real(pReal), dimension(9) :: temp_valueVector !< temporarily from loadcase file when reading in tensors (initialize to 0.0) + 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 @@ -178,11 +174,11 @@ program DAMASK_grid allocate(loadCases(l)%ID(nActiveFields)) field = 1 loadCases(l)%ID(field) = FIELD_MECH_ID ! mechanical active by default - thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then + thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then field = field + 1 loadCases(l)%ID(field) = FIELD_THERMAL_ID endif thermalActive - damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then + damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then field = field + 1 loadCases(l)%ID(field) = FIELD_DAMAGE_ID endif damageActive @@ -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' - if (temp_maskVector(j)) temp_valueVector(j) = step_stress%get_asFloat(j) ! read value where applicable + temp_maskVector(j) = step_deformation%get_asString(j) /= 'x' + 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:' - else if (loadCases(l)%deformation%myType == 'F') then + endif + if (loadCases(l)%deformation%myType == 'F') then print*, ' F:' - else if (loadCases(l)%deformation%myType == 'dot_F') then - print*, ' dot_F:' + else + 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 .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) & - errorID = 838 ! no rotation is allowed by stress BC - print*, ' P / MPa:' - do i = 1, 3; do j = 1, 3 - 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 + 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 + + 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 + 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)%r <= 0.0) errorID = 833 + 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 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 4394b6f81..146f28567 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -31,16 +31,16 @@ module grid_mech_FEM type :: tNumerics integer :: & - itmin, & !< minimum number of iterations - itmax !< maximum number of iterations + itmin, & !< minimum number of iterations + itmax !< maximum number of iterations real(pReal) :: & - eps_div_atol, & !< absolute tolerance for equilibrium - eps_div_rtol, & !< relative tolerance for equilibrium - eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC - eps_stress_rtol !< relative tolerance for fullfillment of stress BC + eps_div_atol, & !< absolute tolerance for equilibrium + eps_div_rtol, & !< relative tolerance for equilibrium + eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC + eps_stress_rtol !< relative tolerance for fullfillment of stress BC end type tNumerics - type(tNumerics) :: num ! numerics parameters. Better name? + type(tNumerics) :: num ! numerics parameters. Better name? logical :: debugRotation @@ -64,7 +64,7 @@ module grid_mech_FEM real(pReal), dimension(3,3) :: & 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 + F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information @@ -93,10 +93,8 @@ contains !-------------------------------------------------------------------------------------------------- subroutine grid_mech_FEM_init - real(pReal) :: HGCoeff = 0.0e-2_pReal - real(pReal), dimension(3,3) :: & - temp33_Real = 0.0_pReal - real(pReal), dimension(4,8) :: & + real(pReal), parameter :: HGCoeff = 0.0e-2_pReal + real(pReal), parameter, 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%itmax = num_grid%get_asInt ('itmax', defaultVal=250) + num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) + 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') @@ -238,9 +238,9 @@ subroutine grid_mech_FEM_init F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + 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,34 +303,26 @@ 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) :: & + logical, intent(in) :: & cutBack, & guess - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - type(tBoundaryCondition), intent(in) :: & + real(pReal), intent(in) :: & + Delta_t_old, & + Delta_t, & + t_remaining !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr 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 - F_aimDot = F_aimDot & + 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 - if (stress_BC%myType=='P') then - P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) - elseif (stress_BC%myType=='dot_P') then !UNTESTED - P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) - endif + F_aim = F_aim_lastInc + F_aimDot * Delta_t + if (stress_BC%myType=='P') P_aim = P_aim & + + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + if (stress_BC%myType=='dot_P') P_aim = P_aim & + + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t - 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_lastInc,u_lastInc,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) + +!-------------------------------------------------------------------------------------------------- +! 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') @@ -441,11 +440,11 @@ end subroutine grid_mech_FEM_restartWrite subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr) SNES :: snes_local - PetscInt, intent(in) :: PETScIter - PetscReal, intent(in) :: & - devNull1, & - devNull2, & - fnorm + PetscInt, intent(in) :: PETScIter + PetscReal, intent(in) :: & + devNull1, & + devNull2, & + fnorm SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: ierr @@ -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. & - all([ err_div/divTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then + if (terminallyIll .or. & + (totalIter >= num%itmin .and. & + all([ err_div/divTol, & + err_BC /BCTol ] < 1.0_pReal))) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -666,16 +665,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 DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) + call DMGetGlobalVector(da_local,coordinates,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 MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes - call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr) + 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 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) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 563b25162..4f5ceff61 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/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,20 +116,20 @@ 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') - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks 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_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%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + 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_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%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') @@ -139,7 +137,7 @@ subroutine grid_mech_spectral_basic_init if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') - + !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) @@ -149,14 +147,14 @@ subroutine grid_mech_spectral_basic_init !-------------------------------------------------------------------------------------------------- ! allocate global fields - 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(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) !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 + localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3d(PETSC_COMM_WORLD, & @@ -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') @@ -200,9 +199,9 @@ subroutine grid_mech_spectral_basic_init F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + 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,32 +269,25 @@ 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) :: & + logical, intent(in) :: & cutBack, & guess - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - type(tBoundaryCondition), intent(in) :: & + real(pReal), intent(in) :: & + Delta_t_old, & + Delta_t, & + t_remaining !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC 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 - F_aimDot = F_aimDot & + 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 - if (stress_BC%myType=='P') then - P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) - elseif (stress_BC%myType=='dot_P') then !UNTESTED - P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) - 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_aim = F_aim_lastInc + F_aimDot * Delta_t + if (stress_BC%myType=='P') P_aim = P_aim & + + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + if (stress_BC%myType=='dot_P') P_aim = P_aim & + + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + + 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') @@ -463,7 +462,7 @@ subroutine formResidual(in, F, & PetscErrorCode :: ierr call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) + call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment @@ -472,7 +471,7 @@ subroutine formResidual(in, F, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax - if (debugRotation) & + if(debugRotation) & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 03780f2e0..d09b7fcb2 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/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 @@ -147,16 +145,16 @@ subroutine grid_mech_spectral_polarisation_init num%alpha = num_grid%get_asFloat('alpha', defaultVal=1.0_pReal) num%beta = num_grid%get_asFloat('beta', defaultVal=1.0_pReal) - 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') - if (num%eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol') - if (num%eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol') - if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') - if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') - if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') - if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') - if (num%alpha <= 0.0_pReal .or. num%alpha > 2.0_pReal) call IO_error(301,ext_msg='alpha') - if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal) call IO_error(301,ext_msg='beta') + 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') + if (num%eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol') + if (num%eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol') + if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol') + if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') + if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') + if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') + if (num%alpha <= 0.0_pReal .or. num%alpha > 2.0_pReal) call IO_error(301,ext_msg='alpha') + if (num%beta < 0.0_pReal .or. num%beta > 2.0_pReal) call IO_error(301,ext_msg='beta') !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc @@ -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') @@ -226,9 +225,9 @@ subroutine grid_mech_spectral_polarisation_init F_tau_lastInc = 2.0_pReal*F_lastInc endif restartRead - homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent + 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,34 +301,27 @@ 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) :: & + logical, intent(in) :: & cutBack, & guess - real(pReal), intent(in) :: & - timeinc_old, & - timeinc, & - loadCaseTime !< remaining time of current load case - type(tBoundaryCondition), intent(in) :: & + real(pReal), intent(in) :: & + Delta_t_old, & + Delta_t, & + t_remaining !< remaining time of current load case + type(tBoundaryCondition), intent(in) :: & stress_BC, & deformation_BC - type(rotation), intent(in) :: & + type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau 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 - F_aimDot = F_aimDot & + 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 - if (stress_BC%myType=='P') then - P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) - elseif (stress_BC%myType=='dot_P') then !UNTESTED - P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) - endif + F_aim = F_aim_lastInc + F_aimDot * Delta_t + if(stress_BC%myType=='P') P_aim = P_aim & + + merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + if(stress_BC%myType=='dot_P') P_aim = P_aim & + + merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t - 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, & - math_mul3333xx33(C_scale,& - matmul(transpose(F_lambda33),& - F_lambda33)-math_I3))*0.5_pReal) & - + math_I3 + F_lambda33 = math_I3 & + + math_mul3333xx33(S_scale,0.5_pReal*matmul(F_lambda33, & + math_mul3333xx33(C_scale,matmul(transpose(F_lambda33),F_lambda33)-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. & - all([ err_div /divTol, & - err_curl/curlTol, & - err_BC /BCTol ] < 1.0_pReal)) & - .or. terminallyIll) then + if (terminallyIll .or. & + (totalIter >= num%itmin .and. & + all([ err_div /divTol, & + err_curl/curlTol, & + err_BC /BCTol ] < 1.0_pReal))) then reason = 1 elseif (totalIter >= num%itmax) then reason = -1 @@ -555,7 +552,7 @@ subroutine formResidual(in, FandF_tau, & newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax - if(debugRotation) & + if (debugRotation) & write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 1e1608d7c..6d6e26cae 100644 --- a/src/grid/spectral_utilities.f90 +++ b/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 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 8aa084ac8..b6ce1e175 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/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 @@ -302,7 +301,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)