223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'polishing-for-release' into 'development' 

Polishing for release

See merge request damask/DAMASK!334
This commit is contained in:
Franz Roters 2021-02-12 18:46:04 +00:00
parent 992b4a7e6d 958936c449
commit b7b764ab5f
45 changed files with 791 additions and 775 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 5d27d879abaeb1542e3f9f3065172be740af2899 Subproject commit e74cf00628285a587ced1e551cc9837c1011ca1c

3
examples/SpectralMethod/Polycrystal/shearXY.yaml
View File

@ -3,9 +3,6 @@
solver: solver:
mechanical: spectral_basic mechanical: spectral_basic
initial_conditions:
T: 300 #in Kelvin
loadstep: loadstep:
- boundary_conditions: - boundary_conditions:
mechanical: mechanical:

3
examples/SpectralMethod/Polycrystal/shearZX.yaml
View File

@ -3,9 +3,6 @@
solver: solver:
mechanical: spectral_basic mechanical: spectral_basic
initial_conditions:
T: 300 #in Kelvin
loadstep: loadstep:
- boundary_conditions: - boundary_conditions:
mechanical: mechanical:

3
examples/SpectralMethod/Polycrystal/tensionX.yaml
View File

@ -3,9 +3,6 @@
solver: solver:
mechanical: spectral_basic mechanical: spectral_basic
initial_conditions:
T: 300 #in Kelvin
loadstep: loadstep:
- boundary_conditions: - boundary_conditions:
mechanical: mechanical:

14
python/damask/_configmaterial.py
View File

@ -168,11 +168,11 @@ class ConfigMaterial(Config):
ok = False ok = False
if 'material' in self: if 'material' in self:
for i,v in enumerate(self['material']): for i,m in enumerate(self['material']):
if 'constituents' in v: if 'constituents' in m:
f = 0.0 v = 0.0
for c in v['constituents']: for c in m['constituents']:
f+= float(c['v']) v+= float(c['v'])
if 'O' in c: if 'O' in c:
try: try:
Rotation.from_quaternion(c['O']) Rotation.from_quaternion(c['O'])
@ -180,8 +180,8 @@ class ConfigMaterial(Config):
o = c['O'] o = c['O']
print(f"Invalid orientation: '{o}' in material '{i}'") print(f"Invalid orientation: '{o}' in material '{i}'")
ok = False ok = False
if not np.isclose(f,1.0): if not np.isclose(v,1.0):
print(f"Invalid total fraction '{f}' in material '{i}'") print(f"Invalid total fraction (v) '{v}' in material '{i}'")
ok = False ok = False
return ok return ok

30
python/damask/_vtk.py
View File

@ -347,23 +347,21 @@ class VTK:
See http://compilatrix.com/article/vtk-1 for further ideas. See http://compilatrix.com/article/vtk-1 for further ideas.
""" """
def screen_size(): try:
import wx
_ = wx.App(False) # noqa
width, height = wx.GetDisplaySize()
except ImportError:
try: try:
import wx import tkinter
_ = wx.App(False) # noqa tk = tkinter.Tk()
width, height = wx.GetDisplaySize() width = tk.winfo_screenwidth()
except ImportError: height = tk.winfo_screenheight()
try: tk.destroy()
import tkinter except Exception as e:
tk = tkinter.Tk() width = 1024
width = tk.winfo_screenwidth() height = 768
height = tk.winfo_screenheight()
tk.destroy()
except Exception as e:
width = 1024
height = 768
return (width,height)
mapper = vtk.vtkDataSetMapper() mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(self.vtk_data) mapper.SetInputData(self.vtk_data)
actor = vtk.vtkActor() actor = vtk.vtkActor()
@ -377,7 +375,7 @@ class VTK:
ren.AddActor(actor) ren.AddActor(actor)
ren.SetBackground(0.2,0.2,0.2) ren.SetBackground(0.2,0.2,0.2)
window.SetSize(screen_size[0],screen_size[1]) window.SetSize(width,height)
iren = vtk.vtkRenderWindowInteractor() iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(window) iren.SetRenderWindow(window)

6
src/CPFEM.f90
View File

@ -85,7 +85,7 @@ subroutine CPFEM_initAll
call discretization_marc_init call discretization_marc_init
call lattice_init call lattice_init
call material_init(.false.) call material_init(.false.)
call constitutive_init call phase_init
call homogenization_init call homogenization_init
call crystallite_init call crystallite_init
call CPFEM_init call CPFEM_init
@ -257,7 +257,7 @@ end subroutine CPFEM_general
subroutine CPFEM_forward subroutine CPFEM_forward
call homogenization_forward call homogenization_forward
call constitutive_forward call phase_forward
end subroutine CPFEM_forward end subroutine CPFEM_forward
@ -272,7 +272,7 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call phase_results
call homogenization_results call homogenization_results
call discretization_results call discretization_results
call results_finalizeIncrement call results_finalizeIncrement

10
src/CPFEM2.f90
View File

@ -60,7 +60,7 @@ subroutine CPFEM_initAll
call discretization_grid_init(restart=interface_restartInc>0) call discretization_grid_init(restart=interface_restartInc>0)
#endif #endif
call material_init(restart=interface_restartInc>0) call material_init(restart=interface_restartInc>0)
call constitutive_init call phase_init
call homogenization_init call homogenization_init
call crystallite_init call crystallite_init
call CPFEM_init call CPFEM_init
@ -87,7 +87,7 @@ subroutine CPFEM_init
fileHandle = HDF5_openFile(fileName) fileHandle = HDF5_openFile(fileName)
call homogenization_restartRead(fileHandle) call homogenization_restartRead(fileHandle)
call constitutive_restartRead(fileHandle) call phase_restartRead(fileHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
endif endif
@ -110,7 +110,7 @@ subroutine CPFEM_restartWrite
fileHandle = HDF5_openFile(fileName,'a') fileHandle = HDF5_openFile(fileName,'a')
call homogenization_restartWrite(fileHandle) call homogenization_restartWrite(fileHandle)
call constitutive_restartWrite(fileHandle) call phase_restartWrite(fileHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
@ -123,7 +123,7 @@ end subroutine CPFEM_restartWrite
subroutine CPFEM_forward subroutine CPFEM_forward
call homogenization_forward call homogenization_forward
call constitutive_forward call phase_forward
end subroutine CPFEM_forward end subroutine CPFEM_forward
@ -138,7 +138,7 @@ subroutine CPFEM_results(inc,time)
call results_openJobFile call results_openJobFile
call results_addIncrement(inc,time) call results_addIncrement(inc,time)
call constitutive_results call phase_results
call homogenization_results call homogenization_results
call discretization_results call discretization_results
call results_finalizeIncrement call results_finalizeIncrement

34
src/commercialFEM_fileList.f90
View File

@ -20,19 +20,19 @@
#include "material.f90" #include "material.f90"
#include "lattice.f90" #include "lattice.f90"
#include "phase.f90" #include "phase.f90"
#include "phase_mechanics.f90" #include "phase_mechanical.f90"
#include "phase_mechanics_plastic.f90" #include "phase_mechanical_plastic.f90"
#include "phase_mechanics_plastic_none.f90" #include "phase_mechanical_plastic_none.f90"
#include "phase_mechanics_plastic_isotropic.f90" #include "phase_mechanical_plastic_isotropic.f90"
#include "phase_mechanics_plastic_phenopowerlaw.f90" #include "phase_mechanical_plastic_phenopowerlaw.f90"
#include "phase_mechanics_plastic_kinehardening.f90" #include "phase_mechanical_plastic_kinehardening.f90"
#include "phase_mechanics_plastic_dislotwin.f90" #include "phase_mechanical_plastic_dislotwin.f90"
#include "phase_mechanics_plastic_dislotungsten.f90" #include "phase_mechanical_plastic_dislotungsten.f90"
#include "phase_mechanics_plastic_nonlocal.f90" #include "phase_mechanical_plastic_nonlocal.f90"
#include "phase_mechanics_eigendeformation.f90" #include "phase_mechanical_eigen.f90"
#include "phase_mechanics_eigendeformation_cleavageopening.f90" #include "phase_mechanical_eigen_cleavageopening.f90"
#include "phase_mechanics_eigendeformation_slipplaneopening.f90" #include "phase_mechanical_eigen_slipplaneopening.f90"
#include "phase_mechanics_eigendeformation_thermalexpansion.f90" #include "phase_mechanical_eigen_thermalexpansion.f90"
#include "phase_thermal.f90" #include "phase_thermal.f90"
#include "phase_thermal_dissipation.f90" #include "phase_thermal_dissipation.f90"
#include "phase_thermal_externalheat.f90" #include "phase_thermal_externalheat.f90"
@ -44,10 +44,10 @@
#include "damage_none.f90" #include "damage_none.f90"
#include "damage_nonlocal.f90" #include "damage_nonlocal.f90"
#include "homogenization.f90" #include "homogenization.f90"
#include "homogenization_mechanics.f90" #include "homogenization_mechanical.f90"
#include "homogenization_mechanics_none.f90" #include "homogenization_mechanical_pass.f90"
#include "homogenization_mechanics_isostrain.f90" #include "homogenization_mechanical_isostrain.f90"
#include "homogenization_mechanics_RGC.f90" #include "homogenization_mechanical_RGC.f90"
#include "homogenization_thermal.f90" #include "homogenization_thermal.f90"
#include "homogenization_damage.f90" #include "homogenization_damage.f90"
#include "CPFEM.f90" #include "CPFEM.f90"

66
src/grid/DAMASK_grid.f90
View File

@ -18,9 +18,9 @@ program DAMASK_grid
use CPFEM2 use CPFEM2
use material use material
use spectral_utilities use spectral_utilities
use grid_mech_spectral_basic use grid_mechanical_spectral_basic
use grid_mech_spectral_polarisation use grid_mechanical_spectral_polarisation
use grid_mech_FEM use grid_mechanical_FEM
use grid_damage_spectral use grid_damage_spectral
use grid_thermal_spectral use grid_thermal_spectral
use results use results
@ -83,16 +83,16 @@ program DAMASK_grid
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
procedure(grid_mech_spectral_basic_init), pointer :: & procedure(grid_mechanical_spectral_basic_init), pointer :: &
mech_init mechanical_init
procedure(grid_mech_spectral_basic_forward), pointer :: & procedure(grid_mechanical_spectral_basic_forward), pointer :: &
mech_forward mechanical_forward
procedure(grid_mech_spectral_basic_solution), pointer :: & procedure(grid_mechanical_spectral_basic_solution), pointer :: &
mech_solution mechanical_solution
procedure(grid_mech_spectral_basic_updateCoords), pointer :: & procedure(grid_mechanical_spectral_basic_updateCoords), pointer :: &
mech_updateCoords mechanical_updateCoords
procedure(grid_mech_spectral_basic_restartWrite), pointer :: & procedure(grid_mechanical_spectral_basic_restartWrite), pointer :: &
mech_restartWrite mechanical_restartWrite
external :: & external :: &
quit quit
@ -139,25 +139,25 @@ program DAMASK_grid
debug_grid => config_debug%get('grid',defaultVal=emptyList) debug_grid => config_debug%get('grid',defaultVal=emptyList)
select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic')))
case ('Basic') case ('Basic')
mech_init => grid_mech_spectral_basic_init mechanical_init => grid_mechanical_spectral_basic_init
mech_forward => grid_mech_spectral_basic_forward mechanical_forward => grid_mechanical_spectral_basic_forward
mech_solution => grid_mech_spectral_basic_solution mechanical_solution => grid_mechanical_spectral_basic_solution
mech_updateCoords => grid_mech_spectral_basic_updateCoords mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords
mech_restartWrite => grid_mech_spectral_basic_restartWrite mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite
case ('Polarisation') case ('Polarisation')
mech_init => grid_mech_spectral_polarisation_init mechanical_init => grid_mechanical_spectral_polarisation_init
mech_forward => grid_mech_spectral_polarisation_forward mechanical_forward => grid_mechanical_spectral_polarisation_forward
mech_solution => grid_mech_spectral_polarisation_solution mechanical_solution => grid_mechanical_spectral_polarisation_solution
mech_updateCoords => grid_mech_spectral_polarisation_updateCoords mechanical_updateCoords => grid_mechanical_spectral_polarisation_updateCoords
mech_restartWrite => grid_mech_spectral_polarisation_restartWrite mechanical_restartWrite => grid_mechanical_spectral_polarisation_restartWrite
case ('FEM') case ('FEM')
mech_init => grid_mech_FEM_init mechanical_init => grid_mechanical_FEM_init
mech_forward => grid_mech_FEM_forward mechanical_forward => grid_mechanical_FEM_forward
mech_solution => grid_mech_FEM_solution mechanical_solution => grid_mechanical_FEM_solution
mech_updateCoords => grid_mech_FEM_updateCoords mechanical_updateCoords => grid_mechanical_FEM_updateCoords
mech_restartWrite => grid_mech_FEM_restartWrite mechanical_restartWrite => grid_mechanical_FEM_restartWrite
case default case default
call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver')))
@ -303,7 +303,7 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%ID(field)) select case (loadCases(1)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mech_init call mechanical_init
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
call grid_thermal_spectral_init call grid_thermal_spectral_init
@ -379,7 +379,7 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(l)%ID(field)) select case(loadCases(l)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mech_forward (& call mechanical_forward (&
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
deformation_BC = loadCases(l)%deformation, & deformation_BC = loadCases(l)%deformation, &
stress_BC = loadCases(l)%stress, & stress_BC = loadCases(l)%stress, &
@ -399,7 +399,7 @@ program DAMASK_grid
do field = 1, nActiveFields do field = 1, nActiveFields
select case(loadCases(l)%ID(field)) select case(loadCases(l)%ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = mech_solution(incInfo) solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(timeinc) solres(field) = grid_thermal_spectral_solution(timeinc)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
@ -420,7 +420,7 @@ program DAMASK_grid
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found .and. .not. solres(1)%termIll) then ! and acceptable solution found
call mech_updateCoords call mechanical_updateCoords
timeIncOld = timeinc timeIncOld = timeinc
cutBack = .false. cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
@ -463,7 +463,7 @@ program DAMASK_grid
call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error' if (ierr /= 0) error stop 'MPI error'
if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then
call mech_restartWrite call mechanical_restartWrite
call CPFEM_restartWrite call CPFEM_restartWrite
endif endif
if(signal) call interface_setSIGUSR2(.false.) if(signal) call interface_setSIGUSR2(.false.)

135
src/grid/grid_mech_FEM.f90
View File

@ -4,7 +4,7 @@
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: FEM !> @brief Grid solver for mechanics: FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mech_FEM module grid_mechanical_FEM
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -45,8 +45,8 @@ module grid_mech_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: mech_grid DM :: mechanical_grid
SNES :: mech_snes SNES :: mechanical_snes
Vec :: solution_current, solution_lastInc, solution_rate Vec :: solution_current, solution_lastInc, solution_rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -79,18 +79,18 @@ module grid_mech_FEM
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_FEM_init, & grid_mechanical_FEM_init, &
grid_mech_FEM_solution, & grid_mechanical_FEM_solution, &
grid_mech_FEM_forward, & grid_mechanical_FEM_forward, &
grid_mech_FEM_updateCoords, & grid_mechanical_FEM_updateCoords, &
grid_mech_FEM_restartWrite grid_mechanical_FEM_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_init subroutine grid_mechanical_FEM_init
real(pReal), parameter :: HGCoeff = 0.0e-2_pReal real(pReal), parameter :: HGCoeff = 0.0e-2_pReal
real(pReal), parameter, dimension(4,8) :: & real(pReal), parameter, dimension(4,8) :: &
@ -114,7 +114,7 @@ subroutine grid_mech_FEM_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! debugging options ! debugging options
@ -141,8 +141,11 @@ subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & call PetscOptionsInsertString(PETSC_NULL_OPTIONS, &
&-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) '-mechanical_snes_type newtonls -mechanical_ksp_type fgmres &
&-mechanical_ksp_max_it 25 -mechanical_pc_type ml &
&-mechanical_mg_levels_ksp_type chebyshev', &
ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -155,8 +158,10 @@ subroutine grid_mech_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr)
CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -167,34 +172,44 @@ subroutine grid_mech_FEM_init
1, 1, worldsize, & 1, 1, worldsize, &
3, 1, & 3, 1, &
[grid(1)],[grid(2)],localK, & [grid(1)],[grid(2)],localK, &
mech_grid,ierr) mechanical_grid,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr) call SNESSetDM(mechanical_snes,mechanical_grid,ierr)
call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr)
call DMsetUp(mech_grid,ierr); CHKERRQ(ierr)
call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr) call DMsetFromOptions(mechanical_grid,ierr)
call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr)
call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) call DMsetUp(mechanical_grid,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged" call DMDASetUniformCoordinates(mechanical_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(mechanical_grid,solution_current,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(mechanical_grid,solution_lastInc,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(mechanical_grid,solution_rate ,ierr)
CHKERRQ(ierr)
call DMSNESSetFunctionLocal(mechanical_grid,formResidual,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mechanical_grid,formJacobian,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call SNESSetConvergenceTest(mechanical_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), ierr) ! ignore linear solve failures
CHKERRQ(ierr)
call SNESSetFromOptions(mechanical_snes,ierr) ! pull it all together with additional cli arguments
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr) call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr)
call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr) call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr)
call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr) call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr)
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent call DMDAGetCorners(mechanical_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent
CHKERRQ(ierr) CHKERRQ(ierr)
xend = xstart+xend-1 xend = xstart+xend-1
yend = ystart+yend-1 yend = ystart+yend-1
@ -242,9 +257,9 @@ subroutine grid_mech_FEM_init
call utilities_constitutiveResponse(P_current,P_av,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 F, & ! target F
0.0_pReal) ! time increment 0.0_pReal) ! time increment
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
restartRead2: if (interface_restartInc > 0) then restartRead2: if (interface_restartInc > 0) then
@ -257,13 +272,13 @@ subroutine grid_mech_FEM_init
endif restartRead2 endif restartRead2
end subroutine grid_mech_FEM_init end subroutine grid_mechanical_FEM_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM scheme with internal iterations !> @brief solution for the FEM scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_FEM_solution(incInfoIn) result(solution) function grid_mechanical_FEM_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -284,11 +299,13 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr); CHKERRQ(ierr) call SNESsolve(mechanical_snes,PETSC_NULL_VEC,solution_current,ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(mechanical_snes,reason,ierr)
CHKERRQ(ierr)
solution%converged = reason > 0 solution%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
@ -296,14 +313,14 @@ function grid_mech_FEM_solution(incInfoIn) result(solution)
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_FEM_solution end function grid_mechanical_FEM_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -323,8 +340,10 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
u_current,u_lastInc u_current,u_lastInc
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
if (cutBack) then if (cutBack) then
C_volAvg = C_volAvgLastInc C_volAvg = C_volAvgLastInc
@ -371,8 +390,10 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
call VecAXPY(solution_current,Delta_t,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(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide available data ! set module wide available data
@ -380,31 +401,33 @@ subroutine grid_mech_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%timeinc = Delta_t
end subroutine grid_mech_FEM_forward end subroutine grid_mechanical_FEM_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_updateCoords subroutine grid_mechanical_FEM_updateCoords
call utilities_updateCoords(F) call utilities_updateCoords(F)
end subroutine grid_mech_FEM_updateCoords end subroutine grid_mechanical_FEM_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_restartWrite subroutine grid_mechanical_FEM_restartWrite
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc
character(len=pStringLen) :: fileName character(len=pStringLen) :: fileName
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT) print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
@ -427,10 +450,12 @@ subroutine grid_mech_FEM_restartWrite
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,ierr)
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr) CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
end subroutine grid_mech_FEM_restartWrite end subroutine grid_mechanical_FEM_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -498,8 +523,10 @@ subroutine formResidual(da_local,x_local, &
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal), dimension(3,3,3,3) :: devNull real(pReal), dimension(3,3,3,3) :: devNull
call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(mechanical_snes,nfuncs,ierr)
call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) CHKERRQ(ierr)
call SNESGetIterationNumber(mechanical_snes,PETScIter,ierr)
CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
@ -679,4 +706,4 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
end subroutine formJacobian end subroutine formJacobian
end module grid_mech_FEM end module grid_mechanical_FEM

40
src/grid/grid_mech_spectral_basic.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: Spectral basic !> @brief Grid solver for mechanics: Spectral basic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mech_spectral_basic module grid_mechanical_spectral_basic
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -79,18 +79,18 @@ module grid_mech_spectral_basic
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_spectral_basic_init, & grid_mechanical_spectral_basic_init, &
grid_mech_spectral_basic_solution, & grid_mechanical_spectral_basic_solution, &
grid_mech_spectral_basic_forward, & grid_mechanical_spectral_basic_forward, &
grid_mech_spectral_basic_updateCoords, & grid_mechanical_spectral_basic_updateCoords, &
grid_mech_spectral_basic_restartWrite grid_mechanical_spectral_basic_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_init subroutine grid_mechanical_spectral_basic_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -105,7 +105,7 @@ subroutine grid_mech_spectral_basic_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_spectral_basic init -+>>>'; flush(IO_STDOUT)
print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013' print*, 'Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
@ -139,7 +139,7 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -152,7 +152,7 @@ subroutine grid_mech_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mechanical_',ierr);CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -222,13 +222,13 @@ subroutine grid_mech_spectral_basic_init
call utilities_updateGamma(C_minMaxAvg) call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness call utilities_saveReferenceStiffness
end subroutine grid_mech_spectral_basic_init end subroutine grid_mechanical_spectral_basic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the basic scheme with internal iterations !> @brief solution for the basic scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_basic_solution(incInfoIn) result(solution) function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -262,14 +262,14 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution)
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_spectral_basic_solution end function grid_mechanical_spectral_basic_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -339,13 +339,13 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%timeinc = Delta_t
end subroutine grid_mech_spectral_basic_forward end subroutine grid_mechanical_spectral_basic_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_updateCoords subroutine grid_mechanical_spectral_basic_updateCoords
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F PetscScalar, dimension(:,:,:,:), pointer :: F
@ -354,13 +354,13 @@ subroutine grid_mech_spectral_basic_updateCoords
call utilities_updateCoords(F) call utilities_updateCoords(F)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_updateCoords end subroutine grid_mechanical_spectral_basic_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_restartWrite subroutine grid_mechanical_spectral_basic_restartWrite
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
@ -393,7 +393,7 @@ subroutine grid_mech_spectral_basic_restartWrite
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_basic_restartWrite end subroutine grid_mechanical_spectral_basic_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -506,4 +506,4 @@ subroutine formResidual(in, F, &
end subroutine formResidual end subroutine formResidual
end module grid_mech_spectral_basic end module grid_mechanical_spectral_basic

40
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Grid solver for mechanics: Spectral Polarisation !> @brief Grid solver for mechanics: Spectral Polarisation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module grid_mech_spectral_polarisation module grid_mechanical_spectral_polarisation
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use PETScdmda use PETScdmda
@ -90,18 +90,18 @@ module grid_mech_spectral_polarisation
totalIter = 0 !< total iteration in current increment totalIter = 0 !< total iteration in current increment
public :: & public :: &
grid_mech_spectral_polarisation_init, & grid_mechanical_spectral_polarisation_init, &
grid_mech_spectral_polarisation_solution, & grid_mechanical_spectral_polarisation_solution, &
grid_mech_spectral_polarisation_forward, & grid_mechanical_spectral_polarisation_forward, &
grid_mech_spectral_polarisation_updateCoords, & grid_mechanical_spectral_polarisation_updateCoords, &
grid_mech_spectral_polarisation_restartWrite grid_mechanical_spectral_polarisation_restartWrite
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_init subroutine grid_mechanical_spectral_polarisation_init
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -118,7 +118,7 @@ subroutine grid_mech_spectral_polarisation_init
num_grid, & num_grid, &
debug_grid debug_grid
print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- grid_mechanical_spectral_polarisation init -+>>>'; flush(IO_STDOUT)
print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
@ -157,7 +157,7 @@ subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -172,7 +172,7 @@ subroutine grid_mech_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mechanical_',ierr);CHKERRQ(ierr)
localK = 0 localK = 0
localK(worldrank) = grid3 localK(worldrank) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@ -250,13 +250,13 @@ subroutine grid_mech_spectral_polarisation_init
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
end subroutine grid_mech_spectral_polarisation_init end subroutine grid_mechanical_spectral_polarisation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
@ -294,14 +294,14 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution)
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_aim,P_av,params%stress_mask)
end function grid_mech_spectral_polarisation_solution end function grid_mechanical_spectral_polarisation_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!> @details find new boundary conditions and best F estimate for end of current timestep !> @details find new boundary conditions and best F estimate for end of current timestep
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,& subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t_old,t_remaining,&
deformation_BC,stress_BC,rotation_BC) deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: & logical, intent(in) :: &
@ -393,13 +393,13 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%timeinc = Delta_t
end subroutine grid_mech_spectral_polarisation_forward end subroutine grid_mechanical_spectral_polarisation_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Update coordinates !> @brief Update coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_updateCoords subroutine grid_mechanical_spectral_polarisation_updateCoords
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau
@ -408,13 +408,13 @@ subroutine grid_mech_spectral_polarisation_updateCoords
call utilities_updateCoords(FandF_tau(0:8,:,:,:)) call utilities_updateCoords(FandF_tau(0:8,:,:,:))
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_updateCoords end subroutine grid_mechanical_spectral_polarisation_updateCoords
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file !> @brief Write current solver and constitutive data for restart to file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_restartWrite subroutine grid_mechanical_spectral_polarisation_restartWrite
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
@ -450,7 +450,7 @@ subroutine grid_mech_spectral_polarisation_restartWrite
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_restartWrite end subroutine grid_mechanical_spectral_polarisation_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -618,4 +618,4 @@ subroutine formResidual(in, FandF_tau, &
end subroutine formResidual end subroutine formResidual
end module grid_mech_spectral_polarisation end module grid_mechanical_spectral_polarisation

32
src/homogenization.f90
View File

@ -45,10 +45,10 @@ module homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
interface interface
module subroutine mech_init(num_homog) module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homog !< pointer to mechanical homogenization numerics data num_homog !< pointer to mechanical homogenization numerics data
end subroutine mech_init end subroutine mechanical_init
module subroutine thermal_init module subroutine thermal_init
end subroutine thermal_init end subroutine thermal_init
@ -56,13 +56,13 @@ module homogenization
module subroutine damage_init module subroutine damage_init
end subroutine damage_init end subroutine damage_init
module subroutine mech_partition(subF,ip,el) module subroutine mechanical_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
subF subF
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
end subroutine mech_partition end subroutine mechanical_partition
module subroutine thermal_partition(ce) module subroutine thermal_partition(ce)
integer, intent(in) :: ce integer, intent(in) :: ce
@ -76,19 +76,19 @@ module homogenization
integer, intent(in) :: ip,el integer, intent(in) :: ip,el
end subroutine thermal_homogenize end subroutine thermal_homogenize
module subroutine mech_homogenize(dt,ip,el) module subroutine mechanical_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
end subroutine mech_homogenize end subroutine mechanical_homogenize
module subroutine mech_results(group_base,h) module subroutine mechanical_results(group_base,h)
character(len=*), intent(in) :: group_base character(len=*), intent(in) :: group_base
integer, intent(in) :: h integer, intent(in) :: h
end subroutine mech_results end subroutine mechanical_results
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< current time step subdt !< current time step
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -97,7 +97,7 @@ module homogenization
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
end function mech_updateState end function mechanical_updateState
module function thermal_conduction_getConductivity(ip,el) result(K) module function thermal_conduction_getConductivity(ip,el) result(K)
@ -216,7 +216,7 @@ subroutine homogenization_init()
if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
call mech_init(num_homog) call mechanical_init(num_homog)
call thermal_init() call thermal_init()
call damage_init() call damage_init()
@ -253,7 +253,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
me = material_homogenizationMemberAt2(ce) me = material_homogenizationMemberAt2(ce)
call constitutive_restore(ce,.false.) ! wrong name (is more a forward function) call phase_restore(ce,.false.) ! wrong name (is more a forward function)
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me) if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me) if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me)
@ -267,7 +267,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (.not. doneAndHappy(1)) then if (.not. doneAndHappy(1)) then
call mech_partition(homogenization_F(1:3,1:3,ce),ip,el) call mechanical_partition(homogenization_F(1:3,1:3,ce),ip,el)
converged = .true. converged = .true.
do co = 1, myNgrains do co = 1, myNgrains
converged = converged .and. crystallite_stress(dt,co,ip,el) converged = converged .and. crystallite_stress(dt,co,ip,el)
@ -276,7 +276,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
if (.not. converged) then if (.not. converged) then
doneAndHappy = [.true.,.false.] doneAndHappy = [.true.,.false.]
else else
doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el)
converged = all(doneAndHappy) converged = all(doneAndHappy)
endif endif
endif endif
@ -338,7 +338,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
enddo enddo
call mech_homogenize(dt,ip,el) call mechanical_homogenize(dt,ip,el)
enddo IpLooping3 enddo IpLooping3
enddo elementLooping3 enddo elementLooping3
!$OMP END DO !$OMP END DO
@ -365,7 +365,7 @@ subroutine homogenization_results
group_base = 'current/homogenization/'//trim(material_name_homogenization(ho)) group_base = 'current/homogenization/'//trim(material_name_homogenization(ho))
call results_closeGroup(results_addGroup(group_base)) call results_closeGroup(results_addGroup(group_base))
call mech_results(group_base,ho) call mechanical_results(group_base,ho)
group = trim(group_base)//'/damage' group = trim(group_base)//'/damage'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))

4
src/homogenization_damage.f90
View File

@ -74,7 +74,7 @@ module subroutine damage_partition(ce)
phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce)) phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_damage_set_phi(phi,co,ce) call phase_damage_set_phi(phi,co,ce)
enddo enddo
end subroutine damage_partition end subroutine damage_partition
@ -120,7 +120,7 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)

94
src/homogenization_mechanics.f90 → src/homogenization_mechanical.f90
View File

@ -7,52 +7,52 @@ submodule(homogenization) mechanics
interface interface
module subroutine mech_none_init module subroutine mechanical_pass_init
end subroutine mech_none_init end subroutine mechanical_pass_init
module subroutine mech_isostrain_init module subroutine mechanical_isostrain_init
end subroutine mech_isostrain_init end subroutine mechanical_isostrain_init
module subroutine mech_RGC_init(num_homogMech) module subroutine mechanical_RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
end subroutine mech_RGC_init end subroutine mechanical_RGC_init
module subroutine mech_isostrain_partitionDeformation(F,avgF) module subroutine mechanical_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
end subroutine mech_isostrain_partitionDeformation end subroutine mechanical_isostrain_partitionDeformation
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
of of
end subroutine mech_RGC_partitionDeformation end subroutine mechanical_RGC_partitionDeformation
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance integer, intent(in) :: instance
end subroutine mech_isostrain_averageStressAndItsTangent end subroutine mechanical_isostrain_averageStressAndItsTangent
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance integer, intent(in) :: instance
end subroutine mech_RGC_averageStressAndItsTangent end subroutine mechanical_RGC_averageStressAndItsTangent
module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
@ -63,13 +63,13 @@ submodule(homogenization) mechanics
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
end function mech_RGC_updateState end function mechanical_RGC_updateState
module subroutine mech_RGC_results(instance,group) module subroutine mechanical_RGC_results(instance,group)
integer, intent(in) :: instance !< homogenization instance integer, intent(in) :: instance !< homogenization instance
character(len=*), intent(in) :: group !< group name in HDF5 file character(len=*), intent(in) :: group !< group name in HDF5 file
end subroutine mech_RGC_results end subroutine mechanical_RGC_results
end interface end interface
@ -78,7 +78,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters. !> @brief Allocate variables and set parameters.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_init(num_homog) module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homog num_homog
@ -94,17 +94,17 @@ module subroutine mech_init(num_homog)
allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal) allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict) num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mechanical_pass_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mechanical_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech) if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mechanical_RGC_init(num_homogMech)
end subroutine mech_init end subroutine mechanical_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Partition F onto the individual constituents. !> @brief Partition F onto the individual constituents.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_partition(subF,ip,el) module subroutine mechanical_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
subF subF
@ -122,25 +122,25 @@ module subroutine mech_partition(subF,ip,el)
Fs(1:3,1:3,1) = subF Fs(1:3,1:3,1) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(Fs,subF) call mechanical_isostrain_partitionDeformation(Fs,subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(Fs,subF,ip,el) call mechanical_RGC_partitionDeformation(Fs,subF,ip,el)
end select chosenHomogenization end select chosenHomogenization
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
call constitutive_mech_setF(Fs(1:3,1:3,co),co,ip,el) call phase_mechanical_setF(Fs(1:3,1:3,co),co,ip,el)
enddo enddo
end subroutine mech_partition end subroutine mechanical_partition
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Average P and dPdF from the individual constituents. !> @brief Average P and dPdF from the individual constituents.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_homogenize(dt,ip,el) module subroutine mechanical_homogenize(dt,ip,el)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
@ -156,15 +156,15 @@ module subroutine mech_homogenize(dt,ip,el)
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el) homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el) Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
enddo enddo
call mech_isostrain_averageStressAndItsTangent(& call mechanical_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
@ -172,10 +172,10 @@ module subroutine mech_homogenize(dt,ip,el)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el) Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
enddo enddo
call mech_RGC_averageStressAndItsTangent(& call mechanical_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
@ -183,14 +183,14 @@ module subroutine mech_homogenize(dt,ip,el)
end select chosenHomogenization end select chosenHomogenization
end subroutine mech_homogenize end subroutine mechanical_homogenize
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result !> "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< current time step subdt !< current time step
@ -209,22 +209,22 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el)
Fs(:,:,co) = constitutive_mech_getF(co,ip,el) Fs(:,:,co) = phase_mechanical_getF(co,ip,el)
Ps(:,:,co) = constitutive_mech_getP(co,ip,el) Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
enddo enddo
doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el)
else else
doneAndHappy = .true. doneAndHappy = .true.
endif endif
end function mech_updateState end function mechanical_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Write results to file. !> @brief Write results to file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_results(group_base,h) module subroutine mechanical_results(group_base,h)
use material, only: & use material, only: &
material_homogenization_type => homogenization_type material_homogenization_type => homogenization_type
@ -239,7 +239,7 @@ module subroutine mech_results(group_base,h)
select case(material_homogenization_type(h)) select case(material_homogenization_type(h))
case(HOMOGENIZATION_rgc_ID) case(HOMOGENIZATION_rgc_ID)
call mech_RGC_results(homogenization_typeInstance(h),group) call mechanical_RGC_results(homogenization_typeInstance(h),group)
end select end select
@ -250,7 +250,7 @@ module subroutine mech_results(group_base,h)
!call results_writeDataset(group,temp,'P',& !call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa') ! '1st Piola-Kirchhoff stress','Pa')
end subroutine mech_results end subroutine mechanical_results
end submodule mechanics end submodule mechanics

24
src/homogenization_mechanics_RGC.f90 → src/homogenization_mechanical_RGC.f90
View File

@ -71,7 +71,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_init(num_homogMech) module subroutine mechanical_RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: & class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
@ -155,7 +155,7 @@ module subroutine mech_RGC_init(num_homogMech)
prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) & if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='N_constituents (mech_RGC)') call IO_error(211,ext_msg='N_constituents (mechanical_RGC)')
prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
prm%c_alpha = homogMech%get_asFloat('c_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha')
@ -190,13 +190,13 @@ module subroutine mech_RGC_init(num_homogMech)
enddo enddo
end subroutine mech_RGC_init end subroutine mechanical_RGC_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain
@ -229,14 +229,14 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
end associate end associate
end subroutine mech_RGC_partitionDeformation end subroutine mechanical_RGC_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result ! "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy) module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
real(pReal), dimension(:,:,:), intent(in) :: & real(pReal), dimension(:,:,:), intent(in) :: &
P,& !< partitioned stresses P,& !< partitioned stresses
@ -658,7 +658,7 @@ module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
C = constitutive_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el)) C = phase_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el))
equivalentMu = lattice_equivalent_mu(C,'voigt') equivalentMu = lattice_equivalent_mu(C,'voigt')
end function equivalentMu end function equivalentMu
@ -704,13 +704,13 @@ module function mech_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy
end subroutine grainDeformation end subroutine grainDeformation
end function mech_RGC_updateState end function mechanical_RGC_updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -722,13 +722,13 @@ module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ins
avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal) avgP = sum(P,3) /real(product(param(instance)%N_constituents),pReal)
dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal)
end subroutine mech_RGC_averageStressAndItsTangent end subroutine mechanical_RGC_averageStressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_RGC_results(instance,group) module subroutine mechanical_RGC_results(instance,group)
integer, intent(in) :: instance integer, intent(in) :: instance
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -754,7 +754,7 @@ module subroutine mech_RGC_results(instance,group)
enddo outputsLoop enddo outputsLoop
end associate end associate
end subroutine mech_RGC_results end subroutine mechanical_RGC_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

14
src/homogenization_mechanics_isostrain.f90 → src/homogenization_mechanical_isostrain.f90
View File

@ -26,7 +26,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file !> @brief allocates all neccessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_init module subroutine mechanical_isostrain_init
integer :: & integer :: &
Ninstances, & Ninstances, &
@ -58,7 +58,7 @@ module subroutine mech_isostrain_init
case ('avg') case ('avg')
prm%mapping = average_ID prm%mapping = average_ID
case default case default
call IO_error(211,ext_msg='sum'//' (mech_isostrain)') call IO_error(211,ext_msg='sum'//' (mechanical_isostrain)')
end select end select
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
@ -70,13 +70,13 @@ module subroutine mech_isostrain_init
enddo enddo
end subroutine mech_isostrain_init end subroutine mechanical_isostrain_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents !> @brief partitions the deformation gradient onto the constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_partitionDeformation(F,avgF) module subroutine mechanical_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
@ -84,13 +84,13 @@ module subroutine mech_isostrain_partitionDeformation(F,avgF)
F = spread(avgF,3,size(F,3)) F = spread(avgF,3,size(F,3))
end subroutine mech_isostrain_partitionDeformation end subroutine mechanical_isostrain_partitionDeformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities !> @brief derive average stress and stiffness from constituent quantities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance) module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
@ -112,6 +112,6 @@ module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP
end associate end associate
end subroutine mech_isostrain_averageStressAndItsTangent end subroutine mechanical_isostrain_averageStressAndItsTangent
end submodule isostrain end submodule isostrain

8
src/homogenization_mechanics_none.f90 → src/homogenization_mechanical_pass.f90
View File

@ -11,14 +11,14 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file !> @brief allocates all necessary fields, reads information from material configuration file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_none_init module subroutine mechanical_pass_init
integer :: & integer :: &
Ninstances, & Ninstances, &
h, & h, &
Nmaterialpoints Nmaterialpoints
print'(/,a)', ' <<<+- homogenization:mechanics:none init -+>>>' print'(/,a)', ' <<<+- homogenization:mechanics:pass init -+>>>'
Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
@ -27,7 +27,7 @@ module subroutine mech_none_init
if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Nconstituents(h) /= 1) & if(homogenization_Nconstituents(h) /= 1) &
call IO_error(211,ext_msg='N_constituents (mech_none)') call IO_error(211,ext_msg='N_constituents (mechanical_pass)')
Nmaterialpoints = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
@ -36,6 +36,6 @@ module subroutine mech_none_init
enddo enddo
end subroutine mech_none_init end subroutine mechanical_pass_init
end submodule none end submodule none

6
src/homogenization_thermal.f90
View File

@ -78,7 +78,7 @@ module subroutine thermal_partition(ce)
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce))
dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
call constitutive_thermal_setField(T,dot_T,co,ce) call phase_thermal_setField(T,dot_T,co,ce)
enddo enddo
end subroutine thermal_partition end subroutine thermal_partition
@ -91,7 +91,7 @@ module subroutine thermal_homogenize(ip,el)
integer, intent(in) :: ip,el integer, intent(in) :: ip,el
!call constitutive_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el) !call phase_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el)
end subroutine thermal_homogenize end subroutine thermal_homogenize
@ -235,7 +235,7 @@ module subroutine thermal_conduction_getSource(Tdot, ip,el)
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
call constitutive_thermal_getRate(dot_T_temp, ph,me) call phase_thermal_getRate(dot_T_temp, ph,me)
Tdot = Tdot + dot_T_temp Tdot = Tdot + dot_T_temp
enddo enddo

8
src/mesh/DAMASK_mesh.f90
View File

@ -18,7 +18,7 @@ program DAMASK_mesh
use config use config
use discretization_mesh use discretization_mesh
use FEM_Utilities use FEM_Utilities
use mesh_mech_FEM use mesh_mechanical_FEM
implicit none implicit none
@ -242,7 +242,7 @@ program DAMASK_mesh
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(1)%fieldBC(field)%ID) select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call FEM_mech_init(loadCases(1)%fieldBC(field)) call FEM_mechanical_init(loadCases(1)%fieldBC(field))
end select end select
enddo enddo
@ -306,7 +306,7 @@ program DAMASK_mesh
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call FEM_mech_forward (& call FEM_mechanical_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select end select
@ -320,7 +320,7 @@ program DAMASK_mesh
do field = 1, nActiveFields do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = FEM_mech_solution (& solres(field) = FEM_mechanical_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select end select

12
src/mesh/FEM_utilities.f90
View File

@ -127,12 +127,12 @@ subroutine FEM_utilities_init
CHKERRQ(ierr) CHKERRQ(ierr)
if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls & call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls &
&-mech_snes_linesearch_type cp -mech_snes_ksp_ew & &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew &
&-mech_snes_ksp_ew_rtol0 0.01 -mech_snes_ksp_ew_rtolmax 0.01 & &-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 &
&-mech_ksp_type fgmres -mech_ksp_max_it 25 & &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25 &
&-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev & &-mechanical_pc_type ml -mechanical_mg_levels_ksp_type chebyshev &
&-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr) &-mechanical_mg_levels_pc_type sor -mechanical_pc_ml_nullspace user',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('petsc_options',defaultVal=''),ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('petsc_options',defaultVal=''),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

110
src/mesh/mesh_mech_FEM.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief FEM PETSc solver !> @brief FEM PETSc solver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module mesh_mech_FEM module mesh_mechanical_FEM
#include <petsc/finclude/petscdmplex.h> #include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdm.h> #include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petsc.h> #include <petsc/finclude/petsc.h>
@ -50,7 +50,7 @@ module mesh_mech_FEM
type(tNumerics), private :: num type(tNumerics), private :: num
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
SNES :: mech_snes SNES :: mechanical_snes
Vec :: solution, solution_rate, solution_local Vec :: solution, solution_rate, solution_local
PetscInt :: dimPlex, cellDof, nQuadrature, nBasis PetscInt :: dimPlex, cellDof, nQuadrature, nBasis
PetscReal, allocatable, target :: qPoints(:), qWeights(:) PetscReal, allocatable, target :: qPoints(:), qWeights(:)
@ -65,20 +65,20 @@ module mesh_mech_FEM
real(pReal), parameter :: eps = 1.0e-18_pReal real(pReal), parameter :: eps = 1.0e-18_pReal
public :: & public :: &
FEM_mech_init, & FEM_mechanical_init, &
FEM_mech_solution, & FEM_mechanical_solution, &
FEM_mech_forward FEM_mechanical_forward
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data !> @brief allocates all neccessary fields and fills them with data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_init(fieldBC) subroutine FEM_mechanical_init(fieldBC)
type(tFieldBC), intent(in) :: fieldBC type(tFieldBC), intent(in) :: fieldBC
DM :: mech_mesh DM :: mechanical_mesh
PetscFE :: mechFE PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS PetscDS :: mechDS
@ -126,8 +126,8 @@ subroutine FEM_mech_init(fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh ! Setup FEM mech mesh
call DMClone(geomMesh,mech_mesh,ierr); CHKERRQ(ierr) call DMClone(geomMesh,mechanical_mesh,ierr); CHKERRQ(ierr)
call DMGetDimension(mech_mesh,dimPlex,ierr); CHKERRQ(ierr) call DMGetDimension(mechanical_mesh,dimPlex,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization ! Setup FEM mech discretization
@ -146,22 +146,22 @@ subroutine FEM_mech_init(fieldBC)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc nBasis = nBasis/nc
call DMAddField(mech_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr) call DMAddField(mechanical_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr)
call DMCreateDS(mech_mesh,ierr); CHKERRQ(ierr) call DMCreateDS(mechanical_mesh,ierr); CHKERRQ(ierr)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mechanical_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr) call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr) call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Setup FEM mech boundary conditions ! Setup FEM mech boundary conditions
call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) call DMGetLabel(mechanical_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) call DMPlexLabelComplete(mechanical_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMGetLocalSection(mechanical_mesh,section,ierr); CHKERRQ(ierr)
allocate(pnumComp(1), source=dimPlex) allocate(pnumComp(1), source=dimPlex)
allocate(pnumDof(0:dimPlex), source = 0) allocate(pnumDof(0:dimPlex), source = 0)
do topologDim = 0, dimPlex do topologDim = 0, dimPlex
call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr) call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr) call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
@ -179,10 +179,10 @@ subroutine FEM_mech_init(fieldBC)
numBC = numBC + 1 numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr) call DMGetStratumIS(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr) call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
@ -195,32 +195,32 @@ subroutine FEM_mech_init(fieldBC)
endif endif
endif endif
enddo; enddo enddo; enddo
call DMPlexCreateSection(mech_mesh,nolabel,pNumComp,pNumDof, & call DMPlexCreateSection(mechanical_mesh,nolabel,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr) numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) call DMSetSection(mechanical_mesh,section,ierr); CHKERRQ(ierr)
do faceSet = 1, numBC do faceSet = 1, numBC
call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr) call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr);CHKERRQ(ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr);CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver call SNESSetDM(mechanical_snes,mechanical_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver
call DMCreateGlobalVector(mech_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs call DMCreateGlobalVector(mechanical_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs
call DMCreateGlobalVector(mech_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step call DMCreateGlobalVector(mechanical_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMCreateLocalVector (mech_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step call DMCreateLocalVector (mechanical_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMSNESSetFunctionLocal(mech_mesh,FEM_mech_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces call DMSNESSetFunctionLocal(mechanical_mesh,FEM_mechanical_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix call DMSNESSetJacobianLocal(mechanical_mesh,FEM_mechanical_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr) call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetTolerances(mech_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr) call SNESSetTolerances(mechanical_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) call SNESSetFromOptions(mechanical_snes,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
@ -236,11 +236,11 @@ subroutine FEM_mech_init(fieldBC)
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr) call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr) call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) call DMPlexGetHeightStratum(mechanical_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0_pReal x_scal = 0.0_pReal
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(mechanical_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis*dimPlex-1, dimPlex do basis = 0, nBasis*dimPlex-1, dimPlex
@ -251,17 +251,17 @@ subroutine FEM_mech_init(fieldBC)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
enddo enddo
px_scal => x_scal px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,5,ierr) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
end subroutine FEM_mech_init end subroutine FEM_mechanical_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM load step !> @brief solution for the FEM load step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function FEM_mech_solution( & type(tSolutionState) function FEM_mechanical_solution( &
incInfoIn,timeinc,timeinc_old,fieldBC) incInfoIn,timeinc,timeinc_old,fieldBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -278,35 +278,35 @@ type(tSolutionState) function FEM_mech_solution( &
SNESConvergedReason :: reason SNESConvergedReason :: reason
incInfo = incInfoIn incInfo = incInfoIn
FEM_mech_solution%converged =.false. FEM_mechanical_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%fieldBC = fieldBC params%fieldBC = fieldBC
call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution) call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false. terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mech_solution%converged = .false. FEM_mechanical_solution%converged = .false.
FEM_mech_solution%iterationsNeeded = num%itmax FEM_mechanical_solution%iterationsNeeded = num%itmax
else ! >= 1 proper convergence (or terminally ill) else ! >= 1 proper convergence (or terminally ill)
FEM_mech_solution%converged = .true. FEM_mechanical_solution%converged = .true.
call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr) call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
print'(/,a)', ' ===========================================================================' print'(/,a)', ' ==========================================================================='
flush(IO_STDOUT) flush(IO_STDOUT)
end function FEM_mech_solution end function FEM_mechanical_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the FEM residual vector !> @brief forms the FEM residual vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
DM :: dm_local DM :: dm_local
PetscObject,intent(in) :: dummy PetscObject,intent(in) :: dummy
@ -431,13 +431,13 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formResidual end subroutine FEM_mechanical_formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forms the FEM stiffness matrix !> @brief forms the FEM stiffness matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr) subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
DM :: dm_local DM :: dm_local
@ -575,13 +575,13 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formJacobian end subroutine FEM_mechanical_formJacobian
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
type(tFieldBC), intent(in) :: & type(tFieldBC), intent(in) :: &
fieldBC fieldBC
@ -603,7 +603,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
if (guess .and. .not. cutBack) then if (guess .and. .not. cutBack) then
ForwardData = .True. ForwardData = .True.
homogenization_F0 = homogenization_F homogenization_F0 = homogenization_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mechanical_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr) call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr)
@ -634,13 +634,13 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr) call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
call VecScale(solution,timeinc,ierr); CHKERRQ(ierr) call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_forward end subroutine FEM_mechanical_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reporting !> @brief reporting
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr) subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
SNES :: snes_local SNES :: snes_local
PetscInt :: PETScIter PetscInt :: PETScIter
@ -662,6 +662,6 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal ' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)
end subroutine FEM_mech_converged end subroutine FEM_mechanical_converged
end module mesh_mech_FEM end module mesh_mechanical_FEM

180
src/phase.f90
View File

@ -62,7 +62,7 @@ module phase
phase_Nsources, & !< number of source mechanisms active in each phase phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_plasticityInstance, & !< instance of particular plasticity of each phase phase_plasticInstance, & !< instance of particular plasticity of each phase
phase_elasticityInstance !< instance of particular elasticity of each phase phase_elasticityInstance !< instance of particular elasticity of each phase
logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler)
@ -75,15 +75,15 @@ module phase
integer, public, protected :: & integer, public, protected :: &
constitutive_plasticity_maxSizeDotState, & phase_plasticity_maxSizeDotState, &
constitutive_source_maxSizeDotState phase_source_maxSizeDotState
interface interface
! == cleaned:begin ================================================================================= ! == cleaned:begin =================================================================================
module subroutine mech_init(phases) module subroutine mechanical_init(phases)
class(tNode), pointer :: phases class(tNode), pointer :: phases
end subroutine mech_init end subroutine mechanical_init
module subroutine damage_init module subroutine damage_init
end subroutine damage_init end subroutine damage_init
@ -93,83 +93,83 @@ module phase
end subroutine thermal_init end subroutine thermal_init
module subroutine mech_results(group,ph) module subroutine mechanical_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine mech_results end subroutine mechanical_results
module subroutine damage_results(group,ph) module subroutine damage_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine damage_results end subroutine damage_results
module subroutine mech_windForward(ph,me) module subroutine mechanical_windForward(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
end subroutine mech_windForward end subroutine mechanical_windForward
module subroutine mech_forward() module subroutine mechanical_forward()
end subroutine mech_forward end subroutine mechanical_forward
module subroutine thermal_forward() module subroutine thermal_forward()
end subroutine thermal_forward end subroutine thermal_forward
module subroutine mech_restore(ce,includeL) module subroutine mechanical_restore(ce,includeL)
integer, intent(in) :: ce integer, intent(in) :: ce
logical, intent(in) :: includeL logical, intent(in) :: includeL
end subroutine mech_restore end subroutine mechanical_restore
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
co, & !< counter in constituent loop co, & !< counter in constituent loop
ip, & !< counter in integration point loop ip, & !< counter in integration point loop
el !< counter in element loop el !< counter in element loop
real(pReal), dimension(3,3,3,3) :: dPdF real(pReal), dimension(3,3,3,3) :: dPdF
end function constitutive_mech_dPdF end function phase_mechanical_dPdF
module subroutine mech_restartWrite(groupHandle,ph) module subroutine mechanical_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine mech_restartWrite end subroutine mechanical_restartWrite
module subroutine mech_restartRead(groupHandle,ph) module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
end subroutine mech_restartRead end subroutine mechanical_restartRead
module function mech_S(ph,me) result(S) module function mechanical_S(ph,me) result(S)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: S real(pReal), dimension(3,3) :: S
end function mech_S end function mechanical_S
module function mech_L_p(ph,me) result(L_p) module function mechanical_L_p(ph,me) result(L_p)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
end function mech_L_p end function mechanical_L_p
module function constitutive_mech_getF(co,ip,el) result(F) module function phase_mechanical_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F real(pReal), dimension(3,3) :: F
end function constitutive_mech_getF end function phase_mechanical_getF
module function mech_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
end function mech_F_e end function mechanical_F_e
module function constitutive_mech_getP(co,ip,el) result(P) module function phase_mechanical_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
end function constitutive_mech_getP end function phase_mechanical_getP
module function constitutive_damage_get_phi(co,ip,el) result(phi) module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal) :: phi real(pReal) :: phi
end function constitutive_damage_get_phi end function phase_damage_get_phi
module function thermal_T(ph,me) result(T) module function thermal_T(ph,me) result(T)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
@ -182,20 +182,20 @@ module phase
end function thermal_dot_T end function thermal_dot_T
module subroutine constitutive_mech_setF(F,co,ip,el) module subroutine phase_mechanical_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
end subroutine constitutive_mech_setF end subroutine phase_mechanical_setF
module subroutine constitutive_thermal_setField(T,dot_T, co,ce) module subroutine phase_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
end subroutine constitutive_thermal_setField end subroutine phase_thermal_setField
module subroutine constitutive_damage_set_phi(phi,co,ce) module subroutine phase_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
end subroutine constitutive_damage_set_phi end subroutine phase_damage_set_phi
! == cleaned:end =================================================================================== ! == cleaned:end ===================================================================================
@ -222,13 +222,13 @@ module phase
logical :: converged_ logical :: converged_
end function crystallite_stress end function crystallite_stress
module function constitutive_homogenizedC(ph,me) result(C) module function phase_homogenizedC(ph,me) result(C)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
end function constitutive_homogenizedC end function phase_homogenizedC
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
@ -237,13 +237,13 @@ module phase
real(pReal), intent(inout) :: & real(pReal), intent(inout) :: &
phiDot, & phiDot, &
dPhiDot_dPhi dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents end subroutine phase_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRate(TDot, ph,me) module subroutine phase_thermal_getRate(TDot, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
TDot TDot
end subroutine constitutive_thermal_getRate end subroutine phase_thermal_getRate
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
integer, intent(in) :: & integer, intent(in) :: &
@ -281,39 +281,39 @@ module phase
#endif #endif
public :: & public :: &
constitutive_init, & phase_init, &
constitutive_homogenizedC, & phase_homogenizedC, &
constitutive_damage_getRateAndItsTangents, & phase_damage_getRateAndItsTangents, &
constitutive_thermal_getRate, & phase_thermal_getRate, &
constitutive_results, & phase_results, &
constitutive_allocateState, & phase_allocateState, &
constitutive_forward, & phase_forward, &
constitutive_restore, & phase_restore, &
plastic_nonlocal_updateCompatibility, & plastic_nonlocal_updateCompatibility, &
converged, & converged, &
crystallite_init, & crystallite_init, &
crystallite_stress, & crystallite_stress, &
thermal_stress, & thermal_stress, &
constitutive_mech_dPdF, & phase_mechanical_dPdF, &
crystallite_orientations, & crystallite_orientations, &
crystallite_push33ToRef, & crystallite_push33ToRef, &
constitutive_restartWrite, & phase_restartWrite, &
constitutive_restartRead, & phase_restartRead, &
integrateDamageState, & integrateDamageState, &
constitutive_thermal_setField, & phase_thermal_setField, &
constitutive_damage_set_phi, & phase_damage_set_phi, &
constitutive_damage_get_phi, & phase_damage_get_phi, &
constitutive_mech_getP, & phase_mechanical_getP, &
constitutive_mech_setF, & phase_mechanical_setF, &
constitutive_mech_getF, & phase_mechanical_getF, &
constitutive_windForward phase_windForward
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics !> @brief Initialze constitutive models for individual physics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_init subroutine phase_init
integer :: & integer :: &
ph, & !< counter in phase loop ph, & !< counter in phase loop
@ -336,12 +336,12 @@ subroutine constitutive_init
phases => config_material%get('phase') phases => config_material%get('phase')
call mech_init(phases) call mechanical_init(phases)
call damage_init call damage_init
call thermal_init(phases) call thermal_init(phases)
constitutive_source_maxSizeDotState = 0 phase_source_maxSizeDotState = 0
PhaseLoop2:do ph = 1,phases%length PhaseLoop2:do ph = 1,phases%length
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! partition and initialize state ! partition and initialize state
@ -350,18 +350,18 @@ subroutine constitutive_init
damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0
end forall end forall
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, &
maxval(damageState(ph)%p%sizeDotState)) maxval(damageState(ph)%p%sizeDotState))
enddo PhaseLoop2 enddo PhaseLoop2
constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) phase_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)
end subroutine constitutive_init end subroutine phase_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState(state, & subroutine phase_allocateState(state, &
Nconstituents,sizeState,sizeDotState,sizeDeltaState) Nconstituents,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: & class(tState), intent(out) :: &
@ -387,13 +387,13 @@ subroutine constitutive_allocateState(state, &
allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal)
end subroutine constitutive_allocateState end subroutine phase_allocateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback. !> @brief Restore data after homog cutback.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_restore(ce,includeL) subroutine phase_restore(ce,includeL)
logical, intent(in) :: includeL logical, intent(in) :: includeL
integer, intent(in) :: ce integer, intent(in) :: ce
@ -410,21 +410,21 @@ subroutine constitutive_restore(ce,includeL)
enddo enddo
enddo enddo
call mech_restore(ce,includeL) call mechanical_restore(ce,includeL)
end subroutine constitutive_restore end subroutine phase_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible? ! ToDo: Any guessing for the current states possible?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_forward() subroutine phase_forward()
integer :: ph, so integer :: ph, so
call mech_forward() call mechanical_forward()
call thermal_forward() call thermal_forward()
do ph = 1, size(damageState) do ph = 1, size(damageState)
@ -432,13 +432,13 @@ subroutine constitutive_forward()
damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state
enddo; enddo enddo; enddo
end subroutine constitutive_forward end subroutine phase_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file !> @brief writes constitutive results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_results() subroutine phase_results()
integer :: ph integer :: ph
character(len=:), allocatable :: group character(len=:), allocatable :: group
@ -451,12 +451,12 @@ subroutine constitutive_results()
group = '/current/phase/'//trim(material_name_phase(ph))//'/' group = '/current/phase/'//trim(material_name_phase(ph))//'/'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
call mech_results(group,ph) call mechanical_results(group,ph)
call damage_results(group,ph) call damage_results(group,ph)
enddo enddo
end subroutine constitutive_results end subroutine phase_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -557,7 +557,7 @@ end subroutine crystallite_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward. !> @brief Wind homog inc forward.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_windForward(ip,el) subroutine phase_windForward(ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -572,7 +572,7 @@ subroutine constitutive_windForward(ip,el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call mech_windForward(ph,me) call mechanical_windForward(ph,me)
do so = 1, phase_Nsources(material_phaseAt(co,el)) do so = 1, phase_Nsources(material_phaseAt(co,el))
damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me)
@ -580,7 +580,7 @@ subroutine constitutive_windForward(ip,el)
enddo enddo
end subroutine constitutive_windForward end subroutine phase_windForward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -595,11 +595,11 @@ subroutine crystallite_orientations(co,ip,el)
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&
mech_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))))) mechanical_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))))
if (plasticState(material_phaseAt(1,el))%nonlocal) & if (plasticState(material_phaseAt(1,el))%nonlocal) &
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(1,el)),ip,el) phase_plasticInstance(material_phaseAt(1,el)),ip,el)
end subroutine crystallite_orientations end subroutine crystallite_orientations
@ -620,7 +620,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
real(pReal), dimension(3,3) :: T real(pReal), dimension(3,3) :: T
T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(constitutive_mech_getF(co,ip,el)))) ! ToDo: initial orientation correct? T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ip,el)))) ! ToDo: initial orientation correct?
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
@ -648,7 +648,7 @@ end function converged
!> @brief Write current restart information (Field and constitutive data) to file. !> @brief Write current restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI ! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_restartWrite(fileHandle) subroutine phase_restartWrite(fileHandle)
integer(HID_T), intent(in) :: fileHandle integer(HID_T), intent(in) :: fileHandle
@ -662,7 +662,7 @@ subroutine constitutive_restartWrite(fileHandle)
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph)) groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph))
call mech_restartWrite(groupHandle(2),ph) call mechanical_restartWrite(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
@ -670,14 +670,14 @@ subroutine constitutive_restartWrite(fileHandle)
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
end subroutine constitutive_restartWrite end subroutine phase_restartWrite
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read data for restart !> @brief Read data for restart
! ToDo: Merge data into one file for MPI ! ToDo: Merge data into one file for MPI
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_restartRead(fileHandle) subroutine phase_restartRead(fileHandle)
integer(HID_T), intent(in) :: fileHandle integer(HID_T), intent(in) :: fileHandle
@ -691,7 +691,7 @@ subroutine constitutive_restartRead(fileHandle)
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph)) groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph))
call mech_restartRead(groupHandle(2),ph) call mechanical_restartRead(groupHandle(2),ph)
call HDF5_closeGroup(groupHandle(2)) call HDF5_closeGroup(groupHandle(2))
@ -699,7 +699,7 @@ subroutine constitutive_restartRead(fileHandle)
call HDF5_closeGroup(groupHandle(1)) call HDF5_closeGroup(groupHandle(1))
end subroutine constitutive_restartRead end subroutine phase_restartRead
end module phase end module phase

34
src/phase_damage.f90
View File

@ -196,7 +196,7 @@ end subroutine damage_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force !< @brief returns local part of nonlocal damage driving force
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point number ip, & !< integration point number
@ -246,7 +246,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
enddo enddo
enddo enddo
end subroutine constitutive_damage_getRateAndItsTangents end subroutine phase_damage_getRateAndItsTangents
@ -272,9 +272,9 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
size_so size_so
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(constitutive_source_maxSizeDotState) :: & real(pReal), dimension(phase_source_maxSizeDotState) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState real(pReal), dimension(phase_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
logical :: & logical :: &
converged_ converged_
@ -283,7 +283,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
converged_ = .true. converged_ = .true.
broken = constitutive_damage_collectDotState(co,ip,el,ph,me) broken = phase_damage_collectDotState(co,ip,el,ph,me)
if(broken) return if(broken) return
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -300,7 +300,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me)
enddo enddo
broken = constitutive_damage_collectDotState(co,ip,el,ph,me) broken = phase_damage_collectDotState(co,ip,el,ph,me)
if(broken) exit iteration if(broken) exit iteration
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -320,7 +320,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
enddo enddo
if(converged_) then if(converged_) then
broken = constitutive_damage_deltaState(mech_F_e(ph,me),ph,me) broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
exit iteration exit iteration
endif endif
@ -393,7 +393,7 @@ end subroutine damage_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken) function phase_damage_collectDotState(co,ip,el,ph,me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
co, & !< component-ID me integration point co, & !< component-ID me integration point
@ -419,7 +419,7 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken)
call anisoductile_dotState(co, ip, el) call anisoductile_dotState(co, ip, el)
case (DAMAGE_ANISOBRITTLE_ID) sourceType case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_dotState(mech_S(ph,me),co, ip, el) ! correct stress? call anisobrittle_dotState(mechanical_S(ph,me),co, ip, el) ! correct stress?
end select sourceType end select sourceType
@ -427,7 +427,7 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken)
enddo SourceLoop enddo SourceLoop
end function constitutive_damage_collectDotState end function phase_damage_collectDotState
@ -435,7 +435,7 @@ end function constitutive_damage_collectDotState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> will return false if delta state is not needed/supported by the constitutive model
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_damage_deltaState(Fe, ph, me) result(broken) function phase_damage_deltaState(Fe, ph, me) result(broken)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
@ -457,7 +457,7 @@ function constitutive_damage_deltaState(Fe, ph, me) result(broken)
sourceType: select case (phase_source(so,ph)) sourceType: select case (phase_source(so,ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType case (DAMAGE_ISOBRITTLE_ID) sourceType
call source_damage_isoBrittle_deltaState(constitutive_homogenizedC(ph,me), Fe, ph,me) call source_damage_isoBrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me))) broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,me)))
if(.not. broken) then if(.not. broken) then
myOffset = damageState(ph)%p(so)%offsetDeltaState myOffset = damageState(ph)%p(so)%offsetDeltaState
@ -470,7 +470,7 @@ function constitutive_damage_deltaState(Fe, ph, me) result(broken)
enddo SourceLoop enddo SourceLoop
end function constitutive_damage_deltaState end function phase_damage_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -507,7 +507,7 @@ end function source_active
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set damage parameter !< @brief Set damage parameter
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_damage_set_phi(phi,co,ce) module subroutine phase_damage_set_phi(phi,co,ce)
real(pReal), intent(in) :: phi real(pReal), intent(in) :: phi
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
@ -515,17 +515,17 @@ module subroutine constitutive_damage_set_phi(phi,co,ce)
current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi
end subroutine constitutive_damage_set_phi end subroutine phase_damage_set_phi
module function constitutive_damage_get_phi(co,ip,el) result(phi) module function phase_damage_get_phi(co,ip,el) result(phi)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal) :: phi real(pReal) :: phi
phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el)) phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
end function constitutive_damage_get_phi end function phase_damage_get_phi
end submodule damagee end submodule damagee

2
src/phase_damage_anisobrittle.f90
View File

@ -101,7 +101,7 @@ module function anisobrittle_init(source_length) result(mySources)
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

2
src/phase_damage_anisoductile.f90
View File

@ -87,7 +87,7 @@ module function anisoductile_init(source_length) result(mySources)
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

2
src/phase_damage_isobrittle.f90
View File

@ -74,7 +74,7 @@ module function isobrittle_init(source_length) result(mySources)
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nconstituents = count(material_phaseAt==p) * discretization_nIPs Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

2
src/phase_damage_isoductile.f90
View File

@ -78,7 +78,7 @@ module function isoductile_init(source_length) result(mySources)
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nconstituents=count(material_phaseAt==p) * discretization_nIPs Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) call phase_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0)
damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

384
src/phase_mechanics.f90 → src/phase_mechanical.f90
View File

@ -31,21 +31,21 @@ submodule(phase) mechanics
type(tTensorContainer), dimension(:), allocatable :: & type(tTensorContainer), dimension(:), allocatable :: &
! current value ! current value
constitutive_mech_Fe, & phase_mechanical_Fe, &
constitutive_mech_Fi, & phase_mechanical_Fi, &
constitutive_mech_Fp, & phase_mechanical_Fp, &
constitutive_mech_F, & phase_mechanical_F, &
constitutive_mech_Li, & phase_mechanical_Li, &
constitutive_mech_Lp, & phase_mechanical_Lp, &
constitutive_mech_S, & phase_mechanical_S, &
constitutive_mech_P, & phase_mechanical_P, &
! converged value at end of last solver increment ! converged value at end of last solver increment
constitutive_mech_Fi0, & phase_mechanical_Fi0, &
constitutive_mech_Fp0, & phase_mechanical_Fp0, &
constitutive_mech_F0, & phase_mechanical_F0, &
constitutive_mech_Li0, & phase_mechanical_Li0, &
constitutive_mech_Lp0, & phase_mechanical_Lp0, &
constitutive_mech_S0 phase_mechanical_S0
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: &
@ -97,7 +97,7 @@ submodule(phase) mechanics
broken broken
end function plastic_deltaState end function plastic_deltaState
module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el) S, Fi, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -114,7 +114,7 @@ submodule(phase) mechanics
dLi_dS, & !< derivative of Li with respect to S dLi_dS, & !< derivative of Li with respect to S
dLi_dFi dLi_dFi
end subroutine constitutive_LiAndItsTangents end subroutine phase_LiAndItsTangents
module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
@ -186,7 +186,7 @@ contains
!> @brief Initialize mechanical field related constitutive models !> @brief Initialize mechanical field related constitutive models
!> @details Initialize elasticity, plasticity and stiffness degradation models. !> @details Initialize elasticity, plasticity and stiffness degradation models.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_init(phases) module subroutine mechanical_init(phases)
class(tNode), pointer :: & class(tNode), pointer :: &
phases phases
@ -215,38 +215,38 @@ module subroutine mech_init(phases)
allocate(phase_NstiffnessDegradations(phases%length),source=0) allocate(phase_NstiffnessDegradations(phases%length),source=0)
allocate(output_constituent(phases%length)) allocate(output_constituent(phases%length))
allocate(constitutive_mech_Fe(phases%length)) allocate(phase_mechanical_Fe(phases%length))
allocate(constitutive_mech_Fi(phases%length)) allocate(phase_mechanical_Fi(phases%length))
allocate(constitutive_mech_Fi0(phases%length)) allocate(phase_mechanical_Fi0(phases%length))
allocate(constitutive_mech_Fp(phases%length)) allocate(phase_mechanical_Fp(phases%length))
allocate(constitutive_mech_Fp0(phases%length)) allocate(phase_mechanical_Fp0(phases%length))
allocate(constitutive_mech_F(phases%length)) allocate(phase_mechanical_F(phases%length))
allocate(constitutive_mech_F0(phases%length)) allocate(phase_mechanical_F0(phases%length))
allocate(constitutive_mech_Li(phases%length)) allocate(phase_mechanical_Li(phases%length))
allocate(constitutive_mech_Li0(phases%length)) allocate(phase_mechanical_Li0(phases%length))
allocate(constitutive_mech_Lp0(phases%length)) allocate(phase_mechanical_Lp0(phases%length))
allocate(constitutive_mech_Lp(phases%length)) allocate(phase_mechanical_Lp(phases%length))
allocate(constitutive_mech_S(phases%length)) allocate(phase_mechanical_S(phases%length))
allocate(constitutive_mech_P(phases%length)) allocate(phase_mechanical_P(phases%length))
allocate(constitutive_mech_S0(phases%length)) allocate(phase_mechanical_S0(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents))
allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents)) allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents))
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanics') mech => phase%get('mechanics')
@ -287,17 +287,17 @@ module subroutine mech_init(phases)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) & phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) &
/ math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3 phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3
constitutive_mech_F0(ph)%data(1:3,1:3,me) = math_I3 phase_mechanical_F0(ph)%data(1:3,1:3,me) = math_I3
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration phase_mechanical_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_F(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) phase_mechanical_F(ph)%data(1:3,1:3,me) = phase_mechanical_F0(ph)%data(1:3,1:3,me)
enddo enddo
enddo; enddo enddo; enddo
@ -307,14 +307,14 @@ module subroutine mech_init(phases)
! initialize plasticity ! initialize plasticity
allocate(plasticState(phases%length)) allocate(plasticState(phases%length))
allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID)
allocate(phase_plasticityInstance(phases%length),source = 0) allocate(phase_plasticInstance(phases%length),source = 0)
allocate(phase_localPlasticity(phases%length), source=.true.) allocate(phase_localPlasticity(phases%length), source=.true.)
call plastic_init() call plastic_init()
do ph = 1, phases%length do ph = 1, phases%length
phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph))
phase_plasticityInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph)) phase_plasticInstance(ph) = count(phase_plasticity(1:ph) == phase_plasticity(ph))
enddo enddo
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
@ -345,14 +345,14 @@ module subroutine mech_init(phases)
call eigendeformation_init(phases) call eigendeformation_init(phases)
end subroutine mech_init end subroutine mechanical_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law !> the elastic and intermediate deformation gradients using Hooke's law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, co, ip, el) Fe, Fi, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -376,7 +376,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
i, j, ph, me i, j, ph, me
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
C = math_66toSym3333(constitutive_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))) C = math_66toSym3333(phase_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el)) DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el))) degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el)))
@ -393,10 +393,10 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo; enddo enddo; enddo
end subroutine constitutive_hooke_SandItsTangents end subroutine phase_hooke_SandItsTangents
module subroutine mech_results(group,ph) module subroutine mechanical_results(group,ph)
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
integer, intent(in) :: ph integer, intent(in) :: ph
@ -407,28 +407,28 @@ module subroutine mech_results(group,ph)
select case(phase_plasticity(ph)) select case(phase_plasticity(ph))
case(PLASTICITY_ISOTROPIC_ID) case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_isotropic_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_PHENOPOWERLAW_ID) case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_phenopowerlaw_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_KINEHARDENING_ID) case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_kinehardening_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_DISLOTWIN_ID) case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_dislotwin_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_DISLOTUNGSTEN_ID) case(PLASTICITY_DISLOTUNGSTEN_ID)
call plastic_dislotungsten_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_dislotungsten_results(phase_plasticInstance(ph),group//'plastic/')
case(PLASTICITY_NONLOCAL_ID) case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(ph),group//'plastic/') call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/')
end select end select
call crystallite_results(group,ph) call crystallite_results(group,ph)
end subroutine mech_results end subroutine mechanical_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -503,8 +503,8 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
call plastic_dependentState(co,ip,el) call plastic_dependentState(co,ip,el)
Lpguess = constitutive_mech_Lp(ph)%data(1:3,1:3,me) ! take as first guess Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,me) ! take as first guess
Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me) ! take as first guess Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,me) ! take as first guess
call math_invert33(invFp_current,devNull,error,subFp0) call math_invert33(invFp_current,devNull,error,subFp0)
if (error) return ! error if (error) return ! error
@ -538,7 +538,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
B = math_I3 - Delta_t*Lpguess B = math_I3 - Delta_t*Lpguess
Fe = matmul(matmul(A,B), invFi_new) Fe = matmul(matmul(A,B), invFi_new)
call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, co, ip, el) Fe, Fi_new, co, ip, el)
call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
@ -582,7 +582,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
+ deltaLp * steplengthLp + deltaLp * steplengthLp
enddo LpLoop enddo LpLoop
call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, co, ip, el) S, Fi_new, co, ip, el)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
@ -633,13 +633,13 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
call math_invert33(Fp_new,devNull,error,invFp_new) call math_invert33(Fp_new,devNull,error,invFp_new)
if (error) return ! error if (error) return ! error
constitutive_mech_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) phase_mechanical_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
constitutive_mech_S(ph)%data(1:3,1:3,me) = S phase_mechanical_S(ph)%data(1:3,1:3,me) = S
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = Lpguess phase_mechanical_Lp(ph)%data(1:3,1:3,me) = Lpguess
constitutive_mech_Li(ph)%data(1:3,1:3,me) = Liguess phase_mechanical_Li(ph)%data(1:3,1:3,me) = Liguess
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize phase_mechanical_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = Fi_new phase_mechanical_Fi(ph)%data(1:3,1:3,me) = Fi_new
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new) phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new)
broken = .false. broken = .false.
end function integrateStress end function integrateStress
@ -668,9 +668,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
sizeDotState sizeDotState
real(pReal) :: & real(pReal) :: &
zeta zeta
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & real(pReal), dimension(phase_plasticity_maxSizeDotState) :: &
r ! state residuum r ! state residuum
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: &
dotState dotState
@ -796,7 +796,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
ph, & ph, &
me, & me, &
sizeDotState sizeDotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
@ -914,7 +914,7 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
ph, & ph, &
me, & me, &
sizeDotState sizeDotState
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
@ -987,28 +987,28 @@ subroutine crystallite_results(group,ph)
select case (output_constituent(ph)%label(ou)) select case (output_constituent(ph)%label(ou))
case('F') case('F')
call results_writeDataset(group//'/mechanics/',constitutive_mech_F(ph)%data,'F',& call results_writeDataset(group//'/mechanics/',phase_mechanical_F(ph)%data,'F',&
'deformation gradient','1') 'deformation gradient','1')
case('F_e') case('F_e')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Fe(ph)%data,'F_e',& call results_writeDataset(group//'/mechanics/',phase_mechanical_Fe(ph)%data,'F_e',&
'elastic deformation gradient','1') 'elastic deformation gradient','1')
case('F_p') case('F_p')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,'F_p', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Fp(ph)%data,'F_p', &
'plastic deformation gradient','1') 'plastic deformation gradient','1')
case('F_i') case('F_i')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,'F_i', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Fi(ph)%data,'F_i', &
'inelastic deformation gradient','1') 'inelastic deformation gradient','1')
case('L_p') case('L_p')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Lp(ph)%data,'L_p', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Lp(ph)%data,'L_p', &
'plastic velocity gradient','1/s') 'plastic velocity gradient','1/s')
case('L_i') case('L_i')
call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,'L_i', & call results_writeDataset(group//'/mechanics/',phase_mechanical_Li(ph)%data,'L_i', &
'inelastic velocity gradient','1/s') 'inelastic velocity gradient','1/s')
case('P') case('P')
call results_writeDataset(group//'/mechanics/',constitutive_mech_P(ph)%data,'P', & call results_writeDataset(group//'/mechanics/',phase_mechanical_P(ph)%data,'P', &
'First Piola-Kirchhoff stress','Pa') 'First Piola-Kirchhoff stress','Pa')
case('S') case('S')
call results_writeDataset(group//'/mechanics/',constitutive_mech_S(ph)%data,'S', & call results_writeDataset(group//'/mechanics/',phase_mechanical_S(ph)%data,'S', &
'Second Piola-Kirchhoff stress','Pa') 'Second Piola-Kirchhoff stress','Pa')
case('O') case('O')
select case(lattice_structure(ph)) select case(lattice_structure(ph))
@ -1067,43 +1067,43 @@ end subroutine crystallite_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward. !> @brief Wind homog inc forward.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_windForward(ph,me) module subroutine mechanical_windForward(ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
constitutive_mech_F0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) phase_mechanical_F0(ph)%data(1:3,1:3,me) = phase_mechanical_F(ph)%data(1:3,1:3,me)
constitutive_mech_Li0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) phase_mechanical_Li0(ph)%data(1:3,1:3,me) = phase_mechanical_Li(ph)%data(1:3,1:3,me)
constitutive_mech_Lp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) phase_mechanical_Lp0(ph)%data(1:3,1:3,me) = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
constitutive_mech_S0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) phase_mechanical_S0(ph)%data(1:3,1:3,me) = phase_mechanical_S(ph)%data(1:3,1:3,me)
plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me) plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me)
end subroutine mech_windForward end subroutine mechanical_windForward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment. !> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible? ! ToDo: Any guessing for the current states possible?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_forward() module subroutine mechanical_forward()
integer :: ph integer :: ph
do ph = 1, size(plasticState) do ph = 1, size(plasticState)
constitutive_mech_Fi0(ph) = constitutive_mech_Fi(ph) phase_mechanical_Fi0(ph) = phase_mechanical_Fi(ph)
constitutive_mech_Fp0(ph) = constitutive_mech_Fp(ph) phase_mechanical_Fp0(ph) = phase_mechanical_Fp(ph)
constitutive_mech_F0(ph) = constitutive_mech_F(ph) phase_mechanical_F0(ph) = phase_mechanical_F(ph)
constitutive_mech_Li0(ph) = constitutive_mech_Li(ph) phase_mechanical_Li0(ph) = phase_mechanical_Li(ph)
constitutive_mech_Lp0(ph) = constitutive_mech_Lp(ph) phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph)
constitutive_mech_S0(ph) = constitutive_mech_S(ph) phase_mechanical_S0(ph) = phase_mechanical_S(ph)
plasticState(ph)%state0 = plasticState(ph)%state plasticState(ph)%state0 = plasticState(ph)%state
enddo enddo
end subroutine mech_forward end subroutine mechanical_forward
@ -1111,19 +1111,19 @@ end subroutine mech_forward
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent !> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function constitutive_homogenizedC(ph,me) result(C) module function phase_homogenizedC(ph,me) result(C)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
plasticityType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
C = plastic_dislotwin_homogenizedC(ph,me) C = plastic_dislotwin_homogenizedC(ph,me)
case default plasticityType case default plasticType
C = lattice_C66(1:6,1:6,ph) C = lattice_C66(1:6,1:6,ph)
end select plasticityType end select plasticType
end function constitutive_homogenizedC end function phase_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1158,17 +1158,17 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
sizeDotState = plasticState(ph)%sizeDotState sizeDotState = plasticState(ph)%sizeDotState
subLi0 = constitutive_mech_Li0(ph)%data(1:3,1:3,me) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,me)
subLp0 = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%State0(:,me) subState0 = plasticState(ph)%State0(:,me)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me)
enddo enddo
subFp0 = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
subF0 = constitutive_mech_F0(ph)%data(1:3,1:3,me) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me)
subFrac = 0.0_pReal subFrac = 0.0_pReal
subStep = 1.0_pReal/num%subStepSizeCryst subStep = 1.0_pReal/num%subStepSizeCryst
todo = .true. todo = .true.
@ -1186,10 +1186,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
if (todo) then if (todo) then
subF0 = subF subF0 = subF
subLp0 = constitutive_mech_Lp(ph)%data(1:3,1:3,me) subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,me)
subFp0 = constitutive_mech_Fp(ph)%data(1:3,1:3,me) subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me)
subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me) subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me)
subState0 = plasticState(ph)%state(:,me) subState0 = plasticState(ph)%state(:,me)
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me)
@ -1199,12 +1199,12 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
! cut back (reduced time and restore) ! cut back (reduced time and restore)
else else
subStep = num%subStepSizeCryst * subStep subStep = num%subStepSizeCryst * subStep
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = subFp0 phase_mechanical_Fp(ph)%data(1:3,1:3,me) = subFp0
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = subFi0 phase_mechanical_Fi(ph)%data(1:3,1:3,me) = subFi0
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use? phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me) ! why no subS0 ? is S0 of any use?
if (subStep < 1.0_pReal) then ! actual (not initial) cutback if (subStep < 1.0_pReal) then ! actual (not initial) cutback
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = subLp0 phase_mechanical_Lp(ph)%data(1:3,1:3,me) = subLp0
constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0
endif endif
plasticState(ph)%state(:,me) = subState0 plasticState(ph)%state(:,me) = subState0
do so = 1, phase_Nsources(ph) do so = 1, phase_Nsources(ph)
@ -1218,9 +1218,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
! prepare for integration ! prepare for integration
if (todo) then if (todo) then
subF = subF0 & subF = subF0 &
+ subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_F0(ph)%data(1:3,1:3,me)) + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me))
constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) phase_mechanical_Fp(ph)%data(1:3,1:3,me))))
converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
endif endif
@ -1233,7 +1233,7 @@ end function crystallite_stress
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback. !> @brief Restore data after homog cutback.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine mech_restore(ce,includeL) module subroutine mechanical_restore(ce,includeL)
integer, intent(in) :: ce integer, intent(in) :: ce
logical, intent(in) :: & logical, intent(in) :: &
@ -1247,23 +1247,23 @@ module subroutine mech_restore(ce,includeL)
ph = material_phaseAt2(co,ce) ph = material_phaseAt2(co,ce)
me = material_phaseMemberAt2(co,ce) me = material_phaseMemberAt2(co,ce)
if (includeL) then if (includeL) then
constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me)
endif ! maybe protecting everything from overwriting makes more sense endif ! maybe protecting everything from overwriting makes more sense
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me)
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(ph)%data(1:3,1:3,me)
constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) phase_mechanical_S(ph)%data(1:3,1:3,me) = phase_mechanical_S0(ph)%data(1:3,1:3,me)
plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me) plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me)
enddo enddo
end subroutine mech_restore end subroutine mechanical_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF). !> @brief Calculate tangent (dPdF).
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
real(pReal), intent(in) :: dt real(pReal), intent(in) :: dt
integer, intent(in) :: & integer, intent(in) :: &
@ -1297,18 +1297,18 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
constitutive_mech_Fe(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
constitutive_mech_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me), & phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
co,ip,el) co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me))
invSubFp0 = math_inv33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me)) invSubFp0 = math_inv33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))
invSubFi0 = math_inv33(constitutive_mech_Fi0(ph)%data(1:3,1:3,me)) invSubFi0 = math_inv33(phase_mechanical_Fi0(ph)%data(1:3,1:3,me))
if (sum(abs(dLidS)) < tol_math_check) then if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal dFidS = 0.0_pReal
@ -1334,15 +1334,15 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
endif endif
call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
constitutive_mech_S(ph)%data(1:3,1:3,me), & phase_mechanical_S(ph)%data(1:3,1:3,me), &
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dSdF ! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi)) temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invSubFp0) temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invSubFp0)
temp_33_3 = matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp), invSubFi0) temp_33_3 = matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp), invSubFi0)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
@ -1370,9 +1370,9 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assemble dPdF ! assemble dPdF
temp_33_1 = matmul(constitutive_mech_S(ph)%data(1:3,1:3,me),transpose(invFp)) temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,me),transpose(invFp))
temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp) temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp)
temp_33_3 = matmul(temp_33_2,constitutive_mech_S(ph)%data(1:3,1:3,me)) temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,me))
dPdF = 0.0_pReal dPdF = 0.0_pReal
do p=1,3 do p=1,3
@ -1380,129 +1380,129 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
enddo enddo
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) & + matmul(matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) & + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo enddo; enddo
end function constitutive_mech_dPdF end function phase_mechanical_dPdF
module subroutine mech_restartWrite(groupHandle,ph) module subroutine mechanical_restartWrite(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
call HDF5_write(groupHandle,plasticState(ph)%state,'omega') call HDF5_write(groupHandle,plasticState(ph)%state,'omega')
call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,'F_i') call HDF5_write(groupHandle,phase_mechanical_Fi(ph)%data,'F_i')
call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,'L_i') call HDF5_write(groupHandle,phase_mechanical_Li(ph)%data,'L_i')
call HDF5_write(groupHandle,constitutive_mech_Lp(ph)%data,'L_p') call HDF5_write(groupHandle,phase_mechanical_Lp(ph)%data,'L_p')
call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,'F_p') call HDF5_write(groupHandle,phase_mechanical_Fp(ph)%data,'F_p')
call HDF5_write(groupHandle,constitutive_mech_S(ph)%data,'S') call HDF5_write(groupHandle,phase_mechanical_S(ph)%data,'S')
call HDF5_write(groupHandle,constitutive_mech_F(ph)%data,'F') call HDF5_write(groupHandle,phase_mechanical_F(ph)%data,'F')
end subroutine mech_restartWrite end subroutine mechanical_restartWrite
module subroutine mech_restartRead(groupHandle,ph) module subroutine mechanical_restartRead(groupHandle,ph)
integer(HID_T), intent(in) :: groupHandle integer(HID_T), intent(in) :: groupHandle
integer, intent(in) :: ph integer, intent(in) :: ph
call HDF5_read(groupHandle,plasticState(ph)%state0,'omega') call HDF5_read(groupHandle,plasticState(ph)%state0,'omega')
call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,'F_i') call HDF5_read(groupHandle,phase_mechanical_Fi0(ph)%data,'F_i')
call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,'L_i') call HDF5_read(groupHandle,phase_mechanical_Li0(ph)%data,'L_i')
call HDF5_read(groupHandle,constitutive_mech_Lp0(ph)%data,'L_p') call HDF5_read(groupHandle,phase_mechanical_Lp0(ph)%data,'L_p')
call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,'F_p') call HDF5_read(groupHandle,phase_mechanical_Fp0(ph)%data,'F_p')
call HDF5_read(groupHandle,constitutive_mech_S0(ph)%data,'S') call HDF5_read(groupHandle,phase_mechanical_S0(ph)%data,'S')
call HDF5_read(groupHandle,constitutive_mech_F0(ph)%data,'F') call HDF5_read(groupHandle,phase_mechanical_F0(ph)%data,'F')
end subroutine mech_restartRead end subroutine mechanical_restartRead
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) !< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mech_S(ph,me) result(S) module function mechanical_S(ph,me) result(S)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: S real(pReal), dimension(3,3) :: S
S = constitutive_mech_S(ph)%data(1:3,1:3,me) S = phase_mechanical_S(ph)%data(1:3,1:3,me)
end function mech_S end function mechanical_S
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics) !< @brief Get plastic velocity gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mech_L_p(ph,me) result(L_p) module function mechanical_L_p(ph,me) result(L_p)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: L_p real(pReal), dimension(3,3) :: L_p
L_p = constitutive_mech_Lp(ph)%data(1:3,1:3,me) L_p = phase_mechanical_Lp(ph)%data(1:3,1:3,me)
end function mech_L_p end function mechanical_L_p
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization) !< @brief Get deformation gradient (for use by homogenization)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function constitutive_mech_getF(co,ip,el) result(F) module function phase_mechanical_getF(co,ip,el) result(F)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: F real(pReal), dimension(3,3) :: F
F = constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) F = phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
end function constitutive_mech_getF end function phase_mechanical_getF
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics) !< @brief Get elastic deformation gradient (for use by non-mech physics)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function mech_F_e(ph,me) result(F_e) module function mechanical_F_e(ph,me) result(F_e)
integer, intent(in) :: ph,me integer, intent(in) :: ph,me
real(pReal), dimension(3,3) :: F_e real(pReal), dimension(3,3) :: F_e
F_e = constitutive_mech_Fe(ph)%data(1:3,1:3,me) F_e = phase_mechanical_Fe(ph)%data(1:3,1:3,me)
end function mech_F_e end function mechanical_F_e
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization) !< @brief Get second Piola-Kichhoff stress (for use by homogenization)
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module function constitutive_mech_getP(co,ip,el) result(P) module function phase_mechanical_getP(co,ip,el) result(P)
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
P = constitutive_mech_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) P = phase_mechanical_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
end function constitutive_mech_getP end function phase_mechanical_getP
! setter for homogenization ! setter for homogenization
module subroutine constitutive_mech_setF(F,co,ip,el) module subroutine phase_mechanical_setF(F,co,ip,el)
real(pReal), dimension(3,3), intent(in) :: F real(pReal), dimension(3,3), intent(in) :: F
integer, intent(in) :: co, ip, el integer, intent(in) :: co, ip, el
constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F
end subroutine constitutive_mech_setF end subroutine phase_mechanical_setF
end submodule mechanics end submodule mechanics

14
src/phase_mechanics_eigendeformation.f90 → src/phase_mechanical_eigen.f90
View File

@ -126,7 +126,7 @@ end function kinematics_active
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi? ! ToDo: MD: S is Mi?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, co, ip, el) S, Fi, co, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
@ -159,15 +159,15 @@ module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dS = 0.0_pReal dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_isotropic_ID) plasticityType case (PLASTICITY_isotropic_ID) plasticType
of = material_phasememberAt(co,ip,el) of = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(material_phaseAt(co,el)) instance = phase_plasticInstance(material_phaseAt(co,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType case default plasticType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal my_dLi_dS = 0.0_pReal
end select plasticityType end select plasticType
Li = Li + my_Li Li = Li + my_Li
dLi_dS = dLi_dS + my_dLi_dS dLi_dS = dLi_dS + my_dLi_dS
@ -201,7 +201,7 @@ module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
enddo; enddo enddo; enddo
end subroutine constitutive_LiAndItsTangents end subroutine phase_LiAndItsTangents
end submodule eigendeformation end submodule eigendeformation

0
src/phase_mechanics_eigendeformation_cleavageopening.f90 → src/phase_mechanical_eigen_cleavageopening.f90
View File

0
src/phase_mechanics_eigendeformation_slipplaneopening.f90 → src/phase_mechanical_eigen_slipplaneopening.f90
View File

0
src/phase_mechanics_eigendeformation_thermalexpansion.f90 → src/phase_mechanical_eigen_thermalexpansion.f90
View File

64
src/phase_mechanics_plastic.f90 → src/phase_mechanical_plastic.f90
View File

@ -270,31 +270,31 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_NONE_ID) plasticityType case (PLASTICITY_NONE_ID) plasticType
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticType
call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticType
call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el) call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me) call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp, thermal_T(ph,me),ph,me)
end select plasticityType end select plasticType
do i=1,3; do j=1,3 do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
@ -323,29 +323,29 @@ module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken)
logical :: broken logical :: broken
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me))
plasticityType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticType
call isotropic_dotState(Mp,ph,me) call isotropic_dotState(Mp,ph,me)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticType
call phenopowerlaw_dotState(Mp,ph,me) call phenopowerlaw_dotState(Mp,ph,me)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_dotState(Mp,ph,me) call plastic_kinehardening_dotState(Mp,ph,me)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me) call dislotwin_dotState(Mp,thermal_T(ph,me),ph,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me) call dislotungsten_dotState(Mp,thermal_T(ph,me),ph,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el) call nonlocal_dotState(Mp,thermal_T(ph,me),subdt,ph,me,ip,el)
end select plasticityType end select plasticType
broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,me)))
@ -369,20 +369,20 @@ module subroutine plastic_dependentState(co, ip, el)
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phasememberAt(co,ip,el) me = material_phasememberAt(co,ip,el)
instance = phase_plasticityInstance(ph) instance = phase_plasticInstance(ph)
plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticType
call dislotwin_dependentState(thermal_T(ph,me),instance,me) call dislotwin_dependentState(thermal_T(ph,me),instance,me)
case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dependentState(instance,me) call dislotungsten_dependentState(instance,me)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call nonlocal_dependentState(instance,me,ip,el) call nonlocal_dependentState(instance,me,ip,el)
end select plasticityType end select plasticType
end subroutine plastic_dependentState end subroutine plastic_dependentState
@ -410,24 +410,24 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken)
mySize mySize
Mp = matmul(matmul(transpose(constitutive_mech_Fi(ph)%data(1:3,1:3,me)),& Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,me)),&
constitutive_mech_Fi(ph)%data(1:3,1:3,me)),constitutive_mech_S(ph)%data(1:3,1:3,me)) phase_mechanical_Fi(ph)%data(1:3,1:3,me)),phase_mechanical_S(ph)%data(1:3,1:3,me))
instance = phase_plasticityInstance(ph) instance = phase_plasticInstance(ph)
plasticityType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticType
call plastic_kinehardening_deltaState(Mp,instance,me) call plastic_kinehardening_deltaState(Mp,instance,me)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticType
call plastic_nonlocal_deltaState(Mp,instance,me,ip,el) call plastic_nonlocal_deltaState(Mp,instance,me,ip,el)
broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me)))
case default case default
broken = .false. broken = .false.
end select plasticityType end select plasticType
if(.not. broken) then if(.not. broken) then
select case(phase_plasticity(ph)) select case(phase_plasticity(ph))

16
src/phase_mechanics_plastic_dislotungsten.f90 → src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -226,7 +226,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -289,16 +289,16 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
dot_gamma_pos,dot_gamma_neg, & dot_gamma_pos,dot_gamma_neg, &
ddot_gamma_dtau_pos,ddot_gamma_dtau_neg ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)))
call kinetics(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) call kinetics(Mp,T,phase_plasticInstance(ph),me,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -327,7 +327,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
real(pReal) :: & real(pReal) :: &
VacancyDiffusion VacancyDiffusion
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_pos, gdot_neg,& gdot_pos, gdot_neg,&
tau_pos,& tau_pos,&
tau_neg, & tau_neg, &
@ -336,10 +336,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me)
dot_rho_dip_climb, & dot_rho_dip_climb, &
dip_distance dip_distance
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),&
dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph)))
call kinetics(Mp,T,phase_plasticityInstance(ph),me,& call kinetics(Mp,T,phase_plasticInstance(ph),me,&
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg) tau_pos_out = tau_pos,tau_neg_out = tau_neg)

36
src/phase_mechanics_plastic_dislotwin.f90 → src/phase_mechanical_plastic_dislotwin.f90
View File

@ -415,7 +415,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
@ -496,8 +496,8 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC)
real(pReal) :: f_unrotated real(pReal) :: f_unrotated
associate(prm => param(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)),&
stt => state(phase_plasticityInstance(ph))) stt => state(phase_plasticInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
@ -535,11 +535,11 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
BoltzmannRatio, & BoltzmannRatio, &
ddot_gamma_dtau, & ddot_gamma_dtau, &
tau tau
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
dot_gamma_sl,ddot_gamma_dtau_slip dot_gamma_sl,ddot_gamma_dtau_slip
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
dot_gamma_twin,ddot_gamma_dtau_twin dot_gamma_twin,ddot_gamma_dtau_twin
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: &
dot_gamma_tr,ddot_gamma_dtau_trans dot_gamma_tr,ddot_gamma_dtau_trans
real(pReal):: dot_gamma_sb real(pReal):: dot_gamma_sb
real(pReal), dimension(3,3) :: eigVectors, P_sb real(pReal), dimension(3,3) :: eigVectors, P_sb
@ -564,7 +564,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
0, 1, 1 & 0, 1, 1 &
],pReal),[ 3,6]) ],pReal),[ 3,6])
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
@ -573,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip) call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl,ddot_gamma_dtau_slip)
slipContribution: do i = 1, prm%sum_N_sl slipContribution: do i = 1, prm%sum_N_sl
Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -581,7 +581,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
enddo slipContribution enddo slipContribution
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin) call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin,ddot_gamma_dtau_twin)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -589,7 +589,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
+ ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
enddo twinContibution enddo twinContibution
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans) call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr,ddot_gamma_dtau_trans)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -653,24 +653,24 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
b_d !< ratio of Burgers vector to stacking fault width b_d !< ratio of Burgers vector to stacking fault width
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
rho_dip_distance_min, & rho_dip_distance_min, &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
dot_gamma_twin dot_gamma_twin
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph)))
f_unrotated = 1.0_pReal & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,me)) & - sum(stt%f_tw(1:prm%sum_N_tw,me)) &
- sum(stt%f_tr(1:prm%sum_N_tr,me)) - sum(stt%f_tr(1:prm%sum_N_tr,me))
call kinetics_slip(Mp,T,phase_plasticityInstance(ph),me,dot_gamma_sl) call kinetics_slip(Mp,T,phase_plasticInstance(ph),me,dot_gamma_sl)
dot%gamma_sl(:,me) = abs(dot_gamma_sl) dot%gamma_sl(:,me) = abs(dot_gamma_sl)
rho_dip_distance_min = prm%D_a*prm%b_sl rho_dip_distance_min = prm%D_a*prm%b_sl
@ -721,10 +721,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) & - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
- dot_rho_dip_climb - dot_rho_dip_climb
call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_twin) call kinetics_twin(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_twin)
dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char
call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticityInstance(ph),me,dot_gamma_tr) call kinetics_trans(Mp,T,dot_gamma_sl,phase_plasticInstance(ph),me,dot_gamma_tr)
dot%f_tr(:,me) = f_unrotated*dot_gamma_tr dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
end associate end associate

8
src/phase_mechanics_plastic_isotropic.f90 → src/phase_mechanical_plastic_isotropic.f90
View File

@ -135,7 +135,7 @@ module function plastic_isotropic_init() result(myPlasticity)
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -190,7 +190,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: & integer :: &
k, l, m, n k, l, m, n
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)))
Mp_dev = math_deviatoric33(Mp) Mp_dev = math_deviatoric33(Mp)
squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev) squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
@ -275,8 +275,8 @@ module subroutine isotropic_dotState(Mp,ph,me)
xi_inf_star, & !< saturation xi xi_inf_star, & !< saturation xi
norm_Mp !< norm of the (deviatoric) Mandel stress norm_Mp !< norm of the (deviatoric) Mandel stress
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)))
if (prm%dilatation) then if (prm%dilatation) then
norm_Mp = sqrt(math_tensordot(Mp,Mp)) norm_Mp = sqrt(math_tensordot(Mp,Mp))

16
src/phase_mechanics_plastic_kinehardening.f90 → src/phase_mechanical_plastic_kinehardening.f90
View File

@ -180,7 +180,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -255,16 +255,16 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_pos,gdot_neg, & gdot_pos,gdot_neg, &
dgdot_dtau_pos,dgdot_dtau_neg dgdot_dtau_pos,dgdot_dtau_neg
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)))
call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
do i = 1, prm%sum_N_sl do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i) Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
@ -292,14 +292,14 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,me)
real(pReal) :: & real(pReal) :: &
sumGamma sumGamma
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_pos,gdot_neg gdot_pos,gdot_neg
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)),&
dot => dotState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)))
call kinetics(Mp,phase_plasticityInstance(ph),me,gdot_pos,gdot_neg) call kinetics(Mp,phase_plasticInstance(ph),me,gdot_pos,gdot_neg)
dot%accshear(:,me) = abs(gdot_pos+gdot_neg) dot%accshear(:,me) = abs(gdot_pos+gdot_neg)
sumGamma = sum(stt%accshear(:,me)) sumGamma = sum(stt%accshear(:,me))

2
src/phase_mechanics_plastic_none.f90 → src/phase_mechanical_plastic_none.f90
View File

@ -44,7 +44,7 @@ module function plastic_none_init() result(myPlasticity)
phase => phases%get(p) phase => phases%get(p)
if(.not. myPlasticity(p)) cycle if(.not. myPlasticity(p)) cycle
Nconstituents = count(material_phaseAt2 == p) Nconstituents = count(material_phaseAt2 == p)
call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0) call phase_allocateState(plasticState(p),Nconstituents,0,0,0)
enddo enddo
end function plastic_none_init end function plastic_none_init

64
src/phase_mechanics_plastic_nonlocal.f90 → src/phase_mechanical_plastic_nonlocal.f90
View File

@ -407,7 +407,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = pl%get_asBool('nonlocal') plasticState(p)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) & if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -642,8 +642,8 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
rho0 = getRho0(instance,me,ip,el) rho0 = getRho0(instance,me,ip,el)
if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
ph = material_phaseAt(1,el) ph = material_phaseAt(1,el)
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me))
invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me)) invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,me))
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg) rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg) rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@ -662,7 +662,7 @@ module subroutine nonlocal_dependentState(instance, me, ip, el)
neighbor_ip = IPneighborhood(2,n,ip,el) neighbor_ip = IPneighborhood(2,n,ip,el)
no = material_phasememberAt(1,neighbor_ip,neighbor_el) no = material_phasememberAt(1,neighbor_ip,neighbor_el)
if (neighbor_el > 0 .and. neighbor_ip > 0) then if (neighbor_el > 0 .and. neighbor_ip > 0) then
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el))
if (neighbor_instance == instance) then if (neighbor_instance == instance) then
nRealNeighbors = nRealNeighbors + 1.0_pReal nRealNeighbors = nRealNeighbors + 1.0_pReal
@ -782,25 +782,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
l, & l, &
t, & !< dislocation type t, & !< dislocation type
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: &
rhoSgl !< single dislocation densities (including blocked) rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: &
rho rho
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: &
v, & !< velocity v, & !< velocity
tauNS, & !< resolved shear stress including non Schmid and backstress terms tauNS, & !< resolved shear stress including non Schmid and backstress terms
dv_dtau, & !< velocity derivative with respect to the shear stress dv_dtau, & !< velocity derivative with respect to the shear stress
dv_dtauNS !< velocity derivative with respect to the shear stress dv_dtauNS !< velocity derivative with respect to the shear stress
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
tau, & !< resolved shear stress including backstress terms tau, & !< resolved shear stress including backstress terms
gdotTotal !< shear rate gdotTotal !< shear rate
associate(prm => param(phase_plasticityInstance(ph)),dst=>microstructure(phase_plasticityInstance(ph)),& associate(prm => param(phase_plasticInstance(ph)),dst=>microstructure(phase_plasticInstance(ph)),&
stt=>state(phase_plasticityInstance(ph))) stt=>state(phase_plasticInstance(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
!*** shortcut to state variables !*** shortcut to state variables
rho = getRho(phase_plasticityInstance(ph),me,ip,el) rho = getRho(phase_plasticInstance(ph),me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
do s = 1,ns do s = 1,ns
@ -820,7 +820,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
! edges ! edges
call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), & call kinetics(v(:,1), dv_dtau(:,1), dv_dtauNS(:,1), &
tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticityInstance(ph)) tau, tauNS(:,1), dst%tau_pass(:,me),1,Temperature, phase_plasticInstance(ph))
v(:,2) = v(:,1) v(:,2) = v(:,1)
dv_dtau(:,2) = dv_dtau(:,1) dv_dtau(:,2) = dv_dtau(:,1)
dv_dtauNS(:,2) = dv_dtauNS(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1)
@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
else else
do t = 3,4 do t = 3,4
call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), & call kinetics(v(:,t), dv_dtau(:,t), dv_dtauNS(:,t), &
tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticityInstance(ph)) tau, tauNS(:,t), dst%tau_pass(:,me),2,Temperature, phase_plasticInstance(ph))
enddo enddo
endif endif
@ -992,7 +992,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
c, & !< character of dislocation c, & !< character of dislocation
t, & !< type of dislocation t, & !< type of dislocation
s !< index of my current slip system s !< index of my current slip system
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: &
rho, & rho, &
rho0, & !< dislocation density at beginning of time step rho0, & !< dislocation density at beginning of time step
rhoDot, & !< density evolution rhoDot, & !< density evolution
@ -1000,17 +1000,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation
rhoDotThermalAnnihilation !< density evolution by thermal annihilation rhoDotThermalAnnihilation !< density evolution by thermal annihilation
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,8) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,8) :: &
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) my_rhoSgl0 !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,4) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,4) :: &
v, & !< current dislocation glide velocity v, & !< current dislocation glide velocity
v0, & v0, &
gdot !< shear rates gdot !< shear rates
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
tau, & !< current resolved shear stress tau, & !< current resolved shear stress
vClimb !< climb velocity of edge dipoles vClimb !< climb velocity of edge dipoles
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,2) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,2) :: &
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
dLower, & !< minimum stable dipole distance for edges and screws dLower, & !< minimum stable dipole distance for edges and screws
dUpper !< current maximum stable dipole distance for edges and screws dUpper !< current maximum stable dipole distance for edges and screws
@ -1022,22 +1022,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
return return
endif endif
associate(prm => param(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), &
dst => microstructure(phase_plasticityInstance(ph)), & dst => microstructure(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph)), & dot => dotState(phase_plasticInstance(ph)), &
stt => state(phase_plasticityInstance(ph))) stt => state(phase_plasticInstance(ph)))
ns = prm%sum_N_sl ns = prm%sum_N_sl
tau = 0.0_pReal tau = 0.0_pReal
gdot = 0.0_pReal gdot = 0.0_pReal
rho = getRho(phase_plasticityInstance(ph),me,ip,el) rho = getRho(phase_plasticInstance(ph),me,ip,el)
rhoSgl = rho(:,sgl) rhoSgl = rho(:,sgl)
rhoDip = rho(:,dip) rhoDip = rho(:,dip)
rho0 = getRho0(phase_plasticityInstance(ph),me,ip,el) rho0 = getRho0(phase_plasticInstance(ph),me,ip,el)
my_rhoSgl0 = rho0(:,sgl) my_rhoSgl0 = rho0(:,sgl)
forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticityInstance(ph)),me) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,phase_plasticInstance(ph)),me)
gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4) gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
#ifdef DEBUG #ifdef DEBUG
@ -1086,7 +1086,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
* sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4) * sqrt(stt%rho_forest(:,me)) / prm%i_sl / prm%b_sl, 2, 4)
endif isBCC endif isBCC
forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticityInstance(ph)),me) forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,phase_plasticInstance(ph)),me)
!**************************************************************************** !****************************************************************************
@ -1142,7 +1142,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
rhoDot = rhoDotFlux(timestep, phase_plasticityInstance(ph),me,ip,el) & rhoDot = rhoDotFlux(timestep, phase_plasticInstance(ph),me,ip,el) &
+ rhoDotMultiplication & + rhoDotMultiplication &
+ rhoDotSingle2DipoleGlide & + rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation & + rhoDotAthermalAnnihilation &
@ -1284,8 +1284,8 @@ function rhoDotFlux(timestep,instance,me,ip,el)
m(1:3,:,3) = -prm%slip_transverse m(1:3,:,3) = -prm%slip_transverse
m(1:3,:,4) = prm%slip_transverse m(1:3,:,4) = prm%slip_transverse
my_F = constitutive_mech_F(ph)%data(1:3,1:3,me) my_F = phase_mechanical_F(ph)%data(1:3,1:3,me)
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))) my_Fe = matmul(my_F, math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)))
neighbors: do n = 1,nIPneighbors neighbors: do n = 1,nIPneighbors
@ -1301,9 +1301,9 @@ function rhoDotFlux(timestep,instance,me,ip,el)
opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) neighbor_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el))
neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no) neighbor_F = phase_mechanical_F(np)%data(1:3,1:3,no)
neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no))) neighbor_Fe = matmul(neighbor_F, math_inv33(phase_mechanical_Fp(np)%data(1:3,1:3,no)))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F

22
src/phase_mechanics_plastic_phenopowerlaw.f90 → src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -231,7 +231,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
@ -300,18 +300,18 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
integer :: & integer :: &
i,k,l,m,n i,k,l,m,n
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
gdot_slip_pos,gdot_slip_neg, & gdot_slip_pos,gdot_slip_neg, &
dgdot_dtauslip_pos,dgdot_dtauslip_neg dgdot_dtauslip_pos,dgdot_dtauslip_neg
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tw) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tw) :: &
gdot_twin,dgdot_dtautwin gdot_twin,dgdot_dtautwin
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMp = 0.0_pReal dLp_dMp = 0.0_pReal
associate(prm => param(phase_plasticityInstance(ph))) associate(prm => param(phase_plasticInstance(ph)))
call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1, prm%sum_N_sl slipSystems: do i = 1, prm%sum_N_sl
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -320,7 +320,7 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me)
+ dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
enddo slipSystems enddo slipSystems
call kinetics_twin(Mp,phase_plasticityInstance(ph),me,gdot_twin,dgdot_dtautwin) call kinetics_twin(Mp,phase_plasticInstance(ph),me,gdot_twin,dgdot_dtautwin)
twinSystems: do i = 1, prm%sum_N_tw twinSystems: do i = 1, prm%sum_N_tw
Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -348,12 +348,12 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
c_SlipSlip,c_TwinSlip,c_TwinTwin, & c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset,& xi_slip_sat_offset,&
sumGamma,sumF sumGamma,sumF
real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl) :: & real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl) :: &
left_SlipSlip,right_SlipSlip, & left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg gdot_slip_pos,gdot_slip_neg
associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), &
dot => dotState(phase_plasticityInstance(ph))) dot => dotState(phase_plasticInstance(ph)))
sumGamma = sum(stt%gamma_slip(:,me)) sumGamma = sum(stt%gamma_slip(:,me))
sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char) sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char)
@ -373,9 +373,9 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! shear rates ! shear rates
call kinetics_slip(Mp,phase_plasticityInstance(ph),me,gdot_slip_pos,gdot_slip_neg) call kinetics_slip(Mp,phase_plasticInstance(ph),me,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg) dot%gamma_slip(:,me) = abs(gdot_slip_pos+gdot_slip_neg)
call kinetics_twin(Mp,phase_plasticityInstance(ph),me,dot%gamma_twin(:,me)) call kinetics_twin(Mp,phase_plasticInstance(ph),me,dot%gamma_twin(:,me))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hardening ! hardening

14
src/phase_thermal.f90
View File

@ -124,7 +124,7 @@ end subroutine thermal_init
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate !< @brief calculates thermal dissipation rate
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRate(TDot, ph,me) module subroutine phase_thermal_getRate(TDot, ph,me)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
@ -153,13 +153,13 @@ module subroutine constitutive_thermal_getRate(TDot, ph,me)
enddo enddo
end subroutine constitutive_thermal_getRate end subroutine phase_thermal_getRate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken) function phase_thermal_collectDotState(ph,me) result(broken)
integer, intent(in) :: ph, me integer, intent(in) :: ph, me
logical :: broken logical :: broken
@ -178,7 +178,7 @@ function constitutive_thermal_collectDotState(ph,me) result(broken)
enddo SourceLoop enddo SourceLoop
end function constitutive_thermal_collectDotState end function phase_thermal_collectDotState
module function thermal_stress(Delta_t,ph,me) result(converged_) module function thermal_stress(Delta_t,ph,me) result(converged_)
@ -207,7 +207,7 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
so, & so, &
sizeDotState sizeDotState
broken = constitutive_thermal_collectDotState(ph,me) broken = phase_thermal_collectDotState(ph,me)
if(broken) return if(broken) return
do so = 1, thermal_Nsources(ph) do so = 1, thermal_Nsources(ph)
@ -264,7 +264,7 @@ end function thermal_dot_T
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
!< @brief Set temperature !< @brief Set temperature
!---------------------------------------------------------------------------------------------- !----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_setField(T,dot_T, co,ce) module subroutine phase_thermal_setField(T,dot_T, co,ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
@ -273,7 +273,7 @@ module subroutine constitutive_thermal_setField(T,dot_T, co,ce)
current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T
current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T
end subroutine constitutive_thermal_setField end subroutine phase_thermal_setField

4
src/phase_thermal_dissipation.f90
View File

@ -56,7 +56,7 @@ module function dissipation_init(source_length) result(mySources)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
Nconstituents = count(material_phaseAt2 == ph) Nconstituents = count(material_phaseAt2 == ph)
call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0)
end associate end associate
endif endif
@ -78,7 +78,7 @@ module subroutine dissipation_getRate(TDot, ph,me)
associate(prm => param(ph)) associate(prm => param(ph))
TDot = prm%kappa*sum(abs(mech_S(ph,me)*mech_L_p(ph,me))) TDot = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me)))
end associate end associate
end subroutine dissipation_getRate end subroutine dissipation_getRate

2
src/phase_thermal_externalheat.f90
View File

@ -69,7 +69,7 @@ module function externalheat_init(source_length) result(mySources)
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
Nconstituents = count(material_phaseAt2 == ph) Nconstituents = count(material_phaseAt2 == ph)
call constitutive_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0)
end associate end associate
endif endif
enddo enddo