From 6367cb8fcba017330d762273e4ee38fc514fe187 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 13:29:51 +0200 Subject: [PATCH] consistent handling of boundary conditions A stress boundary condition 'P' indicates the stress at the end of the load case (same as for 'F') 'Pdot' for given increase of (technical) strain is not implemented. Does not change anything for the most common case of zero-stress boundary conditions, but simplifies the specification of stress ramps --- src/grid/DAMASK_grid.f90 | 6 +++--- src/grid/grid_mech_FEM.f90 | 19 +++++++++--------- src/grid/grid_mech_spectral_basic.f90 | 21 ++++++++++++-------- src/grid/grid_mech_spectral_polarisation.f90 | 15 +++++++++----- 4 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 59a5b452c..54d080744 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -107,8 +107,8 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! initialize field solver information nActiveFields = 1 - if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 - if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 + if(any(thermal_type == THERMAL_conduction_ID)) nActiveFields = nActiveFields + 1 + if(any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) allocate(newLoadCase%ID(nActiveFields)) @@ -181,7 +181,7 @@ program DAMASK_grid if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase - newLoadCase%stress%myType='stress' + newLoadCase%stress%myType='p' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index bd6f862b2..9efd88505 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -66,8 +66,8 @@ module grid_mech_FEM F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -327,7 +327,6 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -374,6 +373,12 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + endif + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) @@ -489,8 +494,6 @@ subroutine formResidual(da_local,x_local, & PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal PetscScalar, dimension(8,3) :: x_elem, f_elem PetscInt :: i, ii, j, jj, k, kk, ctr, ele - real(pReal), dimension(3,3) :: & - deltaF_aim PetscInt :: & PETScIter, & nfuncs @@ -539,10 +542,8 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index b199c723f..cd795f005 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -62,8 +62,8 @@ module grid_mech_spectral_basic F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -129,7 +129,7 @@ subroutine grid_mech_spectral_basic_init num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.001_pReal) num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) @@ -292,7 +292,6 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -333,6 +332,12 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + endif + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -484,9 +489,9 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) ! S = 0.0 for no bc - F_aim = F_aim -deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + F_aim = F_aim - deltaF_aim + err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme @@ -499,7 +504,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! constructing residual - residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 393f8f3f7..4dc9d6d7c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -70,8 +70,8 @@ module grid_mech_spectral_polarisation F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -326,7 +326,6 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -373,6 +372,12 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + endif + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) @@ -587,10 +592,10 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, P_av - params%stress_BC) ! S = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & -params%rotation_BC%rotate(F_av)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + params%stress_mask * (P_av-P_aim))) ! mask = 0.0 for no bc ! calculate divergence tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise