diff --git a/PRIVATE b/PRIVATE index 5d27d879a..e74cf0062 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 5d27d879abaeb1542e3f9f3065172be740af2899 +Subproject commit e74cf00628285a587ced1e551cc9837c1011ca1c diff --git a/examples/SpectralMethod/Polycrystal/shearXY.yaml b/examples/SpectralMethod/Polycrystal/shearXY.yaml index fe746e7ac..58471e8f1 100644 --- a/examples/SpectralMethod/Polycrystal/shearXY.yaml +++ b/examples/SpectralMethod/Polycrystal/shearXY.yaml @@ -3,9 +3,6 @@ solver: mechanical: spectral_basic -initial_conditions: - T: 300 #in Kelvin - loadstep: - boundary_conditions: mechanical: diff --git a/examples/SpectralMethod/Polycrystal/shearZX.yaml b/examples/SpectralMethod/Polycrystal/shearZX.yaml index 5d977b4ba..a32fafb85 100644 --- a/examples/SpectralMethod/Polycrystal/shearZX.yaml +++ b/examples/SpectralMethod/Polycrystal/shearZX.yaml @@ -3,9 +3,6 @@ solver: mechanical: spectral_basic -initial_conditions: - T: 300 #in Kelvin - loadstep: - boundary_conditions: mechanical: diff --git a/examples/SpectralMethod/Polycrystal/tensionX.yaml b/examples/SpectralMethod/Polycrystal/tensionX.yaml index 0809fd53d..870755d58 100644 --- a/examples/SpectralMethod/Polycrystal/tensionX.yaml +++ b/examples/SpectralMethod/Polycrystal/tensionX.yaml @@ -3,9 +3,6 @@ solver: mechanical: spectral_basic -initial_conditions: - T: 300 #in Kelvin - loadstep: - boundary_conditions: mechanical: diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index d3ed972e3..0adc494b8 100644 --- a/python/damask/_configmaterial.py +++ b/python/damask/_configmaterial.py @@ -168,11 +168,11 @@ class ConfigMaterial(Config): ok = False if 'material' in self: - for i,v in enumerate(self['material']): - if 'constituents' in v: - f = 0.0 - for c in v['constituents']: - f+= float(c['v']) + for i,m in enumerate(self['material']): + if 'constituents' in m: + v = 0.0 + for c in m['constituents']: + v+= float(c['v']) if 'O' in c: try: Rotation.from_quaternion(c['O']) @@ -180,8 +180,8 @@ class ConfigMaterial(Config): o = c['O'] print(f"Invalid orientation: '{o}' in material '{i}'") ok = False - if not np.isclose(f,1.0): - print(f"Invalid total fraction '{f}' in material '{i}'") + if not np.isclose(v,1.0): + print(f"Invalid total fraction (v) '{v}' in material '{i}'") ok = False return ok diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index ad017ee1e..00a07efa5 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -347,23 +347,21 @@ class VTK: 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: - import wx - _ = wx.App(False) # noqa - width, height = wx.GetDisplaySize() - except ImportError: - try: - import tkinter - tk = tkinter.Tk() - width = tk.winfo_screenwidth() - height = tk.winfo_screenheight() - tk.destroy() - except Exception as e: - width = 1024 - height = 768 + import tkinter + tk = tkinter.Tk() + width = tk.winfo_screenwidth() + height = tk.winfo_screenheight() + tk.destroy() + except Exception as e: + width = 1024 + height = 768 - return (width,height) mapper = vtk.vtkDataSetMapper() mapper.SetInputData(self.vtk_data) actor = vtk.vtkActor() @@ -377,7 +375,7 @@ class VTK: ren.AddActor(actor) ren.SetBackground(0.2,0.2,0.2) - window.SetSize(screen_size[0],screen_size[1]) + window.SetSize(width,height) iren = vtk.vtkRenderWindowInteractor() iren.SetRenderWindow(window) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index e441339cb..f0e510433 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -85,7 +85,7 @@ subroutine CPFEM_initAll call discretization_marc_init call lattice_init call material_init(.false.) - call constitutive_init + call phase_init call homogenization_init call crystallite_init call CPFEM_init @@ -257,7 +257,7 @@ end subroutine CPFEM_general subroutine CPFEM_forward call homogenization_forward - call constitutive_forward + call phase_forward end subroutine CPFEM_forward @@ -272,7 +272,7 @@ subroutine CPFEM_results(inc,time) call results_openJobFile call results_addIncrement(inc,time) - call constitutive_results + call phase_results call homogenization_results call discretization_results call results_finalizeIncrement diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 2b32a0cbb..e57de7f67 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -60,7 +60,7 @@ subroutine CPFEM_initAll call discretization_grid_init(restart=interface_restartInc>0) #endif call material_init(restart=interface_restartInc>0) - call constitutive_init + call phase_init call homogenization_init call crystallite_init call CPFEM_init @@ -87,7 +87,7 @@ subroutine CPFEM_init fileHandle = HDF5_openFile(fileName) call homogenization_restartRead(fileHandle) - call constitutive_restartRead(fileHandle) + call phase_restartRead(fileHandle) call HDF5_closeFile(fileHandle) endif @@ -110,7 +110,7 @@ subroutine CPFEM_restartWrite fileHandle = HDF5_openFile(fileName,'a') call homogenization_restartWrite(fileHandle) - call constitutive_restartWrite(fileHandle) + call phase_restartWrite(fileHandle) call HDF5_closeFile(fileHandle) @@ -123,7 +123,7 @@ end subroutine CPFEM_restartWrite subroutine CPFEM_forward call homogenization_forward - call constitutive_forward + call phase_forward end subroutine CPFEM_forward @@ -138,7 +138,7 @@ subroutine CPFEM_results(inc,time) call results_openJobFile call results_addIncrement(inc,time) - call constitutive_results + call phase_results call homogenization_results call discretization_results call results_finalizeIncrement diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 6480f7382..a191298a3 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -20,19 +20,19 @@ #include "material.f90" #include "lattice.f90" #include "phase.f90" -#include "phase_mechanics.f90" -#include "phase_mechanics_plastic.f90" -#include "phase_mechanics_plastic_none.f90" -#include "phase_mechanics_plastic_isotropic.f90" -#include "phase_mechanics_plastic_phenopowerlaw.f90" -#include "phase_mechanics_plastic_kinehardening.f90" -#include "phase_mechanics_plastic_dislotwin.f90" -#include "phase_mechanics_plastic_dislotungsten.f90" -#include "phase_mechanics_plastic_nonlocal.f90" -#include "phase_mechanics_eigendeformation.f90" -#include "phase_mechanics_eigendeformation_cleavageopening.f90" -#include "phase_mechanics_eigendeformation_slipplaneopening.f90" -#include "phase_mechanics_eigendeformation_thermalexpansion.f90" +#include "phase_mechanical.f90" +#include "phase_mechanical_plastic.f90" +#include "phase_mechanical_plastic_none.f90" +#include "phase_mechanical_plastic_isotropic.f90" +#include "phase_mechanical_plastic_phenopowerlaw.f90" +#include "phase_mechanical_plastic_kinehardening.f90" +#include "phase_mechanical_plastic_dislotwin.f90" +#include "phase_mechanical_plastic_dislotungsten.f90" +#include "phase_mechanical_plastic_nonlocal.f90" +#include "phase_mechanical_eigen.f90" +#include "phase_mechanical_eigen_cleavageopening.f90" +#include "phase_mechanical_eigen_slipplaneopening.f90" +#include "phase_mechanical_eigen_thermalexpansion.f90" #include "phase_thermal.f90" #include "phase_thermal_dissipation.f90" #include "phase_thermal_externalheat.f90" @@ -44,10 +44,10 @@ #include "damage_none.f90" #include "damage_nonlocal.f90" #include "homogenization.f90" -#include "homogenization_mechanics.f90" -#include "homogenization_mechanics_none.f90" -#include "homogenization_mechanics_isostrain.f90" -#include "homogenization_mechanics_RGC.f90" +#include "homogenization_mechanical.f90" +#include "homogenization_mechanical_pass.f90" +#include "homogenization_mechanical_isostrain.f90" +#include "homogenization_mechanical_RGC.f90" #include "homogenization_thermal.f90" #include "homogenization_damage.f90" #include "CPFEM.f90" diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 7e649ff38..ed1ada171 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -18,9 +18,9 @@ program DAMASK_grid use CPFEM2 use material use spectral_utilities - use grid_mech_spectral_basic - use grid_mech_spectral_polarisation - use grid_mech_FEM + use grid_mechanical_spectral_basic + use grid_mechanical_spectral_polarisation + use grid_mechanical_FEM use grid_damage_spectral use grid_thermal_spectral use results @@ -83,16 +83,16 @@ program DAMASK_grid type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tSolutionState), allocatable, dimension(:) :: solres - procedure(grid_mech_spectral_basic_init), pointer :: & - mech_init - procedure(grid_mech_spectral_basic_forward), pointer :: & - mech_forward - procedure(grid_mech_spectral_basic_solution), pointer :: & - mech_solution - procedure(grid_mech_spectral_basic_updateCoords), pointer :: & - mech_updateCoords - procedure(grid_mech_spectral_basic_restartWrite), pointer :: & - mech_restartWrite + procedure(grid_mechanical_spectral_basic_init), pointer :: & + mechanical_init + procedure(grid_mechanical_spectral_basic_forward), pointer :: & + mechanical_forward + procedure(grid_mechanical_spectral_basic_solution), pointer :: & + mechanical_solution + procedure(grid_mechanical_spectral_basic_updateCoords), pointer :: & + mechanical_updateCoords + procedure(grid_mechanical_spectral_basic_restartWrite), pointer :: & + mechanical_restartWrite external :: & quit @@ -139,25 +139,25 @@ program DAMASK_grid debug_grid => config_debug%get('grid',defaultVal=emptyList) select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) case ('Basic') - mech_init => grid_mech_spectral_basic_init - mech_forward => grid_mech_spectral_basic_forward - mech_solution => grid_mech_spectral_basic_solution - mech_updateCoords => grid_mech_spectral_basic_updateCoords - mech_restartWrite => grid_mech_spectral_basic_restartWrite + mechanical_init => grid_mechanical_spectral_basic_init + mechanical_forward => grid_mechanical_spectral_basic_forward + mechanical_solution => grid_mechanical_spectral_basic_solution + mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords + mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite case ('Polarisation') - mech_init => grid_mech_spectral_polarisation_init - mech_forward => grid_mech_spectral_polarisation_forward - mech_solution => grid_mech_spectral_polarisation_solution - mech_updateCoords => grid_mech_spectral_polarisation_updateCoords - mech_restartWrite => grid_mech_spectral_polarisation_restartWrite + mechanical_init => grid_mechanical_spectral_polarisation_init + mechanical_forward => grid_mechanical_spectral_polarisation_forward + mechanical_solution => grid_mechanical_spectral_polarisation_solution + mechanical_updateCoords => grid_mechanical_spectral_polarisation_updateCoords + mechanical_restartWrite => grid_mechanical_spectral_polarisation_restartWrite case ('FEM') - mech_init => grid_mech_FEM_init - mech_forward => grid_mech_FEM_forward - mech_solution => grid_mech_FEM_solution - mech_updateCoords => grid_mech_FEM_updateCoords - mech_restartWrite => grid_mech_FEM_restartWrite + mechanical_init => grid_mechanical_FEM_init + mechanical_forward => grid_mechanical_FEM_forward + mechanical_solution => grid_mechanical_FEM_solution + mechanical_updateCoords => grid_mechanical_FEM_updateCoords + mechanical_restartWrite => grid_mechanical_FEM_restartWrite case default 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 select case (loadCases(1)%ID(field)) case(FIELD_MECH_ID) - call mech_init + call mechanical_init case(FIELD_THERMAL_ID) call grid_thermal_spectral_init @@ -379,7 +379,7 @@ program DAMASK_grid do field = 1, nActiveFields select case(loadCases(l)%ID(field)) case(FIELD_MECH_ID) - call mech_forward (& + call mechanical_forward (& cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & deformation_BC = loadCases(l)%deformation, & stress_BC = loadCases(l)%stress, & @@ -399,7 +399,7 @@ program DAMASK_grid do field = 1, nActiveFields select case(loadCases(l)%ID(field)) case(FIELD_MECH_ID) - solres(field) = mech_solution(incInfo) + solres(field) = mechanical_solution(incInfo) case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(timeinc) case(FIELD_DAMAGE_ID) @@ -420,7 +420,7 @@ program DAMASK_grid if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged .and. .not. solres(1)%termIll) then ! and acceptable solution found - call mech_updateCoords + call mechanical_updateCoords timeIncOld = timeinc cutBack = .false. 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) if (ierr /= 0) error stop 'MPI error' if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then - call mech_restartWrite + call mechanical_restartWrite call CPFEM_restartWrite endif if(signal) call interface_setSIGUSR2(.false.) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 003f568c6..806973a4d 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -4,7 +4,7 @@ !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @brief Grid solver for mechanics: FEM !-------------------------------------------------------------------------------------------------- -module grid_mech_FEM +module grid_mechanical_FEM #include #include use PETScdmda @@ -45,8 +45,8 @@ module grid_mech_FEM !-------------------------------------------------------------------------------------------------- ! PETSc data - DM :: mech_grid - SNES :: mech_snes + DM :: mechanical_grid + SNES :: mechanical_snes Vec :: solution_current, solution_lastInc, solution_rate !-------------------------------------------------------------------------------------------------- @@ -79,18 +79,18 @@ module grid_mech_FEM totalIter = 0 !< total iteration in current increment public :: & - grid_mech_FEM_init, & - grid_mech_FEM_solution, & - grid_mech_FEM_forward, & - grid_mech_FEM_updateCoords, & - grid_mech_FEM_restartWrite + grid_mechanical_FEM_init, & + grid_mechanical_FEM_solution, & + grid_mechanical_FEM_forward, & + grid_mechanical_FEM_updateCoords, & + grid_mechanical_FEM_restartWrite contains !-------------------------------------------------------------------------------------------------- !> @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, dimension(4,8) :: & @@ -114,7 +114,7 @@ subroutine grid_mech_FEM_init num_grid, & debug_grid - print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT) + print'(/,a)', ' <<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT) !------------------------------------------------------------------------------------------------- ! debugging options @@ -141,8 +141,11 @@ subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres & - &-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS, & + '-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) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) @@ -155,8 +158,10 @@ subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr) - call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) + call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr) + CHKERRQ(ierr) + call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr) + CHKERRQ(ierr) localK = 0 localK(worldrank) = grid3 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, & 3, 1, & [grid(1)],[grid(2)],localK, & - mech_grid,ierr) + mechanical_grid,ierr) CHKERRQ(ierr) - call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(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) + call SNESSetDM(mechanical_snes,mechanical_grid,ierr) CHKERRQ(ierr) - call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(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) + call DMsetFromOptions(mechanical_grid,ierr) CHKERRQ(ierr) - call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr) + call DMsetUp(mechanical_grid,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) - 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 call VecSet(solution_current,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 DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,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) xend = xstart+xend-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 F, & ! target F 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) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr) + call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr) CHKERRQ(ierr) restartRead2: if (interface_restartInc > 0) then @@ -257,13 +272,13 @@ subroutine grid_mech_FEM_init endif restartRead2 -end subroutine grid_mech_FEM_init +end subroutine grid_mechanical_FEM_init !-------------------------------------------------------------------------------------------------- !> @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 @@ -284,11 +299,13 @@ function grid_mech_FEM_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! 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 - call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) + call SNESGetConvergedReason(mechanical_snes,reason,ierr) + CHKERRQ(ierr) solution%converged = reason > 0 solution%iterationsNeeded = totalIter @@ -296,14 +313,14 @@ function grid_mech_FEM_solution(incInfoIn) result(solution) terminallyIll = .false. 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 !> @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) 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 - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) if (cutBack) then 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 DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,ierr) + CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr) + CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! 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%timeinc = Delta_t -end subroutine grid_mech_FEM_forward +end subroutine grid_mechanical_FEM_forward !-------------------------------------------------------------------------------------------------- !> @brief Update coordinates !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_FEM_updateCoords +subroutine grid_mechanical_FEM_updateCoords 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 !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_FEM_restartWrite +subroutine grid_mechanical_FEM_restartWrite PetscErrorCode :: ierr integer(HID_T) :: fileHandle, groupHandle PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc character(len=pStringLen) :: fileName - call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) - call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) + call DMDAVecGetArrayF90(mechanical_grid,solution_current,u_current,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) @@ -427,10 +450,12 @@ subroutine grid_mech_FEM_restartWrite call HDF5_closeGroup(groupHandle) call HDF5_closeFile(fileHandle) - call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) - call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr) + call DMDAVecRestoreArrayF90(mechanical_grid,solution_current,u_current,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 real(pReal), dimension(3,3,3,3) :: devNull - call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr) - call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr) + call SNESGetNumberFunctionEvals(mechanical_snes,nfuncs,ierr) + CHKERRQ(ierr) + call SNESGetIterationNumber(mechanical_snes,PETScIter,ierr) + CHKERRQ(ierr) 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 module grid_mech_FEM +end module grid_mechanical_FEM diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 9bc36165f..f3f30c0af 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Grid solver for mechanics: Spectral basic !-------------------------------------------------------------------------------------------------- -module grid_mech_spectral_basic +module grid_mechanical_spectral_basic #include #include use PETScdmda @@ -79,18 +79,18 @@ module grid_mech_spectral_basic totalIter = 0 !< total iteration in current increment public :: & - grid_mech_spectral_basic_init, & - grid_mech_spectral_basic_solution, & - grid_mech_spectral_basic_forward, & - grid_mech_spectral_basic_updateCoords, & - grid_mech_spectral_basic_restartWrite + grid_mechanical_spectral_basic_init, & + grid_mechanical_spectral_basic_solution, & + grid_mechanical_spectral_basic_forward, & + grid_mechanical_spectral_basic_updateCoords, & + grid_mechanical_spectral_basic_restartWrite contains !-------------------------------------------------------------------------------------------------- !> @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 PetscErrorCode :: ierr @@ -105,7 +105,7 @@ subroutine grid_mech_spectral_basic_init num_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*, '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 - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',ierr) CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) @@ -152,7 +152,7 @@ subroutine grid_mech_spectral_basic_init !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc 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(worldrank) = grid3 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_saveReferenceStiffness -end subroutine grid_mech_spectral_basic_init +end subroutine grid_mechanical_spectral_basic_init !-------------------------------------------------------------------------------------------------- !> @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 @@ -262,14 +262,14 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution) terminallyIll = .false. 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 !> @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) 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%timeinc = Delta_t -end subroutine grid_mech_spectral_basic_forward +end subroutine grid_mechanical_spectral_basic_forward !-------------------------------------------------------------------------------------------------- !> @brief Update coordinates !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_spectral_basic_updateCoords +subroutine grid_mechanical_spectral_basic_updateCoords PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F @@ -354,13 +354,13 @@ subroutine grid_mech_spectral_basic_updateCoords call utilities_updateCoords(F) 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 !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_spectral_basic_restartWrite +subroutine grid_mechanical_spectral_basic_restartWrite PetscErrorCode :: ierr integer(HID_T) :: fileHandle, groupHandle @@ -393,7 +393,7 @@ subroutine grid_mech_spectral_basic_restartWrite 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 module grid_mech_spectral_basic +end module grid_mechanical_spectral_basic diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 7160c1adc..633fced7f 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Grid solver for mechanics: Spectral Polarisation !-------------------------------------------------------------------------------------------------- -module grid_mech_spectral_polarisation +module grid_mechanical_spectral_polarisation #include #include use PETScdmda @@ -90,18 +90,18 @@ module grid_mech_spectral_polarisation totalIter = 0 !< total iteration in current increment public :: & - grid_mech_spectral_polarisation_init, & - grid_mech_spectral_polarisation_solution, & - grid_mech_spectral_polarisation_forward, & - grid_mech_spectral_polarisation_updateCoords, & - grid_mech_spectral_polarisation_restartWrite + grid_mechanical_spectral_polarisation_init, & + grid_mechanical_spectral_polarisation_solution, & + grid_mechanical_spectral_polarisation_forward, & + grid_mechanical_spectral_polarisation_updateCoords, & + grid_mechanical_spectral_polarisation_restartWrite contains !-------------------------------------------------------------------------------------------------- !> @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 PetscErrorCode :: ierr @@ -118,7 +118,7 @@ subroutine grid_mech_spectral_polarisation_init num_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*, '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 - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type ngmres',ierr) CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) @@ -172,7 +172,7 @@ subroutine grid_mech_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc 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(worldrank) = grid3 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 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 !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) +function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution @@ -294,14 +294,14 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) terminallyIll = .false. 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 !> @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) 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%timeinc = Delta_t -end subroutine grid_mech_spectral_polarisation_forward +end subroutine grid_mechanical_spectral_polarisation_forward !-------------------------------------------------------------------------------------------------- !> @brief Update coordinates !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_spectral_polarisation_updateCoords +subroutine grid_mechanical_spectral_polarisation_updateCoords PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau @@ -408,13 +408,13 @@ subroutine grid_mech_spectral_polarisation_updateCoords call utilities_updateCoords(FandF_tau(0:8,:,:,:)) 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 !-------------------------------------------------------------------------------------------------- -subroutine grid_mech_spectral_polarisation_restartWrite +subroutine grid_mechanical_spectral_polarisation_restartWrite PetscErrorCode :: ierr 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) -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 module grid_mech_spectral_polarisation +end module grid_mechanical_spectral_polarisation diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 2255e83b0..62ac52f06 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -45,10 +45,10 @@ module homogenization !-------------------------------------------------------------------------------------------------- interface - module subroutine mech_init(num_homog) + module subroutine mechanical_init(num_homog) class(tNode), pointer, intent(in) :: & num_homog !< pointer to mechanical homogenization numerics data - end subroutine mech_init + end subroutine mechanical_init module subroutine thermal_init end subroutine thermal_init @@ -56,13 +56,13 @@ module homogenization module 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) :: & subF integer, intent(in) :: & ip, & !< integration point el !< element number - end subroutine mech_partition + end subroutine mechanical_partition module subroutine thermal_partition(ce) integer, intent(in) :: ce @@ -76,19 +76,19 @@ module homogenization integer, intent(in) :: ip,el end subroutine thermal_homogenize - module subroutine mech_homogenize(dt,ip,el) + module subroutine mechanical_homogenize(dt,ip,el) real(pReal), intent(in) :: dt integer, intent(in) :: & ip, & !< integration point 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 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) :: & subdt !< current time step real(pReal), intent(in), dimension(3,3) :: & @@ -97,7 +97,7 @@ module homogenization ip, & !< integration point el !< element number logical, dimension(2) :: doneAndHappy - end function mech_updateState + end function mechanical_updateState 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') - call mech_init(num_homog) + call mechanical_init(num_homog) call thermal_init() call damage_init() @@ -253,7 +253,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ce = (el-1)*discretization_nIPs + ip 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(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 - 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. do co = 1, myNgrains 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 doneAndHappy = [.true.,.false.] 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) endif endif @@ -338,7 +338,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE do co = 1, homogenization_Nconstituents(ho) call crystallite_orientations(co,ip,el) enddo - call mech_homogenize(dt,ip,el) + call mechanical_homogenize(dt,ip,el) enddo IpLooping3 enddo elementLooping3 !$OMP END DO @@ -365,7 +365,7 @@ subroutine homogenization_results group_base = 'current/homogenization/'//trim(material_name_homogenization(ho)) call results_closeGroup(results_addGroup(group_base)) - call mech_results(group_base,ho) + call mechanical_results(group_base,ho) group = trim(group_base)//'/damage' call results_closeGroup(results_addGroup(group)) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index ba021a0e0..8b73ed54a 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -74,7 +74,7 @@ module subroutine damage_partition(ce) phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(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 end subroutine damage_partition @@ -120,7 +120,7 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p phiDot = 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) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) diff --git a/src/homogenization_mechanics.f90 b/src/homogenization_mechanical.f90 similarity index 79% rename from src/homogenization_mechanics.f90 rename to src/homogenization_mechanical.f90 index 4eae859ee..91e1d3d84 100644 --- a/src/homogenization_mechanics.f90 +++ b/src/homogenization_mechanical.f90 @@ -7,52 +7,52 @@ submodule(homogenization) mechanics interface - module subroutine mech_none_init - end subroutine mech_none_init + module subroutine mechanical_pass_init + end subroutine mechanical_pass_init - module subroutine mech_isostrain_init - end subroutine mech_isostrain_init + module subroutine mechanical_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) :: & 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 (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 (3,3), intent(in) :: avgF !< average deformation gradient at material point integer, intent(in) :: & instance, & 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,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses 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,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses 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 real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -63,13 +63,13 @@ submodule(homogenization) mechanics integer, intent(in) :: & ip, & !< integration point 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 character(len=*), intent(in) :: group !< group name in HDF5 file - end subroutine mech_RGC_results + end subroutine mechanical_RGC_results end interface @@ -78,7 +78,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Allocate variables and set parameters. !-------------------------------------------------------------------------------------------------- -module subroutine mech_init(num_homog) +module subroutine mechanical_init(num_homog) class(tNode), pointer, intent(in) :: & 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) num_homogMech => num_homog%get('mech',defaultVal=emptyDict) - if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init - if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech) + if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mechanical_pass_init + if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mechanical_isostrain_init + 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. !-------------------------------------------------------------------------------------------------- -module subroutine mech_partition(subF,ip,el) +module subroutine mechanical_partition(subF,ip,el) real(pReal), intent(in), dimension(3,3) :: & subF @@ -122,25 +122,25 @@ module subroutine mech_partition(subF,ip,el) Fs(1:3,1:3,1) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mech_isostrain_partitionDeformation(Fs,subF) + call mechanical_isostrain_partitionDeformation(Fs,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mech_RGC_partitionDeformation(Fs,subF,ip,el) + call mechanical_RGC_partitionDeformation(Fs,subF,ip,el) end select chosenHomogenization 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 -end subroutine mech_partition +end subroutine mechanical_partition !-------------------------------------------------------------------------------------------------- !> @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 integer, intent(in) :: & @@ -156,15 +156,15 @@ module subroutine mech_homogenize(dt,ip,el) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el) - homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,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) = phase_mechanical_dPdF(dt,1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) - Ps(:,:,co) = constitutive_mech_getP(co,ip,el) + dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el) + Ps(:,:,co) = phase_mechanical_getP(co,ip,el) enddo - call mech_isostrain_averageStressAndItsTangent(& + call mechanical_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & @@ -172,10 +172,10 @@ module subroutine mech_homogenize(dt,ip,el) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) - Ps(:,:,co) = constitutive_mech_getP(co,ip,el) + dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el) + Ps(:,:,co) = phase_mechanical_getP(co,ip,el) enddo - call mech_RGC_averageStressAndItsTangent(& + call mechanical_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & @@ -183,14 +183,14 @@ module subroutine mech_homogenize(dt,ip,el) 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 !> "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) :: & 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 do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el) - Fs(:,:,co) = constitutive_mech_getF(co,ip,el) - Ps(:,:,co) = constitutive_mech_getP(co,ip,el) + dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el) + Fs(:,:,co) = phase_mechanical_getF(co,ip,el) + Ps(:,:,co) = phase_mechanical_getP(co,ip,el) enddo - doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) + doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) else doneAndHappy = .true. endif -end function mech_updateState +end function mechanical_updateState !-------------------------------------------------------------------------------------------------- !> @brief Write results to file. !-------------------------------------------------------------------------------------------------- -module subroutine mech_results(group_base,h) +module subroutine mechanical_results(group_base,h) use material, only: & material_homogenization_type => homogenization_type @@ -239,7 +239,7 @@ module subroutine mech_results(group_base,h) select case(material_homogenization_type(h)) case(HOMOGENIZATION_rgc_ID) - call mech_RGC_results(homogenization_typeInstance(h),group) + call mechanical_RGC_results(homogenization_typeInstance(h),group) end select @@ -250,7 +250,7 @@ module subroutine mech_results(group_base,h) !call results_writeDataset(group,temp,'P',& ! '1st Piola-Kirchhoff stress','Pa') -end subroutine mech_results +end subroutine mechanical_results end submodule mechanics diff --git a/src/homogenization_mechanics_RGC.f90 b/src/homogenization_mechanical_RGC.f90 similarity index 98% rename from src/homogenization_mechanics_RGC.f90 rename to src/homogenization_mechanical_RGC.f90 index 89a84314b..546c2b65c 100644 --- a/src/homogenization_mechanics_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -71,7 +71,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @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) :: & 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) 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%c_alpha = homogMech%get_asFloat('c_alpha') @@ -190,13 +190,13 @@ module subroutine mech_RGC_init(num_homogMech) enddo -end subroutine mech_RGC_init +end subroutine mechanical_RGC_init !-------------------------------------------------------------------------------------------------- !> @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 @@ -229,14 +229,14 @@ module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of) 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 ! "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 real(pReal), dimension(:,:,:), intent(in) :: & 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 - 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') 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 function mech_RGC_updateState +end function mechanical_RGC_updateState !-------------------------------------------------------------------------------------------------- !> @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,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) 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 !-------------------------------------------------------------------------------------------------- -module subroutine mech_RGC_results(instance,group) +module subroutine mechanical_RGC_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -754,7 +754,7 @@ module subroutine mech_RGC_results(instance,group) enddo outputsLoop end associate -end subroutine mech_RGC_results +end subroutine mechanical_RGC_results !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization_mechanics_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 similarity index 91% rename from src/homogenization_mechanics_isostrain.f90 rename to src/homogenization_mechanical_isostrain.f90 index 9634b2a38..f9616de42 100644 --- a/src/homogenization_mechanics_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -26,7 +26,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mech_isostrain_init +module subroutine mechanical_isostrain_init integer :: & Ninstances, & @@ -58,7 +58,7 @@ module subroutine mech_isostrain_init case ('avg') prm%mapping = average_ID case default - call IO_error(211,ext_msg='sum'//' (mech_isostrain)') + call IO_error(211,ext_msg='sum'//' (mechanical_isostrain)') end select Nmaterialpoints = count(material_homogenizationAt == h) @@ -70,13 +70,13 @@ module subroutine mech_isostrain_init enddo -end subroutine mech_isostrain_init +end subroutine mechanical_isostrain_init !-------------------------------------------------------------------------------------------------- !> @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 @@ -84,13 +84,13 @@ module subroutine mech_isostrain_partitionDeformation(F,avgF) 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 !-------------------------------------------------------------------------------------------------- -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,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 subroutine mech_isostrain_averageStressAndItsTangent +end subroutine mechanical_isostrain_averageStressAndItsTangent end submodule isostrain diff --git a/src/homogenization_mechanics_none.f90 b/src/homogenization_mechanical_pass.f90 similarity index 87% rename from src/homogenization_mechanics_none.f90 rename to src/homogenization_mechanical_pass.f90 index ebe2ea8f1..9e8f3e44c 100644 --- a/src/homogenization_mechanics_none.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -11,14 +11,14 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mech_none_init +module subroutine mechanical_pass_init integer :: & Ninstances, & h, & Nmaterialpoints - print'(/,a)', ' <<<+- homogenization:mechanics:none init -+>>>' + print'(/,a)', ' <<<+- homogenization:mechanics:pass init -+>>>' Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID) 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_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) homogState(h)%sizeState = 0 @@ -36,6 +36,6 @@ module subroutine mech_none_init enddo -end subroutine mech_none_init +end subroutine mechanical_pass_init end submodule none diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index b7170e6ec..67027f52a 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -78,7 +78,7 @@ module subroutine thermal_partition(ce) T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(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 end subroutine thermal_partition @@ -91,7 +91,7 @@ module subroutine thermal_homogenize(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 @@ -235,7 +235,7 @@ module subroutine thermal_conduction_getSource(Tdot, ip,el) do co = 1, homogenization_Nconstituents(ho) ph = material_phaseAt(co,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 enddo diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 7369520c1..5ef0f7a36 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -18,7 +18,7 @@ program DAMASK_mesh use config use discretization_mesh use FEM_Utilities - use mesh_mech_FEM + use mesh_mechanical_FEM implicit none @@ -242,7 +242,7 @@ program DAMASK_mesh do field = 1, nActiveFields select case (loadCases(1)%fieldBC(field)%ID) case(FIELD_MECH_ID) - call FEM_mech_init(loadCases(1)%fieldBC(field)) + call FEM_mechanical_init(loadCases(1)%fieldBC(field)) end select enddo @@ -306,7 +306,7 @@ program DAMASK_mesh do field = 1, nActiveFields select case (loadCases(currentLoadCase)%fieldBC(field)%ID) case(FIELD_MECH_ID) - call FEM_mech_forward (& + call FEM_mechanical_forward (& guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) end select @@ -320,7 +320,7 @@ program DAMASK_mesh do field = 1, nActiveFields select case (loadCases(currentLoadCase)%fieldBC(field)%ID) case(FIELD_MECH_ID) - solres(field) = FEM_mech_solution (& + solres(field) = FEM_mechanical_solution (& incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field)) end select diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 2f3633e11..4b3be8a42 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -127,12 +127,12 @@ subroutine FEM_utilities_init CHKERRQ(ierr) if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr) CHKERRQ(ierr) - call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls & - &-mech_snes_linesearch_type cp -mech_snes_ksp_ew & - &-mech_snes_ksp_ew_rtol0 0.01 -mech_snes_ksp_ew_rtolmax 0.01 & - &-mech_ksp_type fgmres -mech_ksp_max_it 25 & - &-mech_pc_type ml -mech_mg_levels_ksp_type chebyshev & - &-mech_mg_levels_pc_type sor -mech_pc_ml_nullspace user',ierr) + call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls & + &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew & + &-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 & + &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25 & + &-mechanical_pc_type ml -mechanical_mg_levels_ksp_type chebyshev & + &-mechanical_mg_levels_pc_type sor -mechanical_pc_ml_nullspace user',ierr) CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('petsc_options',defaultVal=''),ierr) CHKERRQ(ierr) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index e19c35998..c811d842b 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief FEM PETSc solver !-------------------------------------------------------------------------------------------------- -module mesh_mech_FEM +module mesh_mechanical_FEM #include #include #include @@ -50,7 +50,7 @@ module mesh_mech_FEM type(tNumerics), private :: num !-------------------------------------------------------------------------------------------------- ! PETSc data - SNES :: mech_snes + SNES :: mechanical_snes Vec :: solution, solution_rate, solution_local PetscInt :: dimPlex, cellDof, nQuadrature, nBasis PetscReal, allocatable, target :: qPoints(:), qWeights(:) @@ -65,20 +65,20 @@ module mesh_mech_FEM real(pReal), parameter :: eps = 1.0e-18_pReal public :: & - FEM_mech_init, & - FEM_mech_solution, & - FEM_mech_forward + FEM_mechanical_init, & + FEM_mechanical_solution, & + FEM_mechanical_forward contains !-------------------------------------------------------------------------------------------------- !> @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 - DM :: mech_mesh + DM :: mechanical_mesh PetscFE :: mechFE PetscQuadrature :: mechQuad, functional PetscDS :: mechDS @@ -126,8 +126,8 @@ subroutine FEM_mech_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech mesh - call DMClone(geomMesh,mech_mesh,ierr); CHKERRQ(ierr) - call DMGetDimension(mech_mesh,dimPlex,ierr); CHKERRQ(ierr) + call DMClone(geomMesh,mechanical_mesh,ierr); CHKERRQ(ierr) + call DMGetDimension(mechanical_mesh,dimPlex,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization @@ -146,22 +146,22 @@ subroutine FEM_mech_init(fieldBC) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc - call DMAddField(mech_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr) - call DMCreateDS(mech_mesh,ierr); CHKERRQ(ierr) - call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) + call DMAddField(mechanical_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr) + call DMCreateDS(mechanical_mesh,ierr); CHKERRQ(ierr) + call DMGetDS(mechanical_mesh,mechDS,ierr); CHKERRQ(ierr) call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr) call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech boundary conditions - call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) - call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr) - call DMGetLocalSection(mech_mesh,section,ierr); CHKERRQ(ierr) + call DMGetLabel(mechanical_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr) + call DMPlexLabelComplete(mechanical_mesh,BCLabel,ierr); CHKERRQ(ierr) + call DMGetLocalSection(mechanical_mesh,section,ierr); CHKERRQ(ierr) allocate(pnumComp(1), source=dimPlex) allocate(pnumDof(0:dimPlex), source = 0) do topologDim = 0, dimPlex - call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr) + call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,ierr) CHKERRQ(ierr) call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr) CHKERRQ(ierr) @@ -179,10 +179,10 @@ subroutine FEM_mech_init(fieldBC) numBC = numBC + 1 call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),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) 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) call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(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 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) CHKERRQ(ierr) - call DMSetSection(mech_mesh,section,ierr); CHKERRQ(ierr) + call DMSetSection(mechanical_mesh,section,ierr); CHKERRQ(ierr) do faceSet = 1, numBC call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr) enddo !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc - call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr) - call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) - call SNESSetDM(mech_snes,mech_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(mech_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 DMSNESSetFunctionLocal(mech_mesh,FEM_mech_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces + call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr);CHKERRQ(ierr) + call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr);CHKERRQ(ierr) + call SNESSetDM(mechanical_snes,mechanical_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver + call DMCreateGlobalVector(mechanical_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs + call DMCreateGlobalVector(mechanical_mesh,solution_rate ,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(mechanical_mesh,FEM_mechanical_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces 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) - call SNESSetMaxLinearSolveFailures(mech_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 SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures + call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,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) - call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) + call SNESSetFromOptions(mechanical_snes,ierr); CHKERRQ(ierr) !-------------------------------------------------------------------------------------------------- ! init fields @@ -236,11 +236,11 @@ subroutine FEM_mech_init(fieldBC) call PetscDSGetDiscretization(mechDS,0,mechFE,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) do cell = cellStart, cellEnd-1 !< loop over all elements 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) cellJMat = reshape(pCellJ,shape=[dimPlex,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) enddo 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) enddo -end subroutine FEM_mech_init +end subroutine FEM_mechanical_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the FEM load step !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function FEM_mech_solution( & +type(tSolutionState) function FEM_mechanical_solution( & incInfoIn,timeinc,timeinc_old,fieldBC) !-------------------------------------------------------------------------------------------------- @@ -278,35 +278,35 @@ type(tSolutionState) function FEM_mech_solution( & SNESConvergedReason :: reason incInfo = incInfoIn - FEM_mech_solution%converged =.false. + FEM_mechanical_solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc 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 SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? + call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) + call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? terminallyIll = .false. if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error - FEM_mech_solution%converged = .false. - FEM_mech_solution%iterationsNeeded = num%itmax + FEM_mechanical_solution%converged = .false. + FEM_mechanical_solution%iterationsNeeded = num%itmax else ! >= 1 proper convergence (or terminally ill) - FEM_mech_solution%converged = .true. - call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr) + FEM_mechanical_solution%converged = .true. + call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,ierr) CHKERRQ(ierr) endif print'(/,a)', ' ===========================================================================' flush(IO_STDOUT) -end function FEM_mech_solution +end function FEM_mechanical_solution !-------------------------------------------------------------------------------------------------- !> @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 PetscObject,intent(in) :: dummy @@ -431,13 +431,13 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr) enddo 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 !-------------------------------------------------------------------------------------------------- -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 @@ -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 MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr) -end subroutine FEM_mech_formJacobian +end subroutine FEM_mechanical_formJacobian !-------------------------------------------------------------------------------------------------- !> @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) :: & fieldBC @@ -603,7 +603,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC) if (guess .and. .not. cutBack) then ForwardData = .True. 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 DMGetLocalVector(dm_local,x_local,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 VecScale(solution,timeinc,ierr); CHKERRQ(ierr) -end subroutine FEM_mech_forward +end subroutine FEM_mechanical_forward !-------------------------------------------------------------------------------------------------- !> @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 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 flush(IO_STDOUT) -end subroutine FEM_mech_converged +end subroutine FEM_mechanical_converged -end module mesh_mech_FEM +end module mesh_mechanical_FEM diff --git a/src/phase.f90 b/src/phase.f90 index dd14e0831..628bfc443 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -62,7 +62,7 @@ module phase phase_Nsources, & !< number of source 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_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 logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) @@ -75,15 +75,15 @@ module phase integer, public, protected :: & - constitutive_plasticity_maxSizeDotState, & - constitutive_source_maxSizeDotState + phase_plasticity_maxSizeDotState, & + phase_source_maxSizeDotState interface ! == cleaned:begin ================================================================================= - module subroutine mech_init(phases) + module subroutine mechanical_init(phases) class(tNode), pointer :: phases - end subroutine mech_init + end subroutine mechanical_init module subroutine damage_init end subroutine damage_init @@ -93,83 +93,83 @@ module phase end subroutine thermal_init - module subroutine mech_results(group,ph) + module subroutine mechanical_results(group,ph) character(len=*), intent(in) :: group integer, intent(in) :: ph - end subroutine mech_results + end subroutine mechanical_results module subroutine damage_results(group,ph) character(len=*), intent(in) :: group integer, intent(in) :: ph end subroutine damage_results - module subroutine mech_windForward(ph,me) + module subroutine mechanical_windForward(ph,me) integer, intent(in) :: ph, me - end subroutine mech_windForward + end subroutine mechanical_windForward - module subroutine mech_forward() - end subroutine mech_forward + module subroutine mechanical_forward() + end subroutine mechanical_forward module subroutine thermal_forward() end subroutine thermal_forward - module subroutine mech_restore(ce,includeL) + module subroutine mechanical_restore(ce,includeL) integer, intent(in) :: ce 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 integer, intent(in) :: & co, & !< counter in constituent loop ip, & !< counter in integration point loop el !< counter in element loop 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, 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, 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 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 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 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 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 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 real(pReal) :: phi - end function constitutive_damage_get_phi + end function phase_damage_get_phi module function thermal_T(ph,me) result(T) integer, intent(in) :: ph,me @@ -182,20 +182,20 @@ module phase 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 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 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 integer, intent(in) :: co, ce - end subroutine constitutive_damage_set_phi + end subroutine phase_damage_set_phi ! == cleaned:end =================================================================================== @@ -222,13 +222,13 @@ module phase logical :: converged_ 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 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) :: & ip, & !< integration point number el !< element number @@ -237,13 +237,13 @@ module phase real(pReal), intent(inout) :: & phiDot, & 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 real(pReal), intent(out) :: & TDot - end subroutine constitutive_thermal_getRate + end subroutine phase_thermal_getRate module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) integer, intent(in) :: & @@ -281,39 +281,39 @@ module phase #endif public :: & - constitutive_init, & - constitutive_homogenizedC, & - constitutive_damage_getRateAndItsTangents, & - constitutive_thermal_getRate, & - constitutive_results, & - constitutive_allocateState, & - constitutive_forward, & - constitutive_restore, & + phase_init, & + phase_homogenizedC, & + phase_damage_getRateAndItsTangents, & + phase_thermal_getRate, & + phase_results, & + phase_allocateState, & + phase_forward, & + phase_restore, & plastic_nonlocal_updateCompatibility, & converged, & crystallite_init, & crystallite_stress, & thermal_stress, & - constitutive_mech_dPdF, & + phase_mechanical_dPdF, & crystallite_orientations, & crystallite_push33ToRef, & - constitutive_restartWrite, & - constitutive_restartRead, & + phase_restartWrite, & + phase_restartRead, & integrateDamageState, & - constitutive_thermal_setField, & - constitutive_damage_set_phi, & - constitutive_damage_get_phi, & - constitutive_mech_getP, & - constitutive_mech_setF, & - constitutive_mech_getF, & - constitutive_windForward + phase_thermal_setField, & + phase_damage_set_phi, & + phase_damage_get_phi, & + phase_mechanical_getP, & + phase_mechanical_setF, & + phase_mechanical_getF, & + phase_windForward contains !-------------------------------------------------------------------------------------------------- !> @brief Initialze constitutive models for individual physics !-------------------------------------------------------------------------------------------------- -subroutine constitutive_init +subroutine phase_init integer :: & ph, & !< counter in phase loop @@ -336,12 +336,12 @@ subroutine constitutive_init phases => config_material%get('phase') - call mech_init(phases) + call mechanical_init(phases) call damage_init call thermal_init(phases) - constitutive_source_maxSizeDotState = 0 + phase_source_maxSizeDotState = 0 PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state @@ -350,18 +350,18 @@ subroutine constitutive_init damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 end forall - constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & + phase_source_maxSizeDotState = max(phase_source_maxSizeDotState, & maxval(damageState(ph)%p%sizeDotState)) 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 !-------------------------------------------------------------------------------------------------- -subroutine constitutive_allocateState(state, & +subroutine phase_allocateState(state, & Nconstituents,sizeState,sizeDotState,sizeDeltaState) class(tState), intent(out) :: & @@ -387,13 +387,13 @@ subroutine constitutive_allocateState(state, & allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) -end subroutine constitutive_allocateState +end subroutine phase_allocateState !-------------------------------------------------------------------------------------------------- !> @brief Restore data after homog cutback. !-------------------------------------------------------------------------------------------------- -subroutine constitutive_restore(ce,includeL) +subroutine phase_restore(ce,includeL) logical, intent(in) :: includeL integer, intent(in) :: ce @@ -410,21 +410,21 @@ subroutine constitutive_restore(ce,includeL) 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. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -subroutine constitutive_forward() +subroutine phase_forward() integer :: ph, so - call mech_forward() + call mechanical_forward() call thermal_forward() do ph = 1, size(damageState) @@ -432,13 +432,13 @@ subroutine constitutive_forward() damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state enddo; enddo -end subroutine constitutive_forward +end subroutine phase_forward !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine constitutive_results() +subroutine phase_results() integer :: ph character(len=:), allocatable :: group @@ -451,12 +451,12 @@ subroutine constitutive_results() group = '/current/phase/'//trim(material_name_phase(ph))//'/' call results_closeGroup(results_addGroup(group)) - call mech_results(group,ph) + call mechanical_results(group,ph) call damage_results(group,ph) enddo -end subroutine constitutive_results +end subroutine phase_results !-------------------------------------------------------------------------------------------------- @@ -557,7 +557,7 @@ end subroutine crystallite_init !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- -subroutine constitutive_windForward(ip,el) +subroutine phase_windForward(ip,el) integer, intent(in) :: & ip, & !< integration point number @@ -572,7 +572,7 @@ subroutine constitutive_windForward(ip,el) ph = material_phaseAt(co,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)) damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) @@ -580,7 +580,7 @@ subroutine constitutive_windForward(ip,el) 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(& - 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) & 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 @@ -620,7 +620,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33) 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)) @@ -648,7 +648,7 @@ end function converged !> @brief Write current restart information (Field and constitutive data) to file. ! ToDo: Merge data into one file for MPI !-------------------------------------------------------------------------------------------------- -subroutine constitutive_restartWrite(fileHandle) +subroutine phase_restartWrite(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)) - call mech_restartWrite(groupHandle(2),ph) + call mechanical_restartWrite(groupHandle(2),ph) call HDF5_closeGroup(groupHandle(2)) @@ -670,14 +670,14 @@ subroutine constitutive_restartWrite(fileHandle) call HDF5_closeGroup(groupHandle(1)) -end subroutine constitutive_restartWrite +end subroutine phase_restartWrite !-------------------------------------------------------------------------------------------------- !> @brief Read data for restart ! ToDo: Merge data into one file for MPI !-------------------------------------------------------------------------------------------------- -subroutine constitutive_restartRead(fileHandle) +subroutine phase_restartRead(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)) - call mech_restartRead(groupHandle(2),ph) + call mechanical_restartRead(groupHandle(2),ph) call HDF5_closeGroup(groupHandle(2)) @@ -699,7 +699,7 @@ subroutine constitutive_restartRead(fileHandle) call HDF5_closeGroup(groupHandle(1)) -end subroutine constitutive_restartRead +end subroutine phase_restartRead end module phase diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 014347c1a..57145c550 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -196,7 +196,7 @@ end subroutine damage_init !---------------------------------------------------------------------------------------------- !< @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) :: & ip, & !< integration point number @@ -246,7 +246,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi 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 real(pReal) :: & zeta - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + real(pReal), dimension(phase_source_maxSizeDotState) :: & 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 :: & converged_ @@ -283,7 +283,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) me = material_phaseMemberAt(co,ip,el) converged_ = .true. - broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) return 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) enddo - broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + broken = phase_damage_collectDotState(co,ip,el,ph,me) if(broken) exit iteration do so = 1, phase_Nsources(ph) @@ -320,7 +320,7 @@ module function integrateDamageState(dt,co,ip,el) result(broken) enddo 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 endif @@ -393,7 +393,7 @@ end subroutine damage_results !-------------------------------------------------------------------------------------------------- !> @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) :: & 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) 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 @@ -427,7 +427,7 @@ function constitutive_damage_collectDotState(co,ip,el,ph,me) result(broken) 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 !> 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) :: & ph, & @@ -457,7 +457,7 @@ function constitutive_damage_deltaState(Fe, ph, me) result(broken) sourceType: select case (phase_source(so,ph)) 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))) if(.not. broken) then myOffset = damageState(ph)%p(so)%offsetDeltaState @@ -470,7 +470,7 @@ function constitutive_damage_deltaState(Fe, ph, me) result(broken) enddo SourceLoop -end function constitutive_damage_deltaState +end function phase_damage_deltaState !-------------------------------------------------------------------------------------------------- @@ -507,7 +507,7 @@ end function source_active !---------------------------------------------------------------------------------------------- !< @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 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 -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 real(pReal) :: phi 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 diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 2ae5ca951..51ca7786a 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -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' 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) if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index 8d5904e6b..e54eff201 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -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' 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) if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 091377171..529ecd404 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -74,7 +74,7 @@ module function isobrittle_init(source_length) result(mySources) if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' 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) if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index 6f0a2d0fb..71f99078b 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -78,7 +78,7 @@ module function isoductile_init(source_length) result(mySources) if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' 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) if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' diff --git a/src/phase_mechanics.f90 b/src/phase_mechanical.f90 similarity index 82% rename from src/phase_mechanics.f90 rename to src/phase_mechanical.f90 index a84c1a385..efd4d51d2 100644 --- a/src/phase_mechanics.f90 +++ b/src/phase_mechanical.f90 @@ -31,21 +31,21 @@ submodule(phase) mechanics type(tTensorContainer), dimension(:), allocatable :: & ! current value - constitutive_mech_Fe, & - constitutive_mech_Fi, & - constitutive_mech_Fp, & - constitutive_mech_F, & - constitutive_mech_Li, & - constitutive_mech_Lp, & - constitutive_mech_S, & - constitutive_mech_P, & + phase_mechanical_Fe, & + phase_mechanical_Fi, & + phase_mechanical_Fp, & + phase_mechanical_F, & + phase_mechanical_Li, & + phase_mechanical_Lp, & + phase_mechanical_S, & + phase_mechanical_P, & ! converged value at end of last solver increment - constitutive_mech_Fi0, & - constitutive_mech_Fp0, & - constitutive_mech_F0, & - constitutive_mech_Li0, & - constitutive_mech_Lp0, & - constitutive_mech_S0 + phase_mechanical_Fi0, & + phase_mechanical_Fp0, & + phase_mechanical_F0, & + phase_mechanical_Li0, & + phase_mechanical_Lp0, & + phase_mechanical_S0 integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & @@ -97,7 +97,7 @@ submodule(phase) mechanics broken 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) integer, intent(in) :: & @@ -114,7 +114,7 @@ submodule(phase) mechanics dLi_dS, & !< derivative of Li with respect to S dLi_dFi - end subroutine constitutive_LiAndItsTangents + end subroutine phase_LiAndItsTangents module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & @@ -186,7 +186,7 @@ contains !> @brief Initialize mechanical field related constitutive models !> @details Initialize elasticity, plasticity and stiffness degradation models. !-------------------------------------------------------------------------------------------------- -module subroutine mech_init(phases) +module subroutine mechanical_init(phases) class(tNode), pointer :: & phases @@ -215,38 +215,38 @@ module subroutine mech_init(phases) allocate(phase_NstiffnessDegradations(phases%length),source=0) allocate(output_constituent(phases%length)) - allocate(constitutive_mech_Fe(phases%length)) - allocate(constitutive_mech_Fi(phases%length)) - allocate(constitutive_mech_Fi0(phases%length)) - allocate(constitutive_mech_Fp(phases%length)) - allocate(constitutive_mech_Fp0(phases%length)) - allocate(constitutive_mech_F(phases%length)) - allocate(constitutive_mech_F0(phases%length)) - allocate(constitutive_mech_Li(phases%length)) - allocate(constitutive_mech_Li0(phases%length)) - allocate(constitutive_mech_Lp0(phases%length)) - allocate(constitutive_mech_Lp(phases%length)) - allocate(constitutive_mech_S(phases%length)) - allocate(constitutive_mech_P(phases%length)) - allocate(constitutive_mech_S0(phases%length)) + allocate(phase_mechanical_Fe(phases%length)) + allocate(phase_mechanical_Fi(phases%length)) + allocate(phase_mechanical_Fi0(phases%length)) + allocate(phase_mechanical_Fp(phases%length)) + allocate(phase_mechanical_Fp0(phases%length)) + allocate(phase_mechanical_F(phases%length)) + allocate(phase_mechanical_F0(phases%length)) + allocate(phase_mechanical_Li(phases%length)) + allocate(phase_mechanical_Li0(phases%length)) + allocate(phase_mechanical_Lp0(phases%length)) + allocate(phase_mechanical_Lp(phases%length)) + allocate(phase_mechanical_S(phases%length)) + allocate(phase_mechanical_P(phases%length)) + allocate(phase_mechanical_S0(phases%length)) do ph = 1, phases%length Nconstituents = count(material_phaseAt == ph) * discretization_nIPs - allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents)) + allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents)) phase => phases%get(ph) mech => phase%get('mechanics') @@ -287,17 +287,17 @@ module subroutine mech_init(phases) ph = material_phaseAt(co,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) - constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_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) - constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3 - constitutive_mech_F0(ph)%data(1:3,1:3,me) = math_I3 + 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) + phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) & + / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal) + phase_mechanical_Fi0(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), & - constitutive_mech_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) - constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_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_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,me), & + phase_mechanical_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration + phase_mechanical_Fp(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(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) + phase_mechanical_F(ph)%data(1:3,1:3,me) = phase_mechanical_F0(ph)%data(1:3,1:3,me) enddo enddo; enddo @@ -307,14 +307,14 @@ module subroutine mech_init(phases) ! initialize plasticity allocate(plasticState(phases%length)) 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.) call plastic_init() do ph = 1, phases%length 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 num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -345,14 +345,14 @@ module subroutine mech_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 !> 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) integer, intent(in) :: & @@ -376,7 +376,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j, ph, me 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)) 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 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 integer, intent(in) :: ph @@ -407,28 +407,28 @@ module subroutine mech_results(group,ph) select case(phase_plasticity(ph)) 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) - call plastic_phenopowerlaw_results(phase_plasticityInstance(ph),group//'plastic/') + call plastic_phenopowerlaw_results(phase_plasticInstance(ph),group//'plastic/') 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) - call plastic_dislotwin_results(phase_plasticityInstance(ph),group//'plastic/') + call plastic_dislotwin_results(phase_plasticInstance(ph),group//'plastic/') 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) - call plastic_nonlocal_results(phase_plasticityInstance(ph),group//'plastic/') + call plastic_nonlocal_results(phase_plasticInstance(ph),group//'plastic/') end select 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) - Lpguess = constitutive_mech_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 + Lpguess = phase_mechanical_Lp(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) 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 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) 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 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) !* 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) if (error) return ! error - constitutive_mech_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 - constitutive_mech_Lp(ph)%data(1:3,1:3,me) = Lpguess - constitutive_mech_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 - constitutive_mech_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_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + phase_mechanical_S(ph)%data(1:3,1:3,me) = S + phase_mechanical_Lp(ph)%data(1:3,1:3,me) = Lpguess + phase_mechanical_Li(ph)%data(1:3,1:3,me) = Liguess + phase_mechanical_Fp(ph)%data(1:3,1:3,me) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize + phase_mechanical_Fi(ph)%data(1:3,1:3,me) = Fi_new + phase_mechanical_Fe(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),invFi_new) broken = .false. end function integrateStress @@ -668,9 +668,9 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul sizeDotState real(pReal) :: & zeta - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & + real(pReal), dimension(phase_plasticity_maxSizeDotState) :: & r ! state residuum - real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & + real(pReal), dimension(phase_plasticity_maxSizeDotState,2) :: & dotState @@ -796,7 +796,7 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip ph, & me, & sizeDotState - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic + real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic 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, & me, & 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) @@ -987,28 +987,28 @@ subroutine crystallite_results(group,ph) select case (output_constituent(ph)%label(ou)) 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') 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') 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') 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') 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') 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') 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') 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') case('O') select case(lattice_structure(ph)) @@ -1067,43 +1067,43 @@ end subroutine crystallite_results !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- -module subroutine mech_windForward(ph,me) +module subroutine mechanical_windForward(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) - constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = constitutive_mech_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) - constitutive_mech_Li0(ph)%data(1:3,1:3,me) = constitutive_mech_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) - constitutive_mech_S0(ph)%data(1:3,1:3,me) = constitutive_mech_S(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) + phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = phase_mechanical_Fi(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) + phase_mechanical_Li0(ph)%data(1:3,1:3,me) = phase_mechanical_Li(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) + 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) -end subroutine mech_windForward +end subroutine mechanical_windForward !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -module subroutine mech_forward() +module subroutine mechanical_forward() integer :: ph do ph = 1, size(plasticState) - constitutive_mech_Fi0(ph) = constitutive_mech_Fi(ph) - constitutive_mech_Fp0(ph) = constitutive_mech_Fp(ph) - constitutive_mech_F0(ph) = constitutive_mech_F(ph) - constitutive_mech_Li0(ph) = constitutive_mech_Li(ph) - constitutive_mech_Lp0(ph) = constitutive_mech_Lp(ph) - constitutive_mech_S0(ph) = constitutive_mech_S(ph) + phase_mechanical_Fi0(ph) = phase_mechanical_Fi(ph) + phase_mechanical_Fp0(ph) = phase_mechanical_Fp(ph) + phase_mechanical_F0(ph) = phase_mechanical_F(ph) + phase_mechanical_Li0(ph) = phase_mechanical_Li(ph) + phase_mechanical_Lp0(ph) = phase_mechanical_Lp(ph) + phase_mechanical_S0(ph) = phase_mechanical_S(ph) plasticState(ph)%state0 = plasticState(ph)%state enddo -end subroutine mech_forward +end subroutine mechanical_forward @@ -1111,19 +1111,19 @@ end subroutine mech_forward !> @brief returns the homogenize elasticity matrix !> 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 integer, intent(in) :: ph, me - plasticityType: select case (phase_plasticity(ph)) - case (PLASTICITY_DISLOTWIN_ID) plasticityType + plasticType: select case (phase_plasticity(ph)) + case (PLASTICITY_DISLOTWIN_ID) plasticType C = plastic_dislotwin_homogenizedC(ph,me) - case default plasticityType + case default plasticType 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) sizeDotState = plasticState(ph)%sizeDotState - subLi0 = constitutive_mech_Li0(ph)%data(1:3,1:3,me) - subLp0 = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) + subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,me) + subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%State0(:,me) do so = 1, phase_Nsources(ph) damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) enddo - subFp0 = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - subFi0 = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - subF0 = constitutive_mech_F0(ph)%data(1:3,1:3,me) + subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) + subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,me) + subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,me) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. @@ -1186,10 +1186,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) if (todo) then subF0 = subF - subLp0 = constitutive_mech_Lp(ph)%data(1:3,1:3,me) - subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me) - subFp0 = constitutive_mech_Fp(ph)%data(1:3,1:3,me) - subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + subLp0 = phase_mechanical_Lp(ph)%data(1:3,1:3,me) + subLi0 = phase_mechanical_Li(ph)%data(1:3,1:3,me) + subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,me) + subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) do so = 1, phase_Nsources(ph) 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) else subStep = num%subStepSizeCryst * subStep - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = subFp0 - constitutive_mech_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_Fp(ph)%data(1:3,1:3,me) = subFp0 + phase_mechanical_Fi(ph)%data(1:3,1:3,me) = subFi0 + 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 - constitutive_mech_Lp(ph)%data(1:3,1:3,me) = subLp0 - constitutive_mech_Li(ph)%data(1:3,1:3,me) = subLi0 + phase_mechanical_Lp(ph)%data(1:3,1:3,me) = subLp0 + phase_mechanical_Li(ph)%data(1:3,1:3,me) = subLi0 endif plasticState(ph)%state(:,me) = subState0 do so = 1, phase_Nsources(ph) @@ -1218,9 +1218,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ! prepare for integration if (todo) then subF = subF0 & - + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_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), & - constitutive_mech_Fp(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)) + 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), & + 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_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) endif @@ -1233,7 +1233,7 @@ end function crystallite_stress !-------------------------------------------------------------------------------------------------- !> @brief Restore data after homog cutback. !-------------------------------------------------------------------------------------------------- -module subroutine mech_restore(ce,includeL) +module subroutine mechanical_restore(ce,includeL) integer, intent(in) :: ce logical, intent(in) :: & @@ -1247,23 +1247,23 @@ module subroutine mech_restore(ce,includeL) ph = material_phaseAt2(co,ce) me = material_phaseMemberAt2(co,ce) if (includeL) then - constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_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_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(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 - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_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) - constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(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) + phase_mechanical_Fi(ph)%data(1:3,1:3,me) = phase_mechanical_Fi0(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) enddo -end subroutine mech_restore +end subroutine mechanical_restore !-------------------------------------------------------------------------------------------------- !> @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 integer, intent(in) :: & @@ -1297,18 +1297,18 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & - constitutive_mech_Fe(ph)%data(1:3,1:3,me), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) - call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & - constitutive_mech_S(ph)%data(1:3,1:3,me), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me), & + call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & + phase_mechanical_Fe(ph)%data(1:3,1:3,me), & + phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el) + call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & + phase_mechanical_S(ph)%data(1:3,1:3,me), & + phase_mechanical_Fi(ph)%data(1:3,1:3,me), & co,ip,el) - invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) - invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) - invSubFp0 = math_inv33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me)) - invSubFi0 = math_inv33(constitutive_mech_Fi0(ph)%data(1:3,1:3,me)) + invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,me)) + invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,me)) + invSubFp0 = math_inv33(phase_mechanical_Fp0(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 dFidS = 0.0_pReal @@ -1334,15 +1334,15 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) endif call plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & - constitutive_mech_S(ph)%data(1:3,1:3,me), & - constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el) + phase_mechanical_S(ph)%data(1:3,1:3,me), & + phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el) dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- ! calculate dSdF temp_33_1 = transpose(matmul(invFp,invFi)) - temp_33_2 = matmul(constitutive_mech_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_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invSubFp0) + 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 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 - temp_33_1 = matmul(constitutive_mech_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_3 = matmul(temp_33_2,constitutive_mech_S(ph)%data(1:3,1:3,me)) + temp_33_1 = matmul(phase_mechanical_S(ph)%data(1:3,1:3,me),transpose(invFp)) + temp_33_2 = matmul(phase_mechanical_F(ph)%data(1:3,1:3,me),invFp) + temp_33_3 = matmul(temp_33_2,phase_mechanical_S(ph)%data(1:3,1:3,me)) dPdF = 0.0_pReal do p=1,3 @@ -1380,129 +1380,129 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) enddo do o=1,3; do p=1,3 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(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o))) 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, intent(in) :: ph call HDF5_write(groupHandle,plasticState(ph)%state,'omega') - call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,'F_i') - call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,'L_i') - call HDF5_write(groupHandle,constitutive_mech_Lp(ph)%data,'L_p') - call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,'F_p') - call HDF5_write(groupHandle,constitutive_mech_S(ph)%data,'S') - call HDF5_write(groupHandle,constitutive_mech_F(ph)%data,'F') + call HDF5_write(groupHandle,phase_mechanical_Fi(ph)%data,'F_i') + call HDF5_write(groupHandle,phase_mechanical_Li(ph)%data,'L_i') + call HDF5_write(groupHandle,phase_mechanical_Lp(ph)%data,'L_p') + call HDF5_write(groupHandle,phase_mechanical_Fp(ph)%data,'F_p') + call HDF5_write(groupHandle,phase_mechanical_S(ph)%data,'S') + 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, intent(in) :: ph call HDF5_read(groupHandle,plasticState(ph)%state0,'omega') - call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,'F_i') - call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,'L_i') - call HDF5_read(groupHandle,constitutive_mech_Lp0(ph)%data,'L_p') - call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,'F_p') - call HDF5_read(groupHandle,constitutive_mech_S0(ph)%data,'S') - call HDF5_read(groupHandle,constitutive_mech_F0(ph)%data,'F') + call HDF5_read(groupHandle,phase_mechanical_Fi0(ph)%data,'F_i') + call HDF5_read(groupHandle,phase_mechanical_Li0(ph)%data,'L_i') + call HDF5_read(groupHandle,phase_mechanical_Lp0(ph)%data,'L_p') + call HDF5_read(groupHandle,phase_mechanical_Fp0(ph)%data,'F_p') + call HDF5_read(groupHandle,phase_mechanical_S0(ph)%data,'S') + 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) !---------------------------------------------------------------------------------------------- -module function mech_S(ph,me) result(S) +module function mechanical_S(ph,me) result(S) integer, intent(in) :: ph,me 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) !---------------------------------------------------------------------------------------------- -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 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) !---------------------------------------------------------------------------------------------- -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 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) !---------------------------------------------------------------------------------------------- -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 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) !---------------------------------------------------------------------------------------------- -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 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 -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 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 diff --git a/src/phase_mechanics_eigendeformation.f90 b/src/phase_mechanical_eigen.f90 similarity index 96% rename from src/phase_mechanics_eigendeformation.f90 rename to src/phase_mechanical_eigen.f90 index 45cfd82d3..6df993565 100644 --- a/src/phase_mechanics_eigendeformation.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -126,7 +126,7 @@ end function kinematics_active !> @brief contains the constitutive equation for calculating the velocity gradient ! 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) integer, intent(in) :: & @@ -159,15 +159,15 @@ module subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & dLi_dS = 0.0_pReal dLi_dFi = 0.0_pReal - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - case (PLASTICITY_isotropic_ID) plasticityType + plasticType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_isotropic_ID) plasticType 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) - case default plasticityType + case default plasticType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal - end select plasticityType + end select plasticType Li = Li + my_Li 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) enddo; enddo -end subroutine constitutive_LiAndItsTangents +end subroutine phase_LiAndItsTangents end submodule eigendeformation diff --git a/src/phase_mechanics_eigendeformation_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 similarity index 100% rename from src/phase_mechanics_eigendeformation_cleavageopening.f90 rename to src/phase_mechanical_eigen_cleavageopening.f90 diff --git a/src/phase_mechanics_eigendeformation_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 similarity index 100% rename from src/phase_mechanics_eigendeformation_slipplaneopening.f90 rename to src/phase_mechanical_eigen_slipplaneopening.f90 diff --git a/src/phase_mechanics_eigendeformation_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 similarity index 100% rename from src/phase_mechanics_eigendeformation_thermalexpansion.f90 rename to src/phase_mechanical_eigen_thermalexpansion.f90 diff --git a/src/phase_mechanics_plastic.f90 b/src/phase_mechanical_plastic.f90 similarity index 91% rename from src/phase_mechanics_plastic.f90 rename to src/phase_mechanical_plastic.f90 index fff1fbdcf..ed1dc64fb 100644 --- a/src/phase_mechanics_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -270,31 +270,31 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & me = material_phasememberAt(co,ip,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 dLp_dMp = 0.0_pReal - case (PLASTICITY_ISOTROPIC_ID) plasticityType + case (PLASTICITY_ISOTROPIC_ID) plasticType 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) - case (PLASTICITY_KINEHARDENING_ID) plasticityType + case (PLASTICITY_KINEHARDENING_ID) plasticType 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) - case (PLASTICITY_DISLOTWIN_ID) plasticityType + case (PLASTICITY_DISLOTWIN_ID) plasticType 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) - end select plasticityType + end select plasticType 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))) + & @@ -323,29 +323,29 @@ module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) logical :: broken - Mp = matmul(matmul(transpose(constitutive_mech_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)) + Mp = matmul(matmul(transpose(phase_mechanical_Fi(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) - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + case (PLASTICITY_PHENOPOWERLAW_ID) plasticType call phenopowerlaw_dotState(Mp,ph,me) - case (PLASTICITY_KINEHARDENING_ID) plasticityType + case (PLASTICITY_KINEHARDENING_ID) plasticType 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) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType 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) - end select plasticityType + end select plasticType 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) 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) - case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType + case (PLASTICITY_DISLOTUNGSTEN_ID) plasticType call dislotungsten_dependentState(instance,me) - case (PLASTICITY_NONLOCAL_ID) plasticityType + case (PLASTICITY_NONLOCAL_ID) plasticType call nonlocal_dependentState(instance,me,ip,el) - end select plasticityType + end select plasticType end subroutine plastic_dependentState @@ -410,24 +410,24 @@ module function plastic_deltaState(co, ip, el, ph, me) result(broken) mySize - Mp = matmul(matmul(transpose(constitutive_mech_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)) - instance = phase_plasticityInstance(ph) + Mp = matmul(matmul(transpose(phase_mechanical_Fi(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_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) 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) broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,me))) case default broken = .false. - end select plasticityType + end select plasticType if(.not. broken) then select case(phase_plasticity(ph)) diff --git a/src/phase_mechanics_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 similarity index 97% rename from src/phase_mechanics_plastic_dislotungsten.f90 rename to src/phase_mechanical_plastic_dislotungsten.f90 index a47132c63..a9946b6e7 100644 --- a/src/phase_mechanics_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -226,7 +226,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization @@ -289,16 +289,16 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & integer :: & 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, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg Lp = 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 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) & @@ -327,7 +327,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) real(pReal) :: & 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,& tau_pos,& tau_neg, & @@ -336,10 +336,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,me) dot_rho_dip_climb, & dip_distance - associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& - dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(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, & tau_pos_out = tau_pos,tau_neg_out = tau_neg) diff --git a/src/phase_mechanics_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 similarity index 97% rename from src/phase_mechanics_plastic_dislotwin.f90 rename to src/phase_mechanical_plastic_dislotwin.f90 index 9f7464323..e1259fbd2 100644 --- a/src/phase_mechanics_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -415,7 +415,7 @@ module function plastic_dislotwin_init() result(myPlasticity) 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 @@ -496,8 +496,8 @@ module function plastic_dislotwin_homogenizedC(ph,me) result(homogenizedC) real(pReal) :: f_unrotated - 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 & - 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, & ddot_gamma_dtau, & 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 - 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 - 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 real(pReal):: dot_gamma_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 & ],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 & - 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 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 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) & @@ -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) 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 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) & @@ -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) 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 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) & @@ -653,24 +653,24 @@ module subroutine dislotwin_dotState(Mp,T,ph,me) tau, & sigma_cl, & !< climb stress 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_climb, & rho_dip_distance_min, & 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 - real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_tr) :: & + real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_tr) :: & dot_gamma_tr - associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & - dot => dotState(phase_plasticityInstance(ph)), dst => dependentState(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), & + dot => dotState(phase_plasticInstance(ph)), dst => dependentState(phase_plasticInstance(ph))) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,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) 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) & - 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 - 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 end associate diff --git a/src/phase_mechanics_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 similarity index 97% rename from src/phase_mechanics_plastic_isotropic.f90 rename to src/phase_mechanical_plastic_isotropic.f90 index 6789e74b4..1a4e1449a 100644 --- a/src/phase_mechanics_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -135,7 +135,7 @@ module function plastic_isotropic_init() result(myPlasticity) sizeDotState = size(['xi ','gamma']) sizeState = sizeDotState - call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization @@ -190,7 +190,7 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) integer :: & 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) 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 norm_Mp !< norm of the (deviatoric) Mandel stress - associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & - dot => dotState(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), & + dot => dotState(phase_plasticInstance(ph))) if (prm%dilatation) then norm_Mp = sqrt(math_tensordot(Mp,Mp)) diff --git a/src/phase_mechanics_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 similarity index 97% rename from src/phase_mechanics_plastic_kinehardening.f90 rename to src/phase_mechanical_plastic_kinehardening.f90 index 919302bc1..f00e4171b 100644 --- a/src/phase_mechanics_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -180,7 +180,7 @@ module function plastic_kinehardening_init() result(myPlasticity) sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names 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 @@ -255,16 +255,16 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) integer :: & 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, & dgdot_dtau_pos,dgdot_dtau_neg Lp = 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 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) :: & 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 - associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)),& - dot => dotState(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(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) sumGamma = sum(stt%accshear(:,me)) diff --git a/src/phase_mechanics_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 similarity index 96% rename from src/phase_mechanics_plastic_none.f90 rename to src/phase_mechanical_plastic_none.f90 index b8d1678eb..b95eaa039 100644 --- a/src/phase_mechanics_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -44,7 +44,7 @@ module function plastic_none_init() result(myPlasticity) phase => phases%get(p) if(.not. myPlasticity(p)) cycle Nconstituents = count(material_phaseAt2 == p) - call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0) + call phase_allocateState(plasticState(p),Nconstituents,0,0,0) enddo end function plastic_none_init diff --git a/src/phase_mechanics_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 similarity index 97% rename from src/phase_mechanics_plastic_nonlocal.f90 rename to src/phase_mechanical_plastic_nonlocal.f90 index 693ffcf93..724086916 100644 --- a/src/phase_mechanics_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -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 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') 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) if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then ph = material_phaseAt(1,el) - invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) - invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me)) + invFp = math_inv33(phase_mechanical_Fp(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_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) no = material_phasememberAt(1,neighbor_ip,neighbor_el) 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 nRealNeighbors = nRealNeighbors + 1.0_pReal @@ -782,25 +782,25 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & l, & t, & !< dislocation type 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) - real(pReal), dimension(param(phase_plasticityInstance(ph))%sum_N_sl,10) :: & + real(pReal), dimension(param(phase_plasticInstance(ph))%sum_N_sl,10) :: & 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 tauNS, & !< resolved shear stress including non Schmid and backstress terms dv_dtau, & !< 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 gdotTotal !< shear rate - associate(prm => param(phase_plasticityInstance(ph)),dst=>microstructure(phase_plasticityInstance(ph)),& - stt=>state(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)),dst=>microstructure(phase_plasticInstance(ph)),& + stt=>state(phase_plasticInstance(ph))) ns = prm%sum_N_sl !*** shortcut to state variables - rho = getRho(phase_plasticityInstance(ph),me,ip,el) + rho = getRho(phase_plasticInstance(ph),me,ip,el) rhoSgl = rho(:,sgl) do s = 1,ns @@ -820,7 +820,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & ! edges 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) dv_dtau(:,2) = dv_dtau(:,1) dv_dtauNS(:,2) = dv_dtauNS(:,1) @@ -833,7 +833,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, & else do t = 3,4 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 endif @@ -992,7 +992,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & c, & !< character of dislocation t, & !< type of dislocation 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, & rho0, & !< dislocation density at beginning of time step rhoDot, & !< density evolution @@ -1000,17 +1000,17 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, & rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) rhoDotAthermalAnnihilation, & !< density evolution by athermal 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) 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 v0, & 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 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) dLower, & !< minimum 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 endif - associate(prm => param(phase_plasticityInstance(ph)), & - dst => microstructure(phase_plasticityInstance(ph)), & - dot => dotState(phase_plasticityInstance(ph)), & - stt => state(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)), & + dst => microstructure(phase_plasticInstance(ph)), & + dot => dotState(phase_plasticInstance(ph)), & + stt => state(phase_plasticInstance(ph))) ns = prm%sum_N_sl tau = 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) rhoDip = rho(:,dip) - rho0 = getRho0(phase_plasticityInstance(ph),me,ip,el) + rho0 = getRho0(phase_plasticInstance(ph),me,ip,el) 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) #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) 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) & - 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 & + rhoDotSingle2DipoleGlide & + rhoDotAthermalAnnihilation & @@ -1284,8 +1284,8 @@ function rhoDotFlux(timestep,instance,me,ip,el) m(1:3,:,3) = -prm%slip_transverse m(1:3,:,4) = prm%slip_transverse - my_F = constitutive_mech_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_F = phase_mechanical_F(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 @@ -1301,9 +1301,9 @@ function rhoDotFlux(timestep,instance,me,ip,el) opposite_n = IPneighborhood(3,opposite_neighbor,ip,el) if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient - neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el)) - neighbor_F = constitutive_mech_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_instance = phase_plasticInstance(material_phaseAt(1,neighbor_el)) + neighbor_F = phase_mechanical_F(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) else ! if no neighbor, take my value as average Favg = my_F diff --git a/src/phase_mechanics_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 similarity index 96% rename from src/phase_mechanics_plastic_phenopowerlaw.f90 rename to src/phase_mechanical_plastic_phenopowerlaw.f90 index ea2f53100..f9791faa6 100644 --- a/src/phase_mechanics_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -231,7 +231,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) sizeState = sizeDotState - call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) + call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization @@ -300,18 +300,18 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,me) integer :: & 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, & 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 Lp = 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 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) & @@ -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) 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 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) & @@ -348,12 +348,12 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me) c_SlipSlip,c_TwinSlip,c_TwinTwin, & xi_slip_sat_offset,& 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, & gdot_slip_pos,gdot_slip_neg - associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph)), & - dot => dotState(phase_plasticityInstance(ph))) + associate(prm => param(phase_plasticInstance(ph)), stt => state(phase_plasticInstance(ph)), & + dot => dotState(phase_plasticInstance(ph))) sumGamma = sum(stt%gamma_slip(:,me)) sumF = sum(stt%gamma_twin(:,me)/prm%gamma_tw_char) @@ -373,9 +373,9 @@ module subroutine phenopowerlaw_dotState(Mp,ph,me) !-------------------------------------------------------------------------------------------------- ! 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) - 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 diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index d5d52b010..85ca2f7a5 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -124,7 +124,7 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @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 real(pReal), intent(out) :: & @@ -153,13 +153,13 @@ module subroutine constitutive_thermal_getRate(TDot, ph,me) enddo -end subroutine constitutive_thermal_getRate +end subroutine phase_thermal_getRate !-------------------------------------------------------------------------------------------------- !> @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 logical :: broken @@ -178,7 +178,7 @@ function constitutive_thermal_collectDotState(ph,me) result(broken) enddo SourceLoop -end function constitutive_thermal_collectDotState +end function phase_thermal_collectDotState module function thermal_stress(Delta_t,ph,me) result(converged_) @@ -207,7 +207,7 @@ function integrateThermalState(Delta_t, ph,me) result(broken) so, & sizeDotState - broken = constitutive_thermal_collectDotState(ph,me) + broken = phase_thermal_collectDotState(ph,me) if(broken) return do so = 1, thermal_Nsources(ph) @@ -264,7 +264,7 @@ end function thermal_dot_T !---------------------------------------------------------------------------------------------- !< @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 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))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T -end subroutine constitutive_thermal_setField +end subroutine phase_thermal_setField diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index 7a043308d..3c0095401 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -56,7 +56,7 @@ module function dissipation_init(source_length) result(mySources) prm%kappa = src%get_asFloat('kappa') 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 endif @@ -78,7 +78,7 @@ module subroutine dissipation_getRate(TDot, ph,me) 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 subroutine dissipation_getRate diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index ddae6763b..d2874ded0 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -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)) 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 endif enddo