From 32bb0d8c6e5b4b223a8591133b0d2283b755dcd6 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 20 Jan 2021 20:54:31 +0100 Subject: [PATCH] new variables for damage --- PRIVATE | 2 +- src/commercialFEM_fileList.f90 | 1 + src/constitutive.f90 | 12 ++++++ src/constitutive_damage.f90 | 42 ++++++++++++++++++- src/grid/grid_damage_spectral.f90 | 69 ++++++++++++++++--------------- src/homogenization.f90 | 11 ++++- src/homogenization_damage.f90 | 40 ++++++++++++++++++ 7 files changed, 141 insertions(+), 36 deletions(-) create mode 100644 src/homogenization_damage.f90 diff --git a/PRIVATE b/PRIVATE index 9e8011876..5fd23c60c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 9e8011876b16896a18902737790eba8018e94f85 +Subproject commit 5fd23c60c721aada27bcd3a4ba96301753b7bd0d diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index e04936a3b..963871c7a 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -52,4 +52,5 @@ #include "homogenization_mech_isostrain.f90" #include "homogenization_mech_RGC.f90" #include "homogenization_thermal.f90" +#include "homogenization_damage.f90" #include "CPFEM.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7ecff0eef..16f3d3b59 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -188,6 +188,11 @@ module constitutive real(pReal), dimension(3,3) :: P end function constitutive_mech_getP + module function constitutive_damage_get_phi(co,ip,el) result(phi) + integer, intent(in) :: co, ip, el + real(pReal) :: phi + end function constitutive_damage_get_phi + module function thermal_T(ph,me) result(T) integer, intent(in) :: ph,me real(pReal) :: T @@ -203,6 +208,11 @@ module constitutive real(pReal), intent(in) :: T integer, intent(in) :: co, ce end subroutine constitutive_thermal_setT + + module subroutine constitutive_damage_set_phi(phi,co,ce) + real(pReal), intent(in) :: phi + integer, intent(in) :: co, ce + end subroutine constitutive_damage_set_phi ! == cleaned:end =================================================================================== @@ -353,6 +363,8 @@ module constitutive constitutive_restartRead, & integrateDamageState, & constitutive_thermal_setT, & + constitutive_damage_set_phi, & + constitutive_damage_get_phi, & constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index a8792d5e2..07bfbeeca 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -10,9 +10,16 @@ submodule(constitutive) constitutive_damage DAMAGE_ANISODUCTILE_ID end enum + + type :: tDataContainer + real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi + end type tDataContainer + integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & phase_source !< active sources mechanisms of each phase + type(tDataContainer), dimension(:), allocatable :: current + interface module function source_damage_anisoBrittle_init(source_length) result(mySources) @@ -152,7 +159,8 @@ contains module subroutine damage_init integer :: & - ph !< counter in phase loop + ph, & !< counter in phase loop + Nconstituents class(tNode), pointer :: & phases, & phase, & @@ -161,10 +169,18 @@ module subroutine damage_init phases => config_material%get('phase') + allocate(current(phases%length)) + allocate(damageState (phases%length)) allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics do ph = 1,phases%length + + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + + allocate(current(ph)%phi(Nconstituents),source=1.0_pReal) + allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal) + phase => phases%get(ph) sources => phase%get('source',defaultVal=emptyList) phase_Nsources(ph) = sources%length @@ -511,4 +527,28 @@ function source_active(source_label,src_length) result(active_source) end function source_active +!---------------------------------------------------------------------------------------------- +!< @brief Set damage parameter +!---------------------------------------------------------------------------------------------- +module subroutine constitutive_damage_set_phi(phi,co,ce) + + real(pReal), intent(in) :: phi + integer, intent(in) :: ce, co + + + current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi + +end subroutine constitutive_damage_set_phi + + +module function constitutive_damage_get_phi(co,ip,el) result(phi) + + integer, intent(in) :: co, ip, el + real(pReal) :: phi + + phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el)) + +end function constitutive_damage_get_phi + + end submodule constitutive_damage diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 79437945b..a891f070d 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -17,8 +17,9 @@ module grid_damage_spectral use discretization_grid use damage_nonlocal use YAML_types + use homogenization use config - + implicit none private @@ -43,13 +44,13 @@ module grid_damage_spectral phi_current, & !< field of current damage phi_lastInc, & !< field of previous damage phi_stagInc !< field of staggered damage - + !-------------------------------------------------------------------------------------------------- -! reference diffusion tensor, mobility etc. +! reference diffusion tensor, mobility etc. integer :: totalIter = 0 !< total iteration in current increment real(pReal), dimension(3,3) :: K_ref real(pReal) :: mu_ref - + public :: & grid_damage_spectral_init, & grid_damage_spectral_solution, & @@ -62,8 +63,8 @@ contains ! ToDo: Restart not implemented !-------------------------------------------------------------------------------------------------- subroutine grid_damage_spectral_init - - PetscInt, dimension(0:worldsize-1) :: localK + + PetscInt, dimension(0:worldsize-1) :: localK DM :: damage_grid Vec :: uBound, lBound PetscErrorCode :: ierr @@ -72,22 +73,22 @@ subroutine grid_damage_spectral_init num_generic character(len=pStringLen) :: & snes_type - + print'(/,a)', ' <<<+- grid_spectral_damage init -+>>>' print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal) num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal) - + num_generic => config_numerics%get('generic',defaultVal=emptyDict) num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) - + if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%eps_damage_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_damage_atol') @@ -104,7 +105,7 @@ subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) @@ -122,7 +123,7 @@ subroutine grid_damage_spectral_init call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector - CHKERRQ(ierr) + CHKERRQ(ierr) call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) if (trim(snes_type) == 'vinewtonrsls' .or. & @@ -137,12 +138,12 @@ subroutine grid_damage_spectral_init endif !-------------------------------------------------------------------------------------------------- -! init fields +! init fields call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) CHKERRQ(ierr) xend = xstart + xend - 1 yend = ystart + yend - 1 - zend = zstart + zend - 1 + zend = zstart + zend - 1 allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) @@ -157,26 +158,26 @@ end subroutine grid_damage_spectral_init !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- function grid_damage_spectral_solution(timeinc) result(solution) - + real(pReal), intent(in) :: & timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull PetscReal :: phi_min, phi_max, stagNorm, solnNorm - + PetscErrorCode :: ierr SNESConvergedReason :: reason solution%converged =.false. - + !-------------------------------------------------------------------------------------------------- -! set module wide availabe data +! set module wide availabe data params%timeinc = timeinc - + call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) - + if (reason < 1) then solution%converged = .false. solution%iterationsNeeded = num%itmax @@ -192,20 +193,21 @@ function grid_damage_spectral_solution(timeinc) result(solution) solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm) !-------------------------------------------------------------------------------------------------- -! updating damage state +! updating damage state cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo - + call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr) if (solution%converged) & print'(/,a)', ' ... nonlocal damage converged .....................................' print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage = ', phi_min, phi_max, stagNorm print'(/,a)', ' ===========================================================================' - flush(IO_STDOUT) + flush(IO_STDOUT) end function grid_damage_spectral_solution @@ -214,31 +216,32 @@ end function grid_damage_spectral_solution !> @brief spectral damage forwarding routine !-------------------------------------------------------------------------------------------------- subroutine grid_damage_spectral_forward(cutBack) - + logical, intent(in) :: cutBack integer :: i, j, k, cell DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: x_scal PetscErrorCode :: ierr - if (cutBack) then + if (cutBack) then phi_current = phi_lastInc phi_stagInc = phi_lastInc !-------------------------------------------------------------------------------------------------- -! reverting damage field state +! reverting damage field state cell = 0 call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) - cell = cell + 1 + cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) + homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo else phi_lastInc = phi_current call updateReference - endif + endif end subroutine grid_damage_spectral_forward @@ -247,7 +250,7 @@ end subroutine grid_damage_spectral_forward !> @brief forms the spectral damage residual vector !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,f_scal,dummy,ierr) - + DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: & in PetscScalar, dimension( & @@ -261,11 +264,11 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) integer :: i, j, k, cell real(pReal) :: phiDot, dPhiDot_dPhi, mobility - phi_current = x_scal + phi_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real = 0.0_pReal - scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current + scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field call utilities_FFTvectorBackward @@ -297,7 +300,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) & scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness - + !-------------------------------------------------------------------------------------------------- ! constructing residual f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current @@ -311,7 +314,7 @@ end subroutine formResidual subroutine updateReference integer :: i,j,k,cell,ierr - + cell = 0 K_ref = 0.0_pReal mu_ref = 0.0_pReal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index d576c5672..ed9b7fead 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -64,6 +64,9 @@ module homogenization module subroutine thermal_init end subroutine thermal_init + module subroutine damage_init + end subroutine damage_init + module subroutine mech_partition(subF,ip,el) real(pReal), intent(in), dimension(3,3) :: & subF @@ -77,6 +80,11 @@ module homogenization integer, intent(in) :: ce end subroutine thermal_partition + module subroutine damage_partition(phi,ce) + real(pReal), intent(in) :: phi + integer, intent(in) :: ce + end subroutine damage_partition + module subroutine thermal_homogenize(ip,el) integer, intent(in) :: ip,el end subroutine thermal_homogenize @@ -120,7 +128,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !-------------------------------------------------------------------------------------------------- -subroutine homogenization_init +subroutine homogenization_init() class (tNode) , pointer :: & num_homog, & @@ -144,6 +152,7 @@ subroutine homogenization_init call mech_init(num_homog) call thermal_init() + call damage_init() if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T) if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 new file mode 100644 index 000000000..479318340 --- /dev/null +++ b/src/homogenization_damage.f90 @@ -0,0 +1,40 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!-------------------------------------------------------------------------------------------------- +submodule(homogenization) homogenization_damage + + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief Allocate variables and set parameters. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_init() + + + print'(/,a)', ' <<<+- homogenization_damage init -+>>>' + + allocate(homogenization_phi(discretization_nIPs*discretization_Nelems)) + allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems)) + +end subroutine damage_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Partition temperature onto the individual constituents. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_partition(phi,ce) + + real(pReal), intent(in) :: phi + integer, intent(in) :: ce + + integer :: co + + do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + call constitutive_damage_set_phi(phi,co,ce) + enddo + +end subroutine damage_partition + + +end submodule homogenization_damage