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

Merge remote-tracking branch 'origin/development' into initial-eigenstrain

This commit is contained in:
Martin Diehl 2021-08-05 16:25:27 +02:00
parent 5503592790 d8c03ad0d5
commit 00a053c4c8
37 changed files with 934 additions and 712 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378 Subproject commit 9699f20f21f8a5f532c735a1aa9daeba395da94d

42
examples/config/Phase_DisloUCLA_Tungsten.config
View File

@ -1,42 +0,0 @@
[Tungsten]
elasticity hooke
plasticity disloucla
(output) edge_density
(output) dipole_density
(output) shear_rate_slip
(output) accumulated_shear_slip
(output) resolved_stress_slip
(output) threshold_stress_slip
grainsize 2.7e-5 # Average grain size [m] 2.0e-5
SolidSolutionStrength 0.0 # Strength due to elements in solid solution
### Dislocation glide parameters ###
#per family
Nslip 12
slipburgers 2.72e-10 # Burgers vector of slip system [m]
rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3]
rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3]
Qedge 2.61154e-19 # Activation energy for dislocation glide [J], 1.63 eV
v0 1 # Initial glide velocity [m/s]
p_slip 0.86 # p-exponent in glide velocity
q_slip 1.69 # q-exponent in glide velocity
tau_peierls 2.03e9 # peierls stress [Pa]
#mobility law
kink_height 2.567e-10 # kink height sqrt(6)/3*lattice_parameter [m]
omega 9.1e11 # attemp frequency (from kMC paper) [s^(-1)]
kink_width 29.95e-10 # kink pair width ~ 11 b (kMC paper) [m]
dislolength 78e-10 # dislocation length (ideally lambda) [m] initial value 11b
friction_coeff 8.3e-5 # friction coeff. B [Pa*s]
#hardening
dipoleformationfactor 0 # to have hardening due to dipole formation off
CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path
D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interaction_slipslip 0.009 0.72 0.009 0.05 0.05 0.06 0.09
nonschmid_coefficients 0.938 0.71 4.43 0.0 0.0 0.0

26
examples/config/Phase_Dislotungsten_W.yaml Normal file
View File

@ -0,0 +1,26 @@
type: dislotungsten
N_sl: [12]
rho_mob_0: [1.0e+9]
rho_dip_0: [1.0]
nu_a: [9.1e+11]
b_sl: [2.72e-10]
Delta_H_kp,0: [2.61154e-19] # 1.63 eV, Delta_H0
tau_Peierls: [2.03e+9]
p_sl: [0.86]
q_sl: [1.69]
h: [2.566e-10]
w: [2.992e-09]
B: [8.3e-5]
D_a: 1.0 # d_edge
# climb (disabled)
D_0: 0.0
Q_cl: 0.0
V_cl: [0.0]
h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09]
a_nonSchmid: [0.938, 0.71, 4.43]

2
examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml
View File

@ -6,7 +6,7 @@ b_sl: [2.56e-10]
rho_mob_0: [1.0e+12] rho_mob_0: [1.0e+12]
rho_dip_0: [1.0] rho_dip_0: [1.0]
v_0: [1.0e+4] v_0: [1.0e+4]
Q_s: [3.7e-19] Q_sl: [3.7e-19]
p_sl: [1.0] p_sl: [1.0]
q_sl: [1.0] q_sl: [1.0]
tau_0: [1.5e+8] tau_0: [1.5e+8]

3
examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
View File

@ -11,7 +11,7 @@ b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e12, 2.8e+12] rho_mob_0: [2.81e12, 2.8e+12]
rho_dip_0: [1.0, 1.0] # not given rho_dip_0: [1.0, 1.0] # not given
v_0: [1.4e+3, 1.4e+3] v_0: [1.4e+3, 1.4e+3]
Q_s: [1.57e-19, 1.57e-19] # Delta_F Q_sl: [1.57e-19, 1.57e-19] # Delta_F
tau_0: [454.e+6, 454.e+6] tau_0: [454.e+6, 454.e+6]
p_sl: [0.325, 0.325] p_sl: [0.325, 0.325]
q_sl: [1.55, 1.55] q_sl: [1.55, 1.55]
@ -19,6 +19,5 @@ i_sl: [23.3, 23.3]
D_a: 7.4 # C_anni D_a: 7.4 # C_anni
B: [0.001, 0.001] B: [0.001, 0.001]
h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073] h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073]
D_0: 4.0e-05
Q_cl: 5.4e-19 # no recovery! Q_cl: 5.4e-19 # no recovery!
D: 40.e-6 # estimated D: 40.e-6 # estimated

33
examples/config/Phase_Phenopowerlaw_Magnesium.yaml → examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
View File

@ -1,20 +1,29 @@
N_sl: [3, 3, 0, 6, 0, 6]
N_tw: [6, 0, 0, 6]
h_0_tw-tw: 50.0e+6
h_0_sl-sl: 500.0e+6
h_0_tw-sl: 150.0e+6
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
output: [xi_sl, xi_tw]
type: phenopowerlaw type: phenopowerlaw
xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] references:
- F. Wang et al.,
Acta Materialia 80:77-93, 2014,
https://doi.org/10.1016/j.actamat.2014.07.048
output: [xi_sl, xi_tw]
N_sl: [3, 3, 0, 6, 0, 6] # basal, 1. prism, -, 1. pyr<a>, -, 2. pyr<c+a>
N_tw: [6, 0, 6] # tension, -, compression
xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6]
xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6] xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6]
xi_0_tw: [40.e+6, 0., 0., 60.e+6] xi_0_tw: [40.e+6, 0., 60.e+6]
a_sl: 2.25 a_sl: 2.25
dot_gamma_0_sl: 0.001 dot_gamma_0_sl: 0.001
dot_gamma_0_tw: 0.001 dot_gamma_0_tw: 0.001
n_sl: 20 n_sl: 20
n_tw: 20 n_tw: 20
f_sat_sl-tw: 10.0 f_sat_sl-tw: 10.0
h_0_sl-sl: 500.0e+6
h_0_tw-tw: 50.0e+6
h_0_tw-sl: 150.0e+6
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

5
examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
View File

@ -7,14 +7,17 @@ references:
Acta Materialia 132:598-610, 2017, Acta Materialia 132:598-610, 2017,
https://doi.org/10.1016/j.actamat.2017.05.015 https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl] output: [gamma_sl]
N_sl: [3, 3, 0, 0, 12]
N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 2. pyr<c+a>
n_sl: 20 n_sl: 20
a_sl: 2.0 a_sl: 2.0
dot_gamma_0_sl: 0.001 dot_gamma_0_sl: 0.001
h_0_sl-sl: 200.e+6 h_0_sl-sl: 200.e+6
# C. Zambaldi et al.: # C. Zambaldi et al.:
xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6] xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6]
xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6] xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6]
# L. Wang et al. : # L. Wang et al. :
# xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6] # xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6]
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

4
examples/mesh/tensionZ.load
View File

@ -1,11 +1,11 @@
#initial elastic step #initial elastic step
$Loadcase 1 time 0.0005 incs 1 frequency 5 $Loadcase 1 t 0.0005 N 1 f_out 5
Face 1 X 0.01 Face 1 X 0.01
Face 2 X 0.0 Face 2 X 0.0
Face 2 Y 0.0 Face 2 Y 0.0
Face 2 Z 0.0 Face 2 Z 0.0
$EndLoadcase $EndLoadcase
$Loadcase 2 time 10.0 incs 200 frequency 5 $Loadcase 2 t 10.0 N 200 f_out 5
Face 1 X 0.01 Face 1 X 0.01
Face 2 X 0.0 Face 2 X 0.0
Face 2 Y 0.0 Face 2 Y 0.0

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha4-137-gb69b85754 v3.0.0-alpha4-221-g4a8c83611

35
python/damask/_result.py
View File

@ -99,8 +99,10 @@ class Result:
self.version_major = f.attrs['DADF5_version_major'] self.version_major = f.attrs['DADF5_version_major']
self.version_minor = f.attrs['DADF5_version_minor'] self.version_minor = f.attrs['DADF5_version_minor']
if self.version_major != 0 or not 12 <= self.version_minor <= 13: if self.version_major != 0 or not 12 <= self.version_minor <= 14:
raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}') raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}')
if self.version_major == 0 and self.version_minor < 14:
self.export_setup = None
self.structured = 'cells' in f['geometry'].attrs.keys() self.structured = 'cells' in f['geometry'].attrs.keys()
@ -1395,7 +1397,7 @@ class Result:
def export_XDMF(self,output='*'): def export_XDMF(self,output='*'):
""" """
Write XDMF file to directly visualize data in DADF5 file. Write XDMF file to directly visualize data from DADF5 file.
The XDMF format is only supported for structured grids The XDMF format is only supported for structured grids
with single phase and single constituent. with single phase and single constituent.
@ -1748,3 +1750,32 @@ class Result:
if flatten: r = util.dict_flatten(r) if flatten: r = util.dict_flatten(r)
return None if (type(r) == dict and r == {}) else r return None if (type(r) == dict and r == {}) else r
def export_setup(self,output='*',overwrite=False):
"""
Export configuration files.
Parameters
----------
output : (list of) str, optional
Names of the datasets to export to the file.
Defaults to '*', in which case all datasets are exported.
overwrite : boolean, optional
Overwrite existing configuration files.
Defaults to False.
"""
def export(name,obj,output,overwrite):
if type(obj) == h5py.Dataset and _match(output,[name]):
d = obj.attrs['description'] if h5py3 else obj.attrs['description'].decode()
if not Path(name).exists() or overwrite:
with open(name,'w') as f_out: f_out.write(obj[()].decode())
print(f"Exported {d} to '{name}'.")
else:
print(f"'{name}' exists, {d} not exported.")
elif type(obj) == h5py.Group:
os.makedirs(name, exist_ok=True)
with h5py.File(self.fname,'r') as f_in:
f_in['setup'].visititems(partial(export,output=output,overwrite=overwrite))

27
src/CPFEM.f90
View File

@ -4,22 +4,23 @@
!> @brief CPFEM engine !> @brief CPFEM engine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module CPFEM module CPFEM
use DAMASK_interface
use prec use prec
use math use IO
use rotations
use YAML_types use YAML_types
use YAML_parse use YAML_parse
use discretization_marc
use material
use config
use homogenization
use IO
use discretization
use DAMASK_interface
use HDF5_utilities use HDF5_utilities
use results use results
use config
use math
use rotations
use lattice use lattice
use material
use phase use phase
use homogenization
use discretization
use discretization_marc
implicit none implicit none
private private
@ -68,7 +69,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief call all module initializations !> @brief Initialize all modules.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll subroutine CPFEM_initAll
@ -77,13 +78,13 @@ subroutine CPFEM_initAll
call IO_init call IO_init
call YAML_types_init call YAML_types_init
call YAML_parse_init call YAML_parse_init
call HDF5_utilities_init
call results_init(.false.)
call config_init call config_init
call math_init call math_init
call rotations_init call rotations_init
call HDF5_utilities_init
call results_init(.false.)
call discretization_marc_init
call lattice_init call lattice_init
call discretization_marc_init
call material_init(.false.) call material_init(.false.)
call phase_init call phase_init
call homogenization_init call homogenization_init

34
src/CPFEM2.f90
View File

@ -4,28 +4,29 @@
!> @brief needs a good name and description !> @brief needs a good name and description
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module CPFEM2 module CPFEM2
use prec
use parallelization use parallelization
use config use DAMASK_interface
use math use prec
use rotations use IO
use YAML_types use YAML_types
use YAML_parse use YAML_parse
use material
use lattice
use IO
use base64
use DAMASK_interface
use discretization
use HDF5 use HDF5
use HDF5_utilities use HDF5_utilities
use results use results
use homogenization use config
use math
use rotations
use lattice
use material
use phase use phase
use homogenization
use discretization
#if defined(MESH) #if defined(MESH)
use FEM_quadrature use FEM_quadrature
use discretization_mesh use discretization_mesh
#elif defined(GRID) #elif defined(GRID)
use base64
use discretization_grid use discretization_grid
#endif #endif
@ -36,7 +37,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief call all module initializations !> @brief Initialize all modules.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine CPFEM_initAll subroutine CPFEM_initAll
@ -44,18 +45,19 @@ subroutine CPFEM_initAll
call DAMASK_interface_init ! Spectral and FEM interface to commandline call DAMASK_interface_init ! Spectral and FEM interface to commandline
call prec_init call prec_init
call IO_init call IO_init
call base64_init #if defined(MESH)
#ifdef MESH
call FEM_quadrature_init call FEM_quadrature_init
#elif defined(GRID)
call base64_init
#endif #endif
call YAML_types_init call YAML_types_init
call YAML_parse_init call YAML_parse_init
call HDF5_utilities_init
call results_init(restart=interface_restartInc>0)
call config_init call config_init
call math_init call math_init
call rotations_init call rotations_init
call lattice_init call lattice_init
call HDF5_utilities_init
call results_init(restart=interface_restartInc>0)
#if defined(MESH) #if defined(MESH)
call discretization_mesh_init(restart=interface_restartInc>0) call discretization_mesh_init(restart=interface_restartInc>0)
#elif defined(GRID) #elif defined(GRID)

6
src/DAMASK_interface.f90
View File

@ -139,10 +139,10 @@ subroutine DAMASK_interface_init
print'(a)', ' Optional arguments:' print'(a)', ' Optional arguments:'
print'(/,a)',' --workingdirectory PathToWorkingDirectory' print'(/,a)',' --workingdirectory PathToWorkingDirectory'
print'(a)', ' Specifies the working directory and overwrites the default ./' print'(a)', ' Specifies the working directory and overwrites the default ./'
print'(a)', ' Make sure the file "material.config" exists in the working' print'(a)', ' Make sure the file "material.yaml" exists in the working'
print'(a)', ' directory.' print'(a)', ' directory.'
print'(a)', ' For further configuration place "numerics.config"' print'(a)', ' For further configuration place "numerics.yaml"'
print'(a)',' and "debug.config" in that directory.' print'(a)',' and "debug.yaml" in that directory.'
print'(/,a)',' --restart N' print'(/,a)',' --restart N'
print'(a)', ' Reads in increment N and continues with calculating' print'(a)', ' Reads in increment N and continues with calculating'
print'(a)', ' increment N+1 based on this.' print'(a)', ' increment N+1 based on this.'

54
src/HDF5_utilities.f90
View File

@ -85,6 +85,7 @@ module HDF5_utilities
HDF5_utilities_init, & HDF5_utilities_init, &
HDF5_read, & HDF5_read, &
HDF5_write, & HDF5_write, &
HDF5_write_str, &
HDF5_addAttribute, & HDF5_addAttribute, &
HDF5_addGroup, & HDF5_addGroup, &
HDF5_openGroup, & HDF5_openGroup, &
@ -128,10 +129,11 @@ end subroutine HDF5_utilities_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief open and initializes HDF5 output file !> @brief open and initializes HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(HID_T) function HDF5_openFile(fileName,mode) integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
character(len=*), intent(in) :: fileName character(len=*), intent(in) :: fileName
character, intent(in), optional :: mode character, intent(in), optional :: mode
logical, intent(in), optional :: parallel
character :: m character :: m
integer(HID_T) :: plist_id integer(HID_T) :: plist_id
@ -148,7 +150,11 @@ integer(HID_T) function HDF5_openFile(fileName,mode)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
#ifdef PETSC #ifdef PETSC
call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr) if (present(parallel)) then
if (parallel) call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
else
call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
endif
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
#endif #endif
@ -1467,6 +1473,48 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel)
end subroutine HDF5_write_real7 end subroutine HDF5_write_real7
!--------------------------------------------------------------------------------------------------
!> @brief Write dataset of type string (scalar).
!> @details Not collective, must be called by one process at at time.
!--------------------------------------------------------------------------------------------------
subroutine HDF5_write_str(dataset,loc_id,datasetName)
character(len=*), intent(in) :: dataset
integer(HID_T), intent(in) :: loc_id
character(len=*), intent(in) :: datasetName !< name of the dataset in the file
integer(HID_T) :: filetype_id, space_id, dataset_id
integer :: hdferr
character(len=len_trim(dataset)+1,kind=C_CHAR), dimension(1), target :: dataset_
type(C_PTR), target, dimension(1) :: ptr
dataset_(1) = trim(dataset)//C_NULL_CHAR
ptr(1) = c_loc(dataset_(1))
call h5tcopy_f(H5T_STRING, filetype_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5tset_size_f(filetype_id, int(len(dataset_),HSIZE_T), hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5screate_f(H5S_SCALAR_F, space_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5dcreate_f(loc_id, datasetName, H5T_STRING, space_id, dataset_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5dwrite_f(dataset_id, H5T_STRING, c_loc(ptr), hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5dclose_f(dataset_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5sclose_f(space_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
call h5tclose_f(filetype_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error'
end subroutine HDF5_write_str
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief write dataset of type integer with 1 dimension !> @brief write dataset of type integer with 1 dimension
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1872,7 +1920,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! creating a property list for transfer properties (is collective when reading in parallel) ! creating a property list for transfer properties (is collective when writing in parallel)
call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr) call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
if(hdferr < 0) error stop 'HDF5 error' if(hdferr < 0) error stop 'HDF5 error'
#ifdef PETSC #ifdef PETSC

11
src/IO.f90
View File

@ -119,27 +119,28 @@ function IO_read(fileName) result(fileContent)
character(len=:), allocatable :: fileContent character(len=:), allocatable :: fileContent
integer :: & integer :: &
fileLength, &
fileUnit, & fileUnit, &
myStat myStat
integer(pI64) :: &
fileLength
inquire(file = fileName, size=fileLength) inquire(file = fileName, size=fileLength)
open(newunit=fileUnit, file=fileName, access='stream',& open(newunit=fileUnit, file=fileName, access='stream',&
status='old', position='rewind', action='read',iostat=myStat) status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName)) if (myStat /= 0) call IO_error(100,ext_msg=trim(fileName))
allocate(character(len=fileLength)::fileContent) allocate(character(len=fileLength)::fileContent)
if(fileLength==0) then if (fileLength==0) then
close(fileUnit) close(fileUnit)
return return
endif endif
read(fileUnit,iostat=myStat) fileContent read(fileUnit,iostat=myStat) fileContent
if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName)) if (myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
close(fileUnit) close(fileUnit)
if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent) if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent)
if(fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF if (fileContent(fileLength:fileLength) /= IO_EOL) fileContent = fileContent//IO_EOL ! ensure EOL@EOF
end function IO_read end function IO_read

8
src/Marc/discretization_Marc.f90
View File

@ -216,7 +216,13 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt)
mapElemSet !< list of elements in elementSet mapElemSet !< list of elements in elementSet
inputFile = IO_readlines(trim(getSolverJobName())//trim(InputFileExtension)) call results_openJobFile
call results_writeDataset_str(IO_read(trim(getSolverJobName())//InputFileExtension), 'setup', &
trim(getSolverJobName())//InputFileExtension, &
'MSC.Marc input deck')
call results_closeJobFile
inputFile = IO_readlines(trim(getSolverJobName())//InputFileExtension)
call inputRead_fileFormat(fileFormatVersion, & call inputRead_fileFormat(fileFormatVersion, &
inputFile) inputFile)
call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, & call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, &

12
src/YAML_parse.f90
View File

@ -14,7 +14,7 @@ module YAML_parse
public :: & public :: &
YAML_parse_init, & YAML_parse_init, &
YAML_parse_file YAML_parse_str
contains contains
@ -29,16 +29,16 @@ end subroutine YAML_parse_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Parse a YAML file into a a structure of nodes. !> @brief Parse a YAML string into a a structure of nodes.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function YAML_parse_file(fname) result(node) function YAML_parse_str(str) result(node)
character(len=*), intent(in) :: fname character(len=*), intent(in) :: str
class (tNode), pointer :: node class (tNode), pointer :: node
node => parse_flow(to_flow(IO_read(fname))) node => parse_flow(to_flow(str))
end function YAML_parse_file end function YAML_parse_str
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

9
src/YAML_types.f90
View File

@ -849,7 +849,7 @@ function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
select type(self) select type(node)
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dFloat = list%as1dFloat() nodeAs1dFloat = list%as1dFloat()
@ -872,11 +872,12 @@ end function tNode_get_byKey_as1dFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Access by key and convert to float array (2D) !> @brief Access by key and convert to float array (2D)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat) function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
real(pReal), intent(in), dimension(:,:), optional :: defaultVal real(pReal), intent(in), dimension(:,:), optional :: defaultVal
integer, intent(in), dimension(2), optional :: requiredShape
real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat
@ -898,6 +899,10 @@ function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat)
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
endif endif
if (present(requiredShape)) then
if (any(requiredShape /= shape(nodeAs2dFloat))) call IO_error(146,ext_msg=k)
endif
end function tNode_get_byKey_as2dFloat end function tNode_get_byKey_as2dFloat

6
src/commercialFEM_fileList.f90
View File

@ -7,18 +7,18 @@
#include "IO.f90" #include "IO.f90"
#include "YAML_types.f90" #include "YAML_types.f90"
#include "YAML_parse.f90" #include "YAML_parse.f90"
#include "HDF5_utilities.f90"
#include "results.f90"
#include "config.f90" #include "config.f90"
#include "LAPACK_interface.f90" #include "LAPACK_interface.f90"
#include "math.f90" #include "math.f90"
#include "rotations.f90" #include "rotations.f90"
#include "lattice.f90"
#include "element.f90" #include "element.f90"
#include "HDF5_utilities.f90"
#include "results.f90"
#include "geometry_plastic_nonlocal.f90" #include "geometry_plastic_nonlocal.f90"
#include "discretization.f90" #include "discretization.f90"
#include "Marc/discretization_Marc.f90" #include "Marc/discretization_Marc.f90"
#include "material.f90" #include "material.f90"
#include "lattice.f90"
#include "phase.f90" #include "phase.f90"
#include "phase_mechanical.f90" #include "phase_mechanical.f90"
#include "phase_mechanical_elastic.f90" #include "phase_mechanical_elastic.f90"

70
src/config.f90
View File

@ -8,7 +8,8 @@ module config
use IO use IO
use YAML_parse use YAML_parse
use YAML_types use YAML_types
use results
use parallelization
implicit none implicit none
private private
@ -31,6 +32,7 @@ subroutine config_init
print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- config init -+>>>'; flush(IO_STDOUT)
call parse_material call parse_material
call parse_numerics call parse_numerics
call parse_debug call parse_debug
@ -41,15 +43,25 @@ end subroutine config_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read material.yaml or <jobname>.yaml. !> @brief Read material.yaml or <jobname>.yaml.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine parse_material subroutine parse_material()
logical :: fileExists logical :: fileExists
character(len=:), allocatable :: fileContent
inquire(file='material.yaml',exist=fileExists) inquire(file='material.yaml',exist=fileExists)
if(.not. fileExists) call IO_error(100,ext_msg='material.yaml') if(.not. fileExists) call IO_error(100,ext_msg='material.yaml')
print*, 'reading material.yaml'; flush(IO_STDOUT)
config_material => YAML_parse_file('material.yaml') if (worldrank == 0) then
print*, 'reading material.yaml'; flush(IO_STDOUT)
fileContent = IO_read('material.yaml')
call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration')
call results_closeJobFile
endif
call parallelization_bcast_str(fileContent)
config_material => YAML_parse_str(fileContent)
end subroutine parse_material end subroutine parse_material
@ -57,16 +69,28 @@ end subroutine parse_material
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read numerics.yaml. !> @brief Read numerics.yaml.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine parse_numerics subroutine parse_numerics()
logical :: fexist logical :: fileExists
character(len=:), allocatable :: fileContent
config_numerics => emptyDict config_numerics => emptyDict
inquire(file='numerics.yaml', exist=fexist)
if (fexist) then inquire(file='numerics.yaml', exist=fileExists)
print*, 'reading numerics.yaml'; flush(IO_STDOUT) if (fileExists) then
config_numerics => YAML_parse_file('numerics.yaml')
if (worldrank == 0) then
print*, 'reading numerics.yaml'; flush(IO_STDOUT)
fileContent = IO_read('numerics.yaml')
call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration')
call results_closeJobFile
endif
call parallelization_bcast_str(fileContent)
config_numerics => YAML_parse_str(fileContent)
endif endif
end subroutine parse_numerics end subroutine parse_numerics
@ -75,17 +99,29 @@ end subroutine parse_numerics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Read debug.yaml. !> @brief Read debug.yaml.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine parse_debug subroutine parse_debug()
logical :: fexist logical :: fileExists
character(len=:), allocatable :: fileContent
config_debug => emptyDict config_debug => emptyDict
inquire(file='debug.yaml', exist=fexist)
fileExists: if (fexist) then inquire(file='debug.yaml', exist=fileExists)
print*, 'reading debug.yaml'; flush(IO_STDOUT) if (fileExists) then
config_debug => YAML_parse_file('debug.yaml')
endif fileExists if (worldrank == 0) then
print*, 'reading debug.yaml'; flush(IO_STDOUT)
fileContent = IO_read('debug.yaml')
call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration')
call results_closeJobFile
endif
call parallelization_bcast_str(fileContent)
config_debug => YAML_parse_str(fileContent)
endif
end subroutine parse_debug end subroutine parse_debug

1
src/element.f90
View File

@ -151,6 +151,7 @@ module element
integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)), parameter :: IPNEIGHBOR1 = & integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)), parameter :: IPNEIGHBOR1 = &
reshape([& reshape([&
-2,-3,-1 & -2,-3,-1 &
! Note: This fix is for gfortran 9 only. gfortran 8 supports neither, gfortran > 9 both variants
#if !defined(__GFORTRAN__) #if !defined(__GFORTRAN__)
],shape(IPNEIGHBOR1)) ],shape(IPNEIGHBOR1))
#else #else

14
src/grid/DAMASK_grid.f90
View File

@ -107,6 +107,8 @@ program DAMASK_grid
step_bc, & step_bc, &
step_mech, & step_mech, &
step_discretization step_discretization
character(len=:), allocatable :: &
fileContent, fname
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
@ -127,7 +129,17 @@ program DAMASK_grid
if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter') if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter')
if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack') if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack')
config_load => YAML_parse_file(trim(interface_loadFile)) if (worldrank == 0) then
fileContent = IO_read(interface_loadFile)
fname = interface_loadFile
if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup',fname,'load case definition (grid solver)')
call results_closeJobFile
endif
call parallelization_bcast_str(fileContent)
config_load => YAML_parse_str(fileContent)
solver => config_load%get('solver') solver => config_load%get('solver')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

32
src/grid/discretization_grid.f90
View File

@ -68,11 +68,21 @@ subroutine discretization_grid_init(restart)
devNull, z, z_offset devNull, z, z_offset
integer, dimension(worldsize) :: & integer, dimension(worldsize) :: &
displs, sendcounts displs, sendcounts
character(len=:), allocatable :: &
fileContent, fname
print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- discretization_grid init -+>>>'; flush(IO_STDOUT)
if(worldrank == 0) then if(worldrank == 0) then
call readVTI(grid,geomSize,origin,materialAt_global) fileContent = IO_read(interface_geomFile)
call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
fname = interface_geomFile
if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
call results_openJobFile(parallel=.false.)
call results_writeDataset_str(fileContent,'setup',fname,'geometry definition (grid solver)')
call results_closeJobFile
else else
allocate(materialAt_global(0)) ! needed for IntelMPI allocate(materialAt_global(0)) ! needed for IntelMPI
endif endif
@ -157,7 +167,8 @@ end subroutine discretization_grid_init
!> @brief Parse vtk image data (.vti) !> @brief Parse vtk image data (.vti)
!> @details https://vtk.org/Wiki/VTK_XML_Formats !> @details https://vtk.org/Wiki/VTK_XML_Formats
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine readVTI(grid,geomSize,origin,material) subroutine readVTI(grid,geomSize,origin,material, &
fileContent)
integer, dimension(3), intent(out) :: & integer, dimension(3), intent(out) :: &
grid ! grid (across all processes!) grid ! grid (across all processes!)
@ -166,28 +177,19 @@ subroutine readVTI(grid,geomSize,origin,material)
origin ! origin (across all processes!) origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: & integer, dimension(:), intent(out), allocatable :: &
material material
character(len=*), intent(in) :: &
fileContent
character(len=:), allocatable :: fileContent, dataType, headerType character(len=:), allocatable :: dataType, headerType
logical :: inFile,inImage,gotCellData,compressed logical :: inFile,inImage,gotCellData,compressed
integer :: fileUnit, myStat
integer(pI64) :: & integer(pI64) :: &
fileLength, & !< length of the geom file (in characters)
startPos, endPos, & startPos, endPos, &
s s
grid = -1 grid = -1
geomSize = -1.0_pReal geomSize = -1.0_pReal
!--------------------------------------------------------------------------------------------------
! read raw data as stream
inquire(file = trim(interface_geomFile), size=fileLength)
open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',&
status='old', position='rewind', action='read',iostat=myStat)
if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile))
allocate(character(len=fileLength)::fileContent)
read(fileUnit) fileContent
close(fileUnit)
inFile = .false. inFile = .false.
inImage = .false. inImage = .false.
gotCelldata = .false. gotCelldata = .false.

93
src/math.f90
View File

@ -90,12 +90,12 @@ subroutine math_init
integer, dimension(:), allocatable :: randInit integer, dimension(:), allocatable :: randInit
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic num_generic
print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- math init -+>>>'; flush(IO_STDOUT)
num_generic => config_numerics%get('generic',defaultVal=emptyDict) num_generic => config_numerics%get('generic',defaultVal=emptyDict)
randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0) randomSeed = num_generic%get_asInt('random_seed', defaultVal = 0)
call random_seed(size=randSize) call random_seed(size=randSize)
allocate(randInit(randSize)) allocate(randInit(randSize))
if (randomSeed > 0) then if (randomSeed > 0) then
@ -233,7 +233,7 @@ end function math_range
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief second rank identity tensor of specified dimension !> @brief Rank two identity tensor of specified dimension.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_eye(d) pure function math_eye(d)
@ -250,21 +250,25 @@ end function math_eye
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief symmetric fourth rank identity tensor of specified dimension !> @brief Symmetric rank four identity tensor.
! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself ! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_identity4th(d) pure function math_identity4th()
real(pReal), dimension(3,3,3,3) :: math_identity4th
integer, intent(in) :: d !< tensor dimension
integer :: i,j,k,l integer :: i,j,k,l
real(pReal), dimension(d,d,d,d) :: math_identity4th
real(pReal), dimension(d,d) :: identity2nd
identity2nd = math_eye(d)
do i=1,d; do j=1,d; do k=1,d; do l=1,d #ifndef __INTEL_COMPILER
math_identity4th(i,j,k,l) = 0.5_pReal & do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k)) math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
enddo
#else
do i=1,3; do j=1,3; do k=1,3; do l=1,3
math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
#endif
end function math_identity4th end function math_identity4th
@ -331,9 +335,16 @@ pure function math_outer(A,B)
real(pReal), dimension(size(A,1),size(B,1)) :: math_outer real(pReal), dimension(size(A,1),size(B,1)) :: math_outer
integer :: i,j integer :: i,j
#ifndef __INTEL_COMPILER
do concurrent(i=1:size(A,1), j=1:size(B,1))
math_outer(i,j) = A(i)*B(j)
enddo
#else
do i=1,size(A,1); do j=1,size(B,1) do i=1,size(A,1); do j=1,size(B,1)
math_outer(i,j) = A(i)*B(j) math_outer(i,j) = A(i)*B(j)
enddo; enddo enddo; enddo
#endif
end function math_outer end function math_outer
@ -373,9 +384,16 @@ pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3) :: math_mul3333xx33 real(pReal), dimension(3,3) :: math_mul3333xx33
integer :: i,j integer :: i,j
#ifndef __INTEL_COMPILER
do concurrent(i=1:3, j=1:3)
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
enddo
#else
do i=1,3; do j=1,3 do i=1,3; do j=1,3
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3)) math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
enddo; enddo enddo; enddo
#endif
end function math_mul3333xx33 end function math_mul3333xx33
@ -390,9 +408,16 @@ pure function math_mul3333xx3333(A,B)
real(pReal), dimension(3,3,3,3), intent(in) :: B real(pReal), dimension(3,3,3,3), intent(in) :: B
real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333 real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
#ifndef __INTEL_COMPILER
do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
enddo
#else
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l)) math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo
#endif
end function math_mul3333xx3333 end function math_mul3333xx3333
@ -545,19 +570,6 @@ pure function math_symmetric33(m)
end function math_symmetric33 end function math_symmetric33
!--------------------------------------------------------------------------------------------------
!> @brief symmetrize a 6x6 matrix
!--------------------------------------------------------------------------------------------------
pure function math_symmetric66(m)
real(pReal), dimension(6,6) :: math_symmetric66
real(pReal), dimension(6,6), intent(in) :: m
math_symmetric66 = 0.5_pReal * (m + transpose(m))
end function math_symmetric66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief skew part of a 3x3 matrix !> @brief skew part of a 3x3 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -712,6 +724,7 @@ pure function math_6toSym33(v6,weighted)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i integer :: i
if(present(weighted)) then if(present(weighted)) then
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
@ -736,9 +749,16 @@ pure function math_3333to99(m3333)
integer :: i,j integer :: i,j
#ifndef __INTEL_COMPILER
do concurrent(i=1:9, j=1:9)
math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
enddo
#else
do i=1,9; do j=1,9 do i=1,9; do j=1,9
math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
enddo; enddo enddo; enddo
#endif
end function math_3333to99 end function math_3333to99
@ -753,9 +773,15 @@ pure function math_99to3333(m99)
integer :: i,j integer :: i,j
#ifndef __INTEL_COMPILER
do concurrent(i=1:9, j=1:9)
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
enddo
#else
do i=1,9; do j=1,9 do i=1,9; do j=1,9
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j) math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
enddo; enddo enddo; enddo
#endif
end function math_99to3333 end function math_99to3333
@ -775,15 +801,22 @@ pure function math_sym3333to66(m3333,weighted)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i,j integer :: i,j
if(present(weighted)) then if(present(weighted)) then
w = merge(NRMMANDEL,1.0_pReal,weighted) w = merge(NRMMANDEL,1.0_pReal,weighted)
else else
w = NRMMANDEL w = NRMMANDEL
endif endif
#ifndef __INTEL_COMPILER
do concurrent(i=1:6, j=1:6)
math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
enddo
#else
do i=1,6; do j=1,6 do i=1,6; do j=1,6
math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j)) math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
enddo; enddo enddo; enddo
#endif
end function math_sym3333to66 end function math_sym3333to66
@ -803,6 +836,7 @@ pure function math_66toSym3333(m66,weighted)
real(pReal), dimension(6) :: w real(pReal), dimension(6) :: w
integer :: i,j integer :: i,j
if(present(weighted)) then if(present(weighted)) then
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
@ -828,6 +862,7 @@ pure function math_Voigt66to3333(m66)
real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix
integer :: i,j integer :: i,j
do i=1,6; do j=1, 6 do i=1,6; do j=1, 6
math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j) math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
@ -881,10 +916,11 @@ subroutine math_eigh(w,v,error,m)
real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues
real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: v !< eigenvectors real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: v !< eigenvectors
logical, intent(out) :: error logical, intent(out) :: error
integer :: ierr integer :: ierr
real(pReal), dimension(size(m,1)**2) :: work real(pReal), dimension(size(m,1)**2) :: work
v = m ! copy matrix to input (doubles as output) array v = m ! copy matrix to input (doubles as output) array
call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr) call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr)
error = (ierr /= 0) error = (ierr /= 0)
@ -1222,7 +1258,7 @@ subroutine selfTest
error stop 'math_sym33to6/math_6toSym33' error stop 'math_sym33to6/math_6toSym33'
call random_number(t66) call random_number(t66)
if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66))) & if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
error stop 'math_sym3333to66/math_66toSym3333' error stop 'math_sym3333to66/math_66toSym3333'
call random_number(v6) call random_number(v6)
@ -1242,6 +1278,9 @@ subroutine selfTest
if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
error stop 'math_det33/math_detSym33' error stop 'math_det33/math_detSym33'
if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) &
error stop 'math_mul3333xx33/math_identity4th'
if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) & if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) &
error stop 'math_inv33(math_I3)' error stop 'math_inv33(math_I3)'

178
src/mesh/DAMASK_mesh.f90
View File

@ -22,6 +22,15 @@ program DAMASK_mesh
implicit none implicit none
type :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1 !< frequency of result writes
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
@ -68,7 +77,7 @@ program DAMASK_mesh
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet, field, dimPlex PetscInt :: faceSet, currentFaceSet, dimPlex
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: & external :: &
@ -91,8 +100,7 @@ program DAMASK_mesh
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D) call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D)
CHKERRA(ierr) CHKERRA(ierr)
nActiveFields = 1 allocate(solres(1))
allocate(solres(nActiveFields))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
@ -103,8 +111,8 @@ program DAMASK_mesh
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_stringValue(line,chunkPos,i))
case('$loadcase') case('$Loadcase')
N_def = N_def + 1 N_def = N_def + 1
end select end select
enddo ! count all identifiers to allocate memory and do sanity check enddo ! count all identifiers to allocate memory and do sanity check
@ -114,32 +122,26 @@ program DAMASK_mesh
allocate(loadCases(N_def)) allocate(loadCases(N_def))
do i = 1, size(loadCases) do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields)) allocate(loadCases(i)%fieldBC(1))
field = 1 loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
enddo enddo
do i = 1, size(loadCases) do i = 1, size(loadCases)
do field = 1, nActiveFields loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements
select case (loadCases(i)%fieldBC(field)%ID) allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents))
case(FIELD_MECH_ID) do component = 1, loadCases(i)%fieldBC(1)%nComponents
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements select case (component)
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents)) case (1)
do component = 1, loadCases(i)%fieldBC(field)%nComponents loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID
select case (component) case (2)
case (1) loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID case (3)
case (2) loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
enddo
end select end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents enddo
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) do component = 1, loadCases(i)%fieldBC(1)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
enddo allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo enddo
enddo enddo
@ -151,52 +153,45 @@ program DAMASK_mesh
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_stringValue(line,chunkPos,i))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loadcase information ! loadcase information
case('$loadcase') case('$Loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1) currentLoadCase = IO_intValue(line,chunkPos,i+1)
case('face') case('Face')
currentFace = IO_intValue(line,chunkPos,i+1) currentFace = IO_intValue(line,chunkPos,i+1)
currentFaceSet = -1 currentFaceSet = -1
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo enddo
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time case('t')
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1) loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments case('N')
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) case('f_out')
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1) loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('guessreset','dropguessing') case('estimate_rate')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! boundary condition information ! boundary condition information
case('x','y','z') case('X','Y','Z')
select case(IO_lc(IO_stringValue(line,chunkPos,i))) select case(IO_stringValue(line,chunkPos,i))
case('x') case('X')
ID = COMPONENT_MECH_X_ID ID = COMPONENT_MECH_X_ID
case('y') case('Y')
ID = COMPONENT_MECH_Y_ID ID = COMPONENT_MECH_Y_ID
case('z') case('Z')
ID = COMPONENT_MECH_Z_ID ID = COMPONENT_MECH_Z_ID
end select end select
do field = 1, nActiveFields do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = &
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then .true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
.true. IO_floatValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif endif
enddo enddo
end select end select
@ -212,21 +207,16 @@ program DAMASK_mesh
print'(a,i0)', ' load case: ', currentLoadCase print'(a,i0)', ' load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
print'(a)', ' drop guessing along trajectory' print'(a)', ' drop guessing along trajectory'
do field = 1, nActiveFields print'(a)', ' Field '//trim(FIELD_MECH_label)
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
print'(a)', ' Field '//trim(FIELD_MECH_label)
end select do faceSet = 1, mesh_Nboundaries
do faceSet = 1, mesh_Nboundaries do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) &
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), &
print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & ' Component ', component, &
' Component ', component, & ' Value ', loadCases(currentLoadCase)%fieldBC(1)% &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% & componentBC(component)%Value(faceSet)
componentBC(component)%Value(faceSet) enddo
enddo
enddo
enddo enddo
print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count
@ -240,12 +230,7 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers ! doing initialization depending on active solvers
call FEM_Utilities_init call FEM_Utilities_init
do field = 1, nActiveFields call FEM_mechanical_init(loadCases(1)%fieldBC(1))
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mechanical_init(loadCases(1)%fieldBC(field))
end select
enddo
if (worldrank == 0) then if (worldrank == 0) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
@ -266,32 +251,14 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step stepFraction = stepFraction + 1 ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
@ -306,33 +273,16 @@ program DAMASK_mesh
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
flush(IO_STDOUT) flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mechanical_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve fields ! solve fields
stagIter = 0 stagIter = 0
stagIterate = .true. stagIterate = .true.
do while (stagIterate) do while (stagIterate)
do field = 1, nActiveFields solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) if(.not. solres(1)%converged) exit
case(FIELD_MECH_ID)
solres(field) = FEM_mechanical_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1 stagIter = stagIter + 1
stagIterate = stagIter < stagItMax & stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) & .and. all(solres(:)%converged) &

340
src/mesh/FEM_quadrature.f90
View File

@ -9,7 +9,7 @@ module FEM_quadrature
private private
integer, parameter :: & integer, parameter :: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary) maxOrder = 5 !< maximum integration order
real(pReal), dimension(2,3), parameter :: & real(pReal), dimension(2,3), parameter :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, & triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, & 1.0_pReal, -1.0_pReal, &
@ -20,8 +20,12 @@ module FEM_quadrature
-1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
type :: group_float !< variable length datatype
real(pReal), dimension(:), allocatable :: p
end type group_float
integer, dimension(2:3,maxOrder), public, protected :: & integer, dimension(2:3,maxOrder), public, protected :: &
FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder) FEM_nQuadrature !< number of quadrature points for spatial dimension(2-3) and interpolation order (1-maxOrder)
type(group_float), dimension(2:3,maxOrder), public, protected :: & type(group_float), dimension(2:3,maxOrder), public, protected :: &
FEM_quadrature_weights, & !< quadrature weights for each quadrature rule FEM_quadrature_weights, & !< quadrature weights for each quadrature rule
FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule
@ -35,145 +39,146 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief initializes FEM interpolation data !> @brief initializes FEM interpolation data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_quadrature_init subroutine FEM_quadrature_init()
print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6)
print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
print*, 'https://www.jstor.org/stable/43693493'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear
FEM_nQuadrature(2,1) = 1 FEM_nQuadrature(2,1) = 1
allocate(FEM_quadrature_weights(2,1)%p(1)) allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1)))
FEM_quadrature_weights(2,1)%p(1) = 1.0_pReal FEM_quadrature_weights(2,1)%p(1) = 1._pReal
allocate(FEM_quadrature_points (2,1)%p(2)) FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal])
FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1.0_pReal/3.0_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D quadratic ! 2D quadratic
FEM_nQuadrature(2,2) = 3 FEM_nQuadrature(2,2) = 3
allocate(FEM_quadrature_weights(2,2)%p(3)) allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2)))
FEM_quadrature_weights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal
allocate(FEM_quadrature_points (2,2)%p(6)) FEM_quadrature_points (2,2)%p = permutationStar21([1._pReal/6._pReal])
FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1.0_pReal/6.0_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D cubic ! 2D cubic
FEM_nQuadrature(2,3) = 6 FEM_nQuadrature(2,3) = 6
allocate(FEM_quadrature_weights(2,3)%p(6)) allocate(FEM_quadrature_weights(2,3)%p(FEM_nQuadrature(2,3)))
FEM_quadrature_weights(2,3)%p(1:3) = 0.22338158967801146570_pReal FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal
FEM_quadrature_weights(2,3)%p(4:6) = 0.10995174365532186764_pReal FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal
allocate(FEM_quadrature_points (2,3)%p(12)) FEM_quadrature_points (2,3)%p = [ &
FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.44594849091596488632_pReal]) permutationStar21([4.4594849091596489e-1_pReal]), &
FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.091576213509770743460_pReal]) permutationStar21([9.157621350977074e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D quartic ! 2D quartic
FEM_nQuadrature(2,4) = 12 FEM_nQuadrature(2,4) = 12
allocate(FEM_quadrature_weights(2,4)%p(12)) allocate(FEM_quadrature_weights(2,4)%p(FEM_nQuadrature(2,4)))
FEM_quadrature_weights(2,4)%p(1:3) = 0.11678627572638_pReal FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal
FEM_quadrature_weights(2,4)%p(4:6) = 0.05084490637021_pReal FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal
FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837_pReal FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal
allocate(FEM_quadrature_points (2,4)%p(24)) FEM_quadrature_points (2,4)%p = [ &
FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.24928674517091_pReal]) permutationStar21([2.4928674517091042e-1_pReal]), &
FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150_pReal]) permutationStar21([6.308901449150223e-2_pReal]), &
FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal]) permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D quintic ! 2D quintic
FEM_nQuadrature(2,5) = 16 FEM_nQuadrature(2,5) = 16
allocate(FEM_quadrature_weights(2,5)%p(16)) allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5)))
FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal FEM_quadrature_weights(2,5)%p(1:1) = 1.4431560767778717e-1_pReal
FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal
FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal
FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762320_pReal FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal
FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal FEM_quadrature_weights(2,5)%p(11:16) = 2.7230314174434994e-2_pReal
allocate(FEM_quadrature_points (2,5)%p(32)) FEM_quadrature_points (2,5)%p = [ &
FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal]) permutationStar3([1._pReal/3._pReal]), &
FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal]) permutationStar21([4.5929258829272316e-1_pReal]), &
FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal]) permutationStar21([1.705693077517602e-1_pReal]), &
FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.05054722831703_pReal]) permutationStar21([5.0547228317030975e-2_pReal]), &
FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal]) permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D linear ! 3D linear
FEM_nQuadrature(3,1) = 1 FEM_nQuadrature(3,1) = 1
allocate(FEM_quadrature_weights(3,1)%p(1)) allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1)))
FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal
allocate(FEM_quadrature_points (3,1)%p(3)) FEM_quadrature_points (3,1)%p = permutationStar4([0.25_pReal])
FEM_quadrature_points (3,1)%p(1:3)= permutationStar4([0.25_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D quadratic ! 3D quadratic
FEM_nQuadrature(3,2) = 4 FEM_nQuadrature(3,2) = 4
allocate(FEM_quadrature_weights(3,2)%p(4)) allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2)))
FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal
allocate(FEM_quadrature_points (3,2)%p(12)) FEM_quadrature_points (3,2)%p = permutationStar31([1.3819660112501052e-1_pReal])
FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.13819660112501051518_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D cubic ! 3D cubic
FEM_nQuadrature(3,3) = 14 FEM_nQuadrature(3,3) = 14
allocate(FEM_quadrature_weights(3,3)%p(14)) allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3)))
FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801585080_pReal FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal
FEM_quadrature_weights(3,3)%p(1:4) = 0.073493043116361949544_pReal FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal
FEM_quadrature_weights(3,3)%p(9:14) = 0.042546020777081466438_pReal FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal
allocate(FEM_quadrature_points (3,3)%p(42)) FEM_quadrature_points (3,3)%p = [ &
FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.092735250310891226402_pReal]) permutationStar31([9.273525031089123e-2_pReal]), &
FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.31088591926330060980_pReal]) permutationStar31([3.108859192633006e-1_pReal]), &
FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.045503704125649649492_pReal]) permutationStar22([4.5503704125649649e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D quartic ! 3D quartic (lower precision/unknown source)
FEM_nQuadrature(3,4) = 35 FEM_nQuadrature(3,4) = 35
allocate(FEM_quadrature_weights(3,4)%p(35)) allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4)))
FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal
FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal
FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal
FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal
FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal
allocate(FEM_quadrature_points (3,4)%p(105)) FEM_quadrature_points (3,4)%p = [ &
FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) permutationStar31([0.0267367755543735_pReal]), &
FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]), &
FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) permutationStar22([0.4547545999844830_pReal]), &
FEM_quadrature_points (3,4)%p(67:102) = permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]) permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]), &
FEM_quadrature_points (3,4)%p(103:105)= permutationStar4([0.25_pReal]) permutationStar4([0.25_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D quintic ! 3D quintic (lower precision/unknown source)
FEM_nQuadrature(3,5) = 56 FEM_nQuadrature(3,5) = 56
allocate(FEM_quadrature_weights(3,5)%p(56)) allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5)))
FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal
FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal
FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal
FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal
FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal
FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal
allocate(FEM_quadrature_points (3,5)%p(168)) FEM_quadrature_points (3,5)%p = [ &
FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) permutationStar31([0.0149520651530592_pReal]), &
FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]), &
FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]), &
FEM_quadrature_points (3,5)%p(85:120) = permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]) permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]), &
FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]) permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]), &
FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal]) permutationStar31([0.1344783347929940_pReal]) ]
call selfTest
end subroutine FEM_quadrature_init end subroutine FEM_quadrature_init
@ -186,11 +191,9 @@ pure function permutationStar3(point) result(qPt)
real(pReal), dimension(2) :: qPt real(pReal), dimension(2) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,1) :: temp
temp(:,1) = [point(1), point(1), point(1)] qPt = pack(matmul(triangle,reshape([ &
point(1), point(1), point(1)],[3,1])),.true.)
qPt = reshape(matmul(triangle, temp),[2])
end function permutationStar3 end function permutationStar3
@ -203,13 +206,11 @@ pure function permutationStar21(point) result(qPt)
real(pReal), dimension(6) :: qPt real(pReal), dimension(6) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,3) :: temp
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)] qPt = pack(matmul(triangle,reshape([ &
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)] point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), &
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)] point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), &
1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],[3,3])),.true.)
qPt = reshape(matmul(triangle, temp),[6])
end function permutationStar21 end function permutationStar21
@ -222,16 +223,14 @@ pure function permutationStar111(point) result(qPt)
real(pReal), dimension(12) :: qPt real(pReal), dimension(12) :: qPt
real(pReal), dimension(2), intent(in) :: point real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(3,6) :: temp
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)] qPt = pack(matmul(triangle,reshape([ &
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)] point(1), point(2), 1.0_pReal - point(1) - point(2), &
temp(:,3) = [point(2), point(1), 1.0_pReal - point(1) - point(2)] point(1), 1.0_pReal - point(1) - point(2), point(2), &
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)] point(2), point(1), 1.0_pReal - point(1) - point(2), &
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)] point(2), 1.0_pReal - point(1) - point(2), point(1), &
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)] 1.0_pReal - point(1) - point(2), point(2), point(1), &
1.0_pReal - point(1) - point(2), point(1), point(2)],[3,6])),.true.)
qPt = reshape(matmul(triangle, temp),[12])
end function permutationStar111 end function permutationStar111
@ -244,11 +243,9 @@ pure function permutationStar4(point) result(qPt)
real(pReal), dimension(3) :: qPt real(pReal), dimension(3) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,1) :: temp
temp(:,1) = [point(1), point(1), point(1), point(1)] qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(1), point(1), point(1)],[4,1])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[3])
end function permutationStar4 end function permutationStar4
@ -261,14 +258,12 @@ pure function permutationStar31(point) result(qPt)
real(pReal), dimension(12) :: qPt real(pReal), dimension(12) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,4) :: temp
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)] point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), &
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)] point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), &
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)] point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), &
1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],[4,4])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[12])
end function permutationStar31 end function permutationStar31
@ -276,21 +271,19 @@ end function permutationStar31
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief star 22 permutation of input !> @brief star 22 permutation of input
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function permutationStar22(point) result(qPt) function permutationStar22(point) result(qPt)
real(pReal), dimension(18) :: qPt real(pReal), dimension(18) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,6) :: temp
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)] point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), &
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)] point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), &
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)] 0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), &
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)] 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), &
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)] 0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), &
point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],[4,6])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[18])
end function permutationStar22 end function permutationStar22
@ -303,22 +296,20 @@ pure function permutationStar211(point) result(qPt)
real(pReal), dimension(36) :: qPt real(pReal), dimension(36) :: qPt
real(pReal), dimension(2), intent(in) :: point real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(4,12) :: temp
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)] point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), &
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)] point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)] point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), &
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), &
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)] point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)] point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), &
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)] 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), &
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)] 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), &
1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],[4,12])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[36])
end function permutationStar211 end function permutationStar211
@ -331,35 +322,60 @@ pure function permutationStar1111(point) result(qPt)
real(pReal), dimension(72) :: qPt real(pReal), dimension(72) :: qPt
real(pReal), dimension(3), intent(in) :: point real(pReal), dimension(3), intent(in) :: point
real(pReal), dimension(4,24) :: temp
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)] point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)] point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)] point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)] point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)] point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), &
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)] point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), &
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)] point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)] point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)] point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)] point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)] point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), &
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)] point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), &
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)] point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)] point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)] point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)] point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)] point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), &
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)] point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), &
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)] 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3), &
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)] 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2), &
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)] 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), &
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)] 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), &
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)] 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), &
1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],[4,24])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[72])
end function permutationStar1111 end function permutationStar1111
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of quadrature weights and points.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
integer :: o, d, n
real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal]
do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1)
do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1)
if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) &
error stop 'quadrature weights'
enddo
enddo
do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1)
do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1)
n = size(FEM_quadrature_points(d,o)%p,1)/d
if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) &
error stop 'quadrature points'
enddo
enddo
end subroutine selfTest
end module FEM_quadrature end module FEM_quadrature

35
src/mesh/FEM_utilities.f90
View File

@ -23,14 +23,8 @@ module FEM_utilities
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill real(pReal), public, protected :: wgt !< weighting factor 1/Nelems
integer, public, parameter :: maxFields = 6
integer, public :: nActiveFields = 0
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -49,10 +43,6 @@ module FEM_utilities
COMPONENT_MECH_Z_ID COMPONENT_MECH_Z_ID
end enum end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical :: &
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
@ -63,27 +53,17 @@ module FEM_utilities
end type tSolutionState end type tSolutionState
type, public :: tComponentBC type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable, dimension(:) :: Value real(pReal), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask logical, allocatable, dimension(:) :: Mask
end type tComponentBC end type tComponentBC
type, public :: tFieldBC type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer :: nComponents = 0 integer :: nComponents = 0
type(tComponentBC), allocatable :: componentBC(:) type(tComponentBC), allocatable, dimension(:) :: componentBC
end type tFieldBC end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
public :: & public :: &
FEM_utilities_init, & FEM_utilities_init, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
@ -109,8 +89,9 @@ subroutine FEM_utilities_init
integer :: structOrder !< order of displacement shape functions integer :: structOrder !< order of displacement shape functions
character(len=*), parameter :: & character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor ' PETSCDEBUG = ' -snes_view -snes_monitor '
PetscErrorCode :: ierr PetscErrorCode :: ierr
logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution
print'(/,a)', ' <<<+- FEM_utilities init -+>>>' print'(/,a)', ' <<<+- FEM_utilities init -+>>>'

16
src/mesh/discretization_mesh.f90
View File

@ -40,6 +40,11 @@ module discretization_mesh
mesh_maxNips !< max number of IPs in any CP element mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!! !!!! BEGIN DEPRECATED !!!!!
DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!) mesh_node0 !< node x,y,z coordinates (initially!)
@ -50,11 +55,6 @@ module discretization_mesh
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
public :: & public :: &
discretization_mesh_init, & discretization_mesh_init, &
mesh_FEM_build_ipVolumes, & mesh_FEM_build_ipVolumes, &
@ -71,16 +71,14 @@ subroutine discretization_mesh_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer, allocatable, dimension(:) :: chunkPos
integer :: dimPlex, & integer :: dimPlex, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
j, l, & j, &
debug_element, debug_ip debug_element, debug_ip
PetscSF :: sf PetscSF :: sf
DM :: globalMesh DM :: globalMesh
PetscInt :: nFaceSets PetscInt :: nFaceSets
PetscInt, pointer, dimension(:) :: pFaceSets PetscInt, pointer, dimension(:) :: pFaceSets
character(len=pStringLen), dimension(:), allocatable :: fileContent
IS :: faceSetIS IS :: faceSetIS
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -88,7 +86,7 @@ subroutine discretization_mesh_init(restart)
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
integer :: integrationOrder !< order of quadrature rule required integer :: integrationOrder !< order of quadrature rule required
type(tvec) :: coords_node0 type(tvec) :: coords_node0
print'(/,a)', ' <<<+- discretization_mesh init -+>>>' print'(/,a)', ' <<<+- discretization_mesh init -+>>>'

54
src/mesh/mesh_mech_FEM.f90
View File

@ -109,7 +109,7 @@ subroutine FEM_mechanical_init(fieldBC)
character(len=*), parameter :: prefix = 'mechFE_' character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal), dimension(3,3) :: devNull
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
@ -258,6 +258,7 @@ subroutine FEM_mechanical_init(fieldBC)
call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
call utilities_constitutiveResponse(0.0_pReal,devNull,.true.)
end subroutine FEM_mechanical_init end subroutine FEM_mechanical_init
@ -288,8 +289,8 @@ type(tSolutionState) function FEM_mechanical_solution( &
params%timeinc = timeinc params%timeinc = timeinc
params%fieldBC = fieldBC params%fieldBC = fieldBC
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false. terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
@ -397,7 +398,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
ForwardData = .false. ForwardData = .false.
@ -670,7 +671,7 @@ end subroutine FEM_mechanical_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment) !> @brief Calculate current coordinates (both nodal and ip coordinates)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_updateCoords() subroutine FEM_mechanical_updateCoords()
@ -678,21 +679,35 @@ subroutine FEM_mechanical_updateCoords()
nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
nodeCoords !< nodal coordinates (3,Nnodes) nodeCoords !< nodal coordinates (3,Nnodes)
real(pReal), pointer, dimension(:,:,:) :: &
ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems)
integer :: &
qPt, &
comp, &
qOffset, &
nOffset
DM :: dm_local DM :: dm_local
Vec :: x_local Vec :: x_local
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscInt :: dimPlex, pStart, pEnd, p, s, e PetscInt :: pStart, pEnd, p, s, e, q, &
cellStart, cellEnd, c, n
PetscSection :: section PetscSection :: section
PetscQuadrature :: mechQuad
PetscReal, dimension(:), pointer :: basisField, basisFieldDer
PetscScalar, dimension(:), pointer :: x_scal
call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,mechQuad,ierr); CHKERRQ(ierr)
call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr) call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr)
! write cell vertex displacements
call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr) call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr)
allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal) allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
do p=pStart, pEnd-1 do p=pStart, pEnd-1
call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr) call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr)
nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e) nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
@ -700,6 +715,31 @@ subroutine FEM_mechanical_updateCoords()
call discretization_setNodeCoords(nodeCoords) call discretization_setNodeCoords(nodeCoords)
call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
! write ip displacements
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr)
allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal)
do c=cellStart,cellEnd-1
qOffset=0
call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element
do qPt=0,nQuadrature-1
qOffset= qPt * (size(basisField)/nQuadrature)
do comp=0,dimPlex-1 !< loop over components
nOffset=0
q = comp
do n=0,nBasis-1
ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+&
sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*&
x_scal(nOffset+1:nOffset+dimPlex))
q = q+dimPlex
nOffset = nOffset+dimPlex
enddo
enddo
enddo
call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr)
end do
call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature]))
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mechanical_updateCoords end subroutine FEM_mechanical_updateCoords

34
src/parallelization.f90
View File

@ -24,9 +24,18 @@ module parallelization
worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only)
worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only)
#ifdef PETSC #ifndef PETSC
public :: parallelization_bcast_str
contains
subroutine parallelization_bcast_str(string)
character(len=:), allocatable, intent(inout) :: string
end subroutine parallelization_bcast_str
#else
public :: & public :: &
parallelization_init parallelization_init, &
parallelization_bcast_str
contains contains
@ -101,6 +110,27 @@ subroutine parallelization_init
!$ call omp_set_num_threads(OMP_NUM_THREADS) !$ call omp_set_num_threads(OMP_NUM_THREADS)
end subroutine parallelization_init end subroutine parallelization_init
!--------------------------------------------------------------------------------------------------
!> @brief Broadcast a string from process 0.
!--------------------------------------------------------------------------------------------------
subroutine parallelization_bcast_str(string)
character(len=:), allocatable, intent(inout) :: string
integer :: strlen, ierr ! pI64 for strlen not supported by MPI
if (worldrank == 0) strlen = len(string)
call MPI_Bcast(strlen,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr)
if (worldrank /= 0) allocate(character(len=strlen)::string)
call MPI_Bcast(string,strlen,MPI_CHARACTER,0,MPI_COMM_WORLD, ierr)
end subroutine parallelization_bcast_str
#endif #endif
end module parallelization end module parallelization

30
src/phase_mechanical.f90
View File

@ -232,8 +232,8 @@ module subroutine mechanical_init(materials,phases)
allocate(phase_mechanical_F0(phases%length)) allocate(phase_mechanical_F0(phases%length))
allocate(phase_mechanical_Li(phases%length)) allocate(phase_mechanical_Li(phases%length))
allocate(phase_mechanical_Li0(phases%length)) allocate(phase_mechanical_Li0(phases%length))
allocate(phase_mechanical_Lp0(phases%length))
allocate(phase_mechanical_Lp(phases%length)) allocate(phase_mechanical_Lp(phases%length))
allocate(phase_mechanical_Lp0(phases%length))
allocate(phase_mechanical_S(phases%length)) allocate(phase_mechanical_S(phases%length))
allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_P(phases%length))
allocate(phase_mechanical_S0(phases%length)) allocate(phase_mechanical_S0(phases%length))
@ -241,20 +241,20 @@ module subroutine mechanical_init(materials,phases)
do ph = 1, phases%length do ph = 1, phases%length
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanical') mech => phase%get('mechanical')
@ -517,7 +517,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
enddo LpLoop enddo LpLoop
call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, ph,en) S, Fi_new, ph,en)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -1030,7 +1030,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%State0(:,en) allocate(subState0,source=plasticState(ph)%State0(:,en))
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@ -1153,12 +1153,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
en = material_phaseEntry(co,ce) en = material_phaseEntry(co,ce)
call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
phase_mechanical_Fe(ph)%data(1:3,1:3,en), & phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
phase_mechanical_S(ph)%data(1:3,1:3,en), & phase_mechanical_S(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en), & phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
ph,en) ph,en)
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en))

77
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -18,7 +18,7 @@ submodule(phase:plastic) dislotungsten
Q_cl = 1.0_pReal !< activation energy for dislocation climb Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
b_sl, & !< magnitude of Burgers vector [m] b_sl, & !< magnitude of Burgers vector [m]
D_a, & d_caron, & !< distance of spontaneous annhihilation
i_sl, & !< Adj. parameter for distance between 2 forest dislocations i_sl, & !< Adj. parameter for distance between 2 forest dislocations
f_at, & !< factor to calculate atomic volume f_at, & !< factor to calculate atomic volume
tau_Peierls, & !< Peierls stress tau_Peierls, & !< Peierls stress
@ -56,7 +56,7 @@ submodule(phase:plastic) dislotungsten
type :: tDisloTungstendependentState type :: tDisloTungstendependentState
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
Lambda_sl, & Lambda_sl, &
threshold_stress tau_pass
end type tDisloTungstendependentState end type tDisloTungstendependentState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -172,7 +172,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%D_0 = pl%get_asFloat('D_0') prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.) prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
@ -191,7 +190,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%f_at = math_expand(prm%f_at, N_sl) prm%f_at = math_expand(prm%f_at, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl) prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl
! sanity checks ! sanity checks
if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
@ -202,12 +201,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls'
if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)'
if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl'
else slipActive else slipActive
rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, & allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, &
prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
source = emptyRealArray) source = emptyRealArray)
allocate(prm%forestProjection(0,0)) allocate(prm%forestProjection(0,0))
@ -246,8 +246,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers), source=0.0_pReal)
end associate end associate
@ -316,8 +316,6 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
ph, & ph, &
en en
real(pReal) :: &
VacancyDiffusion
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos, dot_gamma_neg,& dot_gamma_pos, dot_gamma_neg,&
tau_pos,& tau_pos,&
@ -325,38 +323,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
v_cl, & v_cl, &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
dip_distance d_hat
associate(prm => param(ph), stt => state(ph),& associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
dot => dotState(ph), dst => dependentState(ph))
call kinetics(Mp,T,ph,en,& call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, & dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg) tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg where(dEq0(tau_pos)) ! ToDo: use avg of +/-
dot_rho_dip_formation = 0.0_pReal dot_rho_dip_formation = 0.0_pReal
dot_rho_dip_climb = 0.0_pReal dot_rho_dip_climb = 0.0_pReal
else where else where
dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/-
prm%D_a, & ! lower limit prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper limit dst%Lambda_sl(:,en)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
0.0_pReal, & 0.0_pReal, &
prm%dipoleformation) prm%dipoleformation)
v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) &
* (1.0_pReal/(dip_distance+prm%D_a)) * (1.0_pReal/(d_hat+prm%d_caron))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency?
end where end where
dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
- dot_rho_dip_formation & - dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 single edge dislocations - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges
dot%rho_dip(:,en) = dot_rho_dip_formation & dot%rho_dip(:,en) = dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole
- dot_rho_dip_climb - dot_rho_dip_climb
end associate end associate
@ -377,11 +373,11 @@ module subroutine dislotungsten_dependentState(ph,en)
dislocationSpacing dislocationSpacing
associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%threshold_stress(:,en) = prm%mu*prm%b_sl & dst%tau_pass(:,en) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
@ -416,7 +412,7 @@ module subroutine plastic_dislotungsten_results(ph,group)
if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), &
'mean free path for slip','m') 'mean free path for slip','m')
case('tau_pass') case('tau_pass')
if(prm%sum_N_sl>0) call results_writeDataset(dst%threshold_stress,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(dst%tau_pass,group,trim(prm%output(o)), &
'threshold stress for slip','Pa') 'threshold stress for slip','Pa')
end select end select
enddo outputsLoop enddo outputsLoop
@ -456,8 +452,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
StressRatio_p,StressRatio_pminus1, & StressRatio_p,StressRatio_pminus1, &
dvel, vel, & dvel, vel, &
tau_pos,tau_neg, & tau_pos,tau_neg, &
t_n, t_k, dtk,dtn, & t_n, t_k, dtk,dtn
needsGoodName ! ToDo: @Karo: any idea?
integer :: j integer :: j
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
@ -475,13 +470,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, &
effectiveLength => dst%Lambda_sl(:,en) - prm%w) effectiveLength => dst%Lambda_sl(:,en) - prm%w)
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%h/(t_n + t_k) vel = prm%h/(t_n + t_k)
@ -492,7 +486,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
end where significantPositiveTau end where significantPositiveTau
if (present(ddot_gamma_dtau_pos)) then if (present(ddot_gamma_dtau_pos)) then
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_pos dtk = -1.0_pReal * t_k / tau_pos
@ -505,13 +499,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
end where significantPositiveTau2 end where significantPositiveTau2
endif endif
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength)
t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
vel = prm%h/(t_n + t_k) vel = prm%h/(t_n + t_k)
@ -522,7 +515,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
end where significantNegativeTau end where significantNegativeTau
if (present(ddot_gamma_dtau_neg)) then if (present(ddot_gamma_dtau_neg)) then
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
* (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_neg dtk = -1.0_pReal * t_k / tau_neg

286
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -16,13 +16,11 @@ submodule(phase:plastic) dislotwin
real(pReal) :: & real(pReal) :: &
mu = 1.0_pReal, & !< equivalent shear modulus mu = 1.0_pReal, & !< equivalent shear modulus
nu = 1.0_pReal, & !< equivalent shear Poisson's ratio nu = 1.0_pReal, & !< equivalent shear Poisson's ratio
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal, & !< activation energy for dislocation climb Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size D = 1.0_pReal, & !< grain size
p_sb = 1.0_pReal, & !< p-exponent in shear band velocity p_sb = 1.0_pReal, & !< p-exponent in shear band velocity
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors
L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors
@ -42,7 +40,7 @@ submodule(phase:plastic) dislotwin
b_sl, & !< absolute length of Burgers vector [m] for each slip system b_sl, & !< absolute length of Burgers vector [m] for each slip system
b_tw, & !< absolute length of Burgers vector [m] for each twin system b_tw, & !< absolute length of Burgers vector [m] for each twin system
b_tr, & !< absolute length of Burgers vector [m] for each transformation system b_tr, & !< absolute length of Burgers vector [m] for each transformation system
Q_s,& !< activation energy for glide [J] for each slip system Q_sl,& !< activation energy for glide [J] for each slip system
v_0, & !< dislocation velocity prefactor [m/s] for each slip system v_0, & !< dislocation velocity prefactor [m/s] for each slip system
dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system
dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system
@ -55,7 +53,8 @@ submodule(phase:plastic) dislotwin
s, & !< s-exponent in trans nucleation rate s, & !< s-exponent in trans nucleation rate
tau_0, & !< strength due to elements in solid solution tau_0, & !< strength due to elements in solid solution
gamma_char, & !< characteristic shear for twins gamma_char, & !< characteristic shear for twins
B !< drag coefficient B, & !< drag coefficient
d_caron !< distance of spontaneous annhihilation
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
h_sl_sl, & !< components of slip-slip interaction matrix h_sl_sl, & !< components of slip-slip interaction matrix
h_sl_tw, & !< components of slip-twin interaction matrix h_sl_tw, & !< components of slip-twin interaction matrix
@ -206,7 +205,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl))
prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl))
prm%Q_s = pl%get_as1dFloat('Q_s', requiredSize=size(N_sl)) prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl))
prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl))
@ -214,9 +213,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), & prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal, i=1,size(N_sl))]) defaultVal=[(0.0_pReal, i=1,size(N_sl))])
prm%D_a = pl%get_asFloat('D_a') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.)
prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.)
@ -231,28 +228,29 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_dip_0 = math_expand(rho_dip_0, N_sl) rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl) prm%Q_sl = math_expand(prm%Q_sl, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, N_sl) prm%q = math_expand(prm%q, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl) prm%tau_0 = math_expand(prm%tau_0, N_sl)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl
! sanity checks ! sanity checks
if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0'
if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl'
if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0'
if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0'
if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0'
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%Q_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_sl'
if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)'
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl'
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl'
else slipActive else slipActive
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
endif slipActive endif slipActive
@ -516,7 +514,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
integer :: i,k,l,m,n integer :: i,k,l,m,n
real(pReal) :: & real(pReal) :: &
f_unrotated,StressRatio_p,& f_unrotated,StressRatio_p,&
BoltzmannRatio, & E_kB_T, &
ddot_gamma_dtau, & ddot_gamma_dtau, &
tau tau
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
@ -586,7 +584,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
shearBandingContribution: if(dNeq0(prm%v_sb)) then shearBandingContribution: if(dNeq0(prm%v_sb)) then
BoltzmannRatio = prm%E_sb/(kB*T) E_kB_T = prm%E_sb/(kB*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
do i = 1,6 do i = 1,6
@ -596,8 +594,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
significantShearBandStress: if (abs(tau) > tol_math_check) then significantShearBandStress: if (abs(tau) > tol_math_check) then
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
@ -631,7 +629,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
integer :: i integer :: i
real(pReal) :: & real(pReal) :: &
f_unrotated, & f_unrotated, &
rho_dip_distance, & d_hat, &
v_cl, & !< climb velocity v_cl, & !< climb velocity
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
@ -639,70 +637,67 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
rho_dip_distance_min, &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_tw dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph), dst => dependentState(ph))
f_unrotated = 1.0_pReal & associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en))
call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) f_unrotated = 1.0_pReal &
dot%gamma_sl(:,en) = abs(dot_gamma_sl) - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en))
rho_dip_distance_min = prm%D_a*prm%b_sl call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
dot%gamma_sl(:,en) = abs(dot_gamma_sl)
slipState: do i = 1, prm%sum_N_sl slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_formation(i) = 0.0_pReal
dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress
rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en))
rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i))
dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
* stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else else significantSlipStress
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en))
b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) & d_hat = math_clip(d_hat, left = prm%d_caron(i))
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
1.0_pReal, &
prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) &
/ (rho_dip_distance-rho_dip_distance_min(i)) * stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
endif
endif significantSlipStress
enddo slipState
dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & if (dEq(d_hat,prm%d_caron(i))) then
- dot_rho_dip_formation & dot_rho_dip_climb(i) = 0.0_pReal
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) else
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
1.0_pReal, &
prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot%rho_dip(:,en) = dot_rho_dip_formation & dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & / (d_hat-prm%d_caron(i))
- dot_rho_dip_climb endif
endif significantSlipStress
enddo slipState
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char - dot_rho_dip_formation &
- 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl)
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) dot%rho_dip(:,en) = dot_rho_dip_formation &
dot%f_tr(:,en) = f_unrotated*dot_gamma_tr - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
- dot_rho_dip_climb
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
end associate end associate
@ -763,19 +758,17 @@ module subroutine dislotwin_dependentState(T,ph,en)
dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
!* threshold stress for growing twin/martensite !* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) & dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_tw)
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
if(prm%sum_N_tr == prm%sum_N_sl) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_tr) &
dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & + prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
@ -867,7 +860,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
tau, & tau, &
stressRatio, & stressRatio, &
StressRatio_p, & StressRatio_p, &
BoltzmannRatio, & Q_kB_T, &
v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned)
v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned)
dV_wait_inverse_dTau, & dV_wait_inverse_dTau, &
@ -876,33 +869,34 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
tau_eff !< effective resolved stress tau_eff !< effective resolved stress
integer :: i integer :: i
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
tau_eff = abs(tau)-dst%tau_pass(:,en) tau_eff = abs(tau)-dst%tau_pass(:,en)
significantStress: where(tau_eff > tol_math_check) significantStress: where(tau_eff > tol_math_check)
stressRatio = tau_eff/prm%tau_0 stressRatio = tau_eff/prm%tau_0
StressRatio_p = stressRatio** prm%p StressRatio_p = stressRatio** prm%p
BoltzmannRatio = prm%Q_s/(kB*T) Q_kB_T = prm%Q_sl/(kB*T)
v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q)
v_run_inverse = prm%B/(tau_eff*prm%b_sl) v_run_inverse = prm%B/(tau_eff*prm%b_sl)
dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * Q_kB_T &
* (stressRatio**(prm%p-1.0_pReal)) & * (stressRatio**(prm%p-1.0_pReal)) &
* (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) &
/ prm%tau_0 / prm%tau_0
dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff
dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
/ (v_wait_inverse+v_run_inverse)**2.0_pReal / (v_wait_inverse+v_run_inverse)**2.0_pReal
ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl
else where significantStress else where significantStress
dot_gamma_sl = 0.0_pReal dot_gamma_sl = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate
@ -945,34 +939,35 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
integer :: i,s1,s2 integer :: i,s1,s2
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_tw do i = 1, prm%sum_N_tw
tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
isFCC: if (prm%fccTwinTransNucleation) then isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i) s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i) s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
(prm%L_tw*prm%b_sl(i))*& (prm%L_tw*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) ! P_ncs (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))
else else
Ndot0=0.0_pReal Ndot0=0.0_pReal
end if end if
else isFCC else isFCC
Ndot0=prm%dot_N_0_tw(i) Ndot0=prm%dot_N_0_tw(i)
endif isFCC endif isFCC
enddo enddo
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r
dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r)
ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
else where significantStress else where significantStress
dot_gamma_tw = 0.0_pReal dot_gamma_tw = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate
@ -1011,36 +1006,37 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
Ndot0, & Ndot0, &
stressRatio_s, & stressRatio_s, &
ddot_gamma_dtau ddot_gamma_dtau
integer :: i,s1,s2 integer :: i,s1,s2
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_tr do i = 1, prm%sum_N_tr
tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
isFCC: if (prm%fccTwinTransNucleation) then isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i) s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i) s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
(prm%L_tr*prm%b_sl(i))*& (prm%L_tr*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) ! P_ncs (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))
else else
Ndot0=0.0_pReal Ndot0=0.0_pReal
end if end if
else isFCC else isFCC
Ndot0=prm%dot_N_0_tr(i) Ndot0=prm%dot_N_0_tr(i)
endif isFCC endif isFCC
enddo enddo
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s
dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s)
ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
else where significantStress else where significantStress
dot_gamma_tr = 0.0_pReal dot_gamma_tr = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate

4
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -359,11 +359,11 @@ module subroutine plastic_kinehardening_results(ph,group)
case ('xi') case ('xi')
if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
'resistance against plastic slip','Pa') 'resistance against plastic slip','Pa')
case ('tau_b') !ToDo: chi case ('chi')
if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), &
'back stress','Pa') 'back stress','Pa')
case ('sgn(gamma)') case ('sgn(gamma)')
if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int if(prm%sum_N_sl>0) call results_writeDataset(int(stt%sgn_gamma),group,trim(prm%output(o)), &
'sense of shear','1') 'sense of shear','1')
case ('chi_0') case ('chi_0')
if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), &

6
src/prec.f90
View File

@ -27,11 +27,6 @@ module prec
real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
type :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p
end type group_float
type :: tState type :: tState
integer :: & integer :: &
sizeState = 0, & !< size of state sizeState = 0, & !< size of state
@ -94,6 +89,7 @@ end subroutine prec_init
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ ! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! AlmostEqualRelative ! AlmostEqualRelative
! ToDo: Use 'spacing': https://gcc.gnu.org/onlinedocs/gfortran/SPACING.html#SPACING
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical elemental pure function dEq(a,b,tol) logical elemental pure function dEq(a,b,tol)

5
src/quit.f90
View File

@ -6,7 +6,10 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine quit(stop_id) subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
use PetscSys use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI_f08
#endif
use HDF5 use HDF5
implicit none implicit none

50
src/results.f90
View File

@ -52,6 +52,7 @@ module results
results_openGroup, & results_openGroup, &
results_closeGroup, & results_closeGroup, &
results_writeDataset, & results_writeDataset, &
results_writeDataset_str, &
results_setLink, & results_setLink, &
results_addAttribute, & results_addAttribute, &
results_removeLink, & results_removeLink, &
@ -64,16 +65,20 @@ subroutine results_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
character(len=pPathLen) :: commandLine character(len=pPathLen) :: commandLine
integer :: hdferr
integer(HID_T) :: group_id
character(len=:), allocatable :: date
print'(/,a)', ' <<<+- results init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- results init -+>>>'; flush(IO_STDOUT)
print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017' print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL
if(.not. restart) then if (.not. restart) then
resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w') resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w')
call results_addAttribute('DADF5_version_major',0) call results_addAttribute('DADF5_version_major',0)
call results_addAttribute('DADF5_version_minor',13) call results_addAttribute('DADF5_version_minor',14)
call get_command_argument(0,commandLine) call get_command_argument(0,commandLine)
call results_addAttribute('creator',trim(commandLine)//' '//DAMASKVERSION) call results_addAttribute('creator',trim(commandLine)//' '//DAMASKVERSION)
call results_addAttribute('created',now()) call results_addAttribute('created',now())
@ -81,18 +86,34 @@ subroutine results_init(restart)
call results_addAttribute('call',trim(commandLine)) call results_addAttribute('call',trim(commandLine))
call results_closeGroup(results_addGroup('cell_to')) call results_closeGroup(results_addGroup('cell_to'))
call results_addAttribute('description','mappings to place data in space','cell_to') call results_addAttribute('description','mappings to place data in space','cell_to')
call results_closeJobFile call results_closeGroup(results_addGroup('setup'))
call results_addAttribute('description','input data used to run the simulation','setup')
else
date = now()
call results_openJobFile
call get_command(commandLine)
call results_addAttribute('call (restart at '//date//')',trim(commandLine))
call h5gmove_f(resultsFile,'setup','tmp',hdferr)
call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp')
call results_closeGroup(results_addGroup('setup'))
call results_addAttribute('description','input data used to run the simulation','setup')
call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr)
endif endif
call results_closeJobFile
end subroutine results_init end subroutine results_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief opens the results file to append data !> @brief opens the results file to append data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine results_openJobFile subroutine results_openJobFile(parallel)
resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','a') logical, intent(in), optional :: parallel
resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','a',parallel)
end subroutine results_openJobFile end subroutine results_openJobFile
@ -297,6 +318,25 @@ subroutine results_removeLink(link)
end subroutine results_removeLink end subroutine results_removeLink
!--------------------------------------------------------------------------------------------------
!> @brief Store string dataset.
!> @details Not collective, must be called by one process at at time.
!--------------------------------------------------------------------------------------------------
subroutine results_writeDataset_str(dataset,group,label,description)
character(len=*), intent(in) :: label,group,description,dataset
integer(HID_T) :: groupHandle
groupHandle = results_openGroup(group)
call HDF5_write_str(dataset,groupHandle,label)
call executionStamp(group//'/'//label,description)
call HDF5_closeGroup(groupHandle)
end subroutine results_writeDataset_str
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Store real scalar dataset with associated metadata. !> @brief Store real scalar dataset with associated metadata.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------