diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 99daef5e5..2464ae7d7 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -158,12 +158,12 @@ Post_AverageDown: - master - release -#Post_General: -# stage: postprocessing -# script: PostProcessing/test.py -# except: -# - master -# - release +Post_General: + stage: postprocessing + script: PostProcessing/test.py + except: + - master + - release Post_GeometryReconstruction: stage: postprocessing @@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip: - master - release -#TextureComponents: -# stage: spectral -# script: TextureComponents/test.py -# except: -# - master -# - release +TextureComponents: + stage: spectral + script: TextureComponents/test.py + except: + - master + - release ################################################################################################### @@ -468,27 +468,24 @@ AbaqusStd: script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT abaqus - except: - - master - - release + only: + - development Marc: stage: createDocumentation script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT marc - except: - - master - - release + only: + - development Spectral: stage: createDocumentation script: - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT spectral - except: - - master - - release + only: + - development ################################################################################################## backupData: diff --git a/LICENSE b/LICENSE index 630dc3a84..1ab20178c 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright 2011-18 Max-Planck-Institut für Eisenforschung GmbH +Copyright 2011-19 Max-Planck-Institut für Eisenforschung GmbH DAMASK is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by diff --git a/Makefile b/Makefile index a8a7a6e0f..cd2690cc7 100644 --- a/Makefile +++ b/Makefile @@ -7,11 +7,11 @@ all: spectral FEM processing .PHONY: spectral spectral: build/spectral - @(cd build/spectral;make --no-print-directory -ws all install;) + @(cd build/spectral;make -j4 --no-print-directory -ws all install;) .PHONY: FEM FEM: build/FEM - @(cd build/FEM; make --no-print-directory -ws all install;) + @(cd build/FEM; make -j4 --no-print-directory -ws all install;) .PHONY: build/spectral build/spectral: diff --git a/PRIVATE b/PRIVATE index b9a52a85c..59b0cbe89 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit b9a52a85cd65cc27a8e863302bd984abdcad1455 +Subproject commit 59b0cbe899f272476fb6f00f0f8860428e6ceba3 diff --git a/VERSION b/VERSION index 4fd318a98..2ad825361 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1187-gcd8ee450 +v2.0.2-1236-g1ef82e35 diff --git a/env/DAMASK.csh b/env/DAMASK.csh index 6b6c58d9d..1819dd305 100644 --- a/env/DAMASK.csh +++ b/env/DAMASK.csh @@ -64,7 +64,7 @@ endif setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS if ( ! $?PYTHONPATH ) then - setenv PYTHONPATH $DAMASK_ROOT/lib + setenv PYTHONPATH $DAMASK_ROOT/python else - setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH + setenv PYTHONPATH $DAMASK_ROOT/python:$PYTHONPATH endif diff --git a/env/DAMASK.sh b/env/DAMASK.sh index bd26a3ebb..fa2c8db25 100644 --- a/env/DAMASK.sh +++ b/env/DAMASK.sh @@ -95,7 +95,7 @@ if [ ! -z "$PS1" ]; then fi export DAMASK_NUM_THREADS -export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH +export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do unset "${var}" diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh index 4d5a1e47d..61b9c89f9 100644 --- a/env/DAMASK.zsh +++ b/env/DAMASK.zsh @@ -88,7 +88,7 @@ if [ ! -z "$PS1" ]; then fi export DAMASK_NUM_THREADS -export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH +export PYTHONPATH=$DAMASK_ROOT/python:$PYTHONPATH for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do unset "${var}" diff --git a/processing/post/addEuclideanDistance.py b/processing/post/addEuclideanDistance.py index 2e5794e1b..f759b7a8f 100755 --- a/processing/post/addEuclideanDistance.py +++ b/processing/post/addEuclideanDistance.py @@ -149,7 +149,6 @@ for name in filenames: errors = [] remarks = [] - column = {} if not 3 >= table.label_dimension(options.pos) >= 1: errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos)) diff --git a/processing/post/addMises.py b/processing/post/addMises.py index 55cf6552e..7e757ed9d 100755 --- a/processing/post/addMises.py +++ b/processing/post/addMises.py @@ -4,6 +4,7 @@ import os,sys,math import numpy as np from optparse import OptionParser +from collections import OrderedDict import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] @@ -63,10 +64,10 @@ for name in filenames: # ------------------------------------------ sanity checks ---------------------------------------- - items = { - 'strain': {'dim': 9, 'shape': [3,3], 'labels':options.strain, 'active':[], 'column': []}, - 'stress': {'dim': 9, 'shape': [3,3], 'labels':options.stress, 'active':[], 'column': []}, - } + items = OrderedDict([ + ('strain', {'dim': 9, 'shape': [3,3], 'labels':options.strain, 'active':[], 'column': []}), + ('stress', {'dim': 9, 'shape': [3,3], 'labels':options.stress, 'active':[], 'column': []}) + ]) errors = [] remarks = [] diff --git a/processing/post/groupTable.py b/processing/post/groupTable.py index 67d07a7d1..f78566304 100755 --- a/processing/post/groupTable.py +++ b/processing/post/groupTable.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python2.7 +#!/usr/bin/env python3 # -*- coding: UTF-8 no BOM -*- import os,sys @@ -21,18 +21,19 @@ scriptID = ' '.join([scriptName,damask.version]) # -------------------------------------------------------------------- parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """ -Apply a user-specified function to condense all rows for which column 'label' has identical values into a single row. -Output table will contain as many rows as there are different (unique) values in the grouping column. +Apply a user-specified function to condense into a single row all those rows for which columns 'label' have identical values. +Output table will contain as many rows as there are different (unique) values in the grouping column(s). Periodic domain averaging of coordinate values is supported. Examples: For grain averaged values, replace all rows of particular 'texture' with a single row containing their average. -""", version = scriptID) +{name} --label texture --function np.average data.txt +""".format(name = scriptName), version = scriptID) parser.add_option('-l','--label', dest = 'label', - type = 'string', metavar = 'string', - help = 'column label for grouping rows') + action = 'extend', metavar = '', + help = 'column label(s) for grouping rows') parser.add_option('-f','--function', dest = 'function', type = 'string', metavar = 'string', @@ -40,7 +41,7 @@ parser.add_option('-f','--function', parser.add_option('-a','--all', dest = 'all', action = 'store_true', - help = 'apply mapping function also to grouping column') + help = 'apply mapping function also to grouping column(s)') group = OptionGroup(parser, "periodic averaging", "") @@ -57,6 +58,7 @@ parser.add_option_group(group) parser.set_defaults(function = 'np.average', all = False, + label = [], boundary = [0.0, 1.0]) (options,filenames) = parser.parse_args() @@ -71,7 +73,7 @@ try: except: mapFunction = None -if options.label is None: +if options.label is []: parser.error('no grouping column specified.') if not hasattr(mapFunction,'__call__'): parser.error('function "{}" is not callable.'.format(options.function)) @@ -89,13 +91,20 @@ for name in filenames: # ------------------------------------------ sanity checks --------------------------------------- + remarks = [] + errors = [] + table.head_read() - if table.label_dimension(options.label) != 1: - damask.util.croak('column {} is not of scalar dimension.'.format(options.label)) - table.close(dismiss = True) # close ASCIItable and remove empty file + grpColumns = table.label_index(options.label)[::-1] + grpColumns = grpColumns[np.where(grpColumns>=0)] + + if len(grpColumns) == 0: errors.append('no valid grouping column present.') + + if remarks != []: damask.util.croak(remarks) + if errors != []: + damask.util.croak(errors) + table.close(dismiss=True) continue - else: - grpColumn = table.label_index(options.label) # ------------------------------------------ assemble info --------------------------------------- @@ -108,10 +117,9 @@ for name in filenames: indexrange = table.label_indexrange(options.periodic) if options.periodic is not None else [] rows,cols = table.data.shape - table.data = table.data[np.lexsort([table.data[:,grpColumn]])] # sort data by grpColumn - - values,index = np.unique(table.data[:,grpColumn], return_index = True) # unique grpColumn values and their positions - index = np.append(index,rows) # add termination position + table.data = table.data[np.lexsort(table.data[:,grpColumns].T)] # sort data by grpColumn(s) + values,index = np.unique(table.data[:,grpColumns], axis=0, return_index=True) # unique grpColumn values and their positions + index = sorted(np.append(index,rows)) # add termination position grpTable = np.empty((len(values), cols)) # initialize output for i in range(len(values)): # iterate over groups (unique values in grpColumn) @@ -119,7 +127,7 @@ for name in filenames: grpTable[i,indexrange] = \ periodicAverage(table.data[index[i]:index[i+1],indexrange],options.boundary) # apply periodicAverage mapping function - if not options.all: grpTable[i,grpColumn] = table.data[index[i],grpColumn] # restore grouping column value + if not options.all: grpTable[i,grpColumns] = table.data[index[i],grpColumns] # restore grouping column value table.data = grpTable diff --git a/processing/post/postResults.py b/processing/post/postResults.py index a5a2669d7..c0885a967 100755 --- a/processing/post/postResults.py +++ b/processing/post/postResults.py @@ -831,7 +831,7 @@ if options.info: elementsOfNode = {} Nelems = stat['NumberOfElements'] for e in range(Nelems): - if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Nelems >= 50 and e%(Nelems//50) == 0: # report in 2% steps if possible and avoid modulo by zero damask.util.progressBar(iteration=e,total=Nelems,prefix='1/3: connecting elements') for n in map(p.node_sequence,p.element(e).items): if n not in elementsOfNode: @@ -856,7 +856,7 @@ damask.util.progressBar(iteration=1,total=1,prefix='1/3: connecting elements') if options.nodalScalar: Npoints = stat['NumberOfNodes'] for n in range(Npoints): - if options.verbose and Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Npoints >= 50 and e%(Npoints//50) == 0: # report in 2% steps if possible and avoid modulo by zero damask.util.progressBar(iteration=n,total=Npoints,prefix='2/3: scanning nodes ') myNodeID = p.node_id(n) myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z] @@ -893,7 +893,7 @@ if options.nodalScalar: else: Nelems = stat['NumberOfElements'] for e in range(Nelems): - if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and Nelems >= 50 and e%(Nelems//50) == 0: # report in 2% steps if possible and avoid modulo by zero damask.util.progressBar(iteration=e,total=Nelems,prefix='2/3: scanning elements ') myElemID = p.element_id(e) myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z], @@ -1034,7 +1034,7 @@ for incCount,position in enumerate(locations): # walk through locations Ngroups = len(groups) for j,group in enumerate(groups): f = incCount*Ngroups + j - if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero + if options.verbose and (Ngroups*Nincs) >= 50 and f%((Ngroups*Nincs)//50) == 0: # report in 2% steps if possible and avoid modulo by zero damask.util.progressBar(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ') N = 0 # group member counter for (e,n,i,g,n_local) in group[1:]: # loop over group members @@ -1087,7 +1087,7 @@ for incCount,position in enumerate(locations): # walk through locations ['Crystallite']*len(options.crystalliteResult) + ['Constitutive']*len(options.constitutiveResult) ): - outputIndex = (list(zip(*outputFormat[resultType]['outputs']))[0]).index(label) # find the position of this output in the outputFormat + outputIndex = (list(zip(*outputFormat[resultType]['outputs']))[0]).index(label) # find the position of this output in the outputFormat length = int(outputFormat[resultType]['outputs'][outputIndex][1]) thisHead = heading('_',[[component,''.join( label.split() )] for component in range(int(length>1),length+int(length>1))]) if assembleHeader: header += thisHead diff --git a/processing/post/viewTable.py b/processing/post/viewTable.py index c01bdcddb..309f229e1 100755 --- a/processing/post/viewTable.py +++ b/processing/post/viewTable.py @@ -61,7 +61,14 @@ for name in filenames: table = damask.ASCIItable(name = name, buffered = False, labeled = options.labeled, readonly = True) except: continue - damask.util.report(scriptName,name) + details = ', '.join( + (['header'] if options.table else []) + + (['info'] if options.head or options.info else []) + + (['labels'] if options.head or options.labels else []) + + (['data'] if options.data else []) + + [] + ) + damask.util.report(scriptName,name + ('' if details == '' else ' -- '+details)) # ------------------------------------------ output head --------------------------------------- diff --git a/lib/damask/.gitignore b/python/damask/.gitignore similarity index 100% rename from lib/damask/.gitignore rename to python/damask/.gitignore diff --git a/lib/damask/__init__.py b/python/damask/__init__.py similarity index 100% rename from lib/damask/__init__.py rename to python/damask/__init__.py diff --git a/lib/damask/asciitable.py b/python/damask/asciitable.py similarity index 100% rename from lib/damask/asciitable.py rename to python/damask/asciitable.py diff --git a/lib/damask/colormaps.py b/python/damask/colormaps.py similarity index 100% rename from lib/damask/colormaps.py rename to python/damask/colormaps.py diff --git a/lib/damask/config/__init__.py b/python/damask/config/__init__.py similarity index 100% rename from lib/damask/config/__init__.py rename to python/damask/config/__init__.py diff --git a/lib/damask/config/material.py b/python/damask/config/material.py similarity index 100% rename from lib/damask/config/material.py rename to python/damask/config/material.py diff --git a/lib/damask/environment.py b/python/damask/environment.py similarity index 100% rename from lib/damask/environment.py rename to python/damask/environment.py diff --git a/lib/damask/geometry/__init__.py b/python/damask/geometry/__init__.py similarity index 100% rename from lib/damask/geometry/__init__.py rename to python/damask/geometry/__init__.py diff --git a/lib/damask/geometry/geometry.py b/python/damask/geometry/geometry.py similarity index 100% rename from lib/damask/geometry/geometry.py rename to python/damask/geometry/geometry.py diff --git a/lib/damask/geometry/marc.py b/python/damask/geometry/marc.py similarity index 100% rename from lib/damask/geometry/marc.py rename to python/damask/geometry/marc.py diff --git a/lib/damask/geometry/spectral.py b/python/damask/geometry/spectral.py similarity index 100% rename from lib/damask/geometry/spectral.py rename to python/damask/geometry/spectral.py diff --git a/lib/damask/orientation.py b/python/damask/orientation.py similarity index 100% rename from lib/damask/orientation.py rename to python/damask/orientation.py diff --git a/lib/damask/result.py b/python/damask/result.py similarity index 100% rename from lib/damask/result.py rename to python/damask/result.py diff --git a/lib/damask/result/marc2vtk.py b/python/damask/result/marc2vtk.py similarity index 100% rename from lib/damask/result/marc2vtk.py rename to python/damask/result/marc2vtk.py diff --git a/lib/damask/solver/__init__.py b/python/damask/solver/__init__.py similarity index 100% rename from lib/damask/solver/__init__.py rename to python/damask/solver/__init__.py diff --git a/lib/damask/solver/abaqus.py b/python/damask/solver/abaqus.py similarity index 100% rename from lib/damask/solver/abaqus.py rename to python/damask/solver/abaqus.py diff --git a/lib/damask/solver/marc.py b/python/damask/solver/marc.py similarity index 100% rename from lib/damask/solver/marc.py rename to python/damask/solver/marc.py diff --git a/lib/damask/solver/solver.py b/python/damask/solver/solver.py similarity index 100% rename from lib/damask/solver/solver.py rename to python/damask/solver/solver.py diff --git a/lib/damask/solver/spectral.py b/python/damask/solver/spectral.py similarity index 100% rename from lib/damask/solver/spectral.py rename to python/damask/solver/spectral.py diff --git a/lib/damask/test/__init__.py b/python/damask/test/__init__.py similarity index 100% rename from lib/damask/test/__init__.py rename to python/damask/test/__init__.py diff --git a/lib/damask/test/test.py b/python/damask/test/test.py similarity index 100% rename from lib/damask/test/test.py rename to python/damask/test/test.py diff --git a/lib/damask/util.py b/python/damask/util.py similarity index 94% rename from lib/damask/util.py rename to python/damask/util.py index e7d4a16bc..a88080df5 100644 --- a/lib/damask/util.py +++ b/python/damask/util.py @@ -134,7 +134,7 @@ class extendableOption(Option): # Print iterations progress # from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a -def progressBar(iteration, total, prefix='', suffix='', decimals=1, bar_length=100): +def progressBar(iteration, total, prefix='', bar_length=50): """ Call in a loop to create terminal progress bar @@ -142,16 +142,29 @@ def progressBar(iteration, total, prefix='', suffix='', decimals=1, bar_length=1 iteration - Required : current iteration (Int) total - Required : total iterations (Int) prefix - Optional : prefix string (Str) - suffix - Optional : suffix string (Str) - decimals - Optional : positive number of decimals in percent complete (Int) bar_length - Optional : character length of bar (Int) """ - str_format = "{0:." + str(decimals) + "f}" - percents = str_format.format(100 * (iteration / float(total))) - filled_length = int(round(bar_length * iteration / float(total))) - bar = '█' * filled_length + '-' * (bar_length - filled_length) + fraction = iteration / float(total) + if not hasattr(progressBar, "last_fraction"): # first call to function + progressBar.start_time = time.time() + progressBar.last_fraction = -1.0 + remaining_time = ' n/a' + else: + if fraction <= progressBar.last_fraction or iteration == 0: # reset: called within a new loop + progressBar.start_time = time.time() + progressBar.last_fraction = -1.0 + remaining_time = ' n/a' + else: + progressBar.last_fraction = fraction + remainder = (total - iteration) * (time.time()-progressBar.start_time)/iteration + remaining_time = '{: 3d}:'.format(int( remainder//3600)) + \ + '{:02d}:'.format(int((remainder//60)%60)) + \ + '{:02d}' .format(int( remainder %60)) - sys.stderr.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)), + filled_length = int(round(bar_length * fraction)) + bar = '█' * filled_length + '░' * (bar_length - filled_length) + + sys.stderr.write('\r{} {} {}'.format(prefix, bar, remaining_time)), if iteration == total: sys.stderr.write('\n\n') sys.stderr.flush() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2e4462243..a84609087 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -90,9 +90,7 @@ list(APPEND OBJECTFILES $) add_library (KINEMATICS OBJECT "kinematics_cleavage_opening.f90" "kinematics_slipplane_opening.f90" - "kinematics_thermal_expansion.f90" - "kinematics_vacancy_strain.f90" - "kinematics_hydrogen_strain.f90") + "kinematics_thermal_expansion.f90") add_dependencies(KINEMATICS DAMASK_HELPERS) list(APPEND OBJECTFILES $) @@ -102,10 +100,7 @@ add_library (SOURCE OBJECT "source_damage_isoBrittle.f90" "source_damage_isoDuctile.f90" "source_damage_anisoBrittle.f90" - "source_damage_anisoDuctile.f90" - "source_vacancy_phenoplasticity.f90" - "source_vacancy_irradiation.f90" - "source_vacancy_thermalfluc.f90") + "source_damage_anisoDuctile.f90") add_dependencies(SOURCE DAMASK_HELPERS) list(APPEND OBJECTFILES $) @@ -124,25 +119,6 @@ add_library(HOMOGENIZATION OBJECT add_dependencies(HOMOGENIZATION CRYSTALLITE) list(APPEND OBJECTFILES $) -add_library(HYDROGENFLUX OBJECT - "hydrogenflux_isoconc.f90" - "hydrogenflux_cahnhilliard.f90") -add_dependencies(HYDROGENFLUX CRYSTALLITE) -list(APPEND OBJECTFILES $) - -add_library(POROSITY OBJECT - "porosity_none.f90" - "porosity_phasefield.f90") -add_dependencies(POROSITY CRYSTALLITE) -list(APPEND OBJECTFILES $) - -add_library(VACANCYFLUX OBJECT - "vacancyflux_isoconc.f90" - "vacancyflux_isochempot.f90" - "vacancyflux_cahnhilliard.f90") -add_dependencies(VACANCYFLUX CRYSTALLITE) -list(APPEND OBJECTFILES $) - add_library(DAMAGE OBJECT "damage_none.f90" "damage_local.f90" @@ -158,7 +134,7 @@ add_dependencies(THERMAL CRYSTALLITE) list(APPEND OBJECTFILES $) add_library(DAMASK_ENGINE OBJECT "homogenization.f90") -add_dependencies(DAMASK_ENGINE THERMAL DAMAGE VACANCYFLUX POROSITY HYDROGENFLUX HOMOGENIZATION) +add_dependencies(DAMASK_ENGINE THERMAL DAMAGE HOMOGENIZATION) list(APPEND OBJECTFILES $) add_library(DAMASK_CPFE OBJECT "CPFEM2.f90") diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 674a557b5..847688d57 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -304,8 +304,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt homogState, & thermalState, & damageState, & - vacancyfluxState, & - hydrogenfluxState, & phaseAt, phasememberAt, & material_phase, & phase_plasticity, & @@ -421,8 +419,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt homogState (homog)%state0 = homogState (homog)%state thermalState (homog)%state0 = thermalState (homog)%state damageState (homog)%state0 = damageState (homog)%state - vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state - hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state enddo diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 2aed858a7..29e1ac744 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -203,8 +203,6 @@ subroutine CPFEM_age() homogState, & thermalState, & damageState, & - vacancyfluxState, & - hydrogenfluxState, & material_phase, & phase_plasticity, & phase_Nsources @@ -268,8 +266,6 @@ if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & homogState (homog)%state0 = homogState (homog)%state thermalState (homog)%state0 = thermalState (homog)%state damageState (homog)%state0 = damageState (homog)%state - vacancyfluxState (homog)%state0 = vacancyfluxState (homog)%state - hydrogenfluxState(homog)%state0 = hydrogenfluxState(homog)%state enddo if (restartWrite) then diff --git a/src/IO.f90 b/src/IO.f90 index af59b11b9..193580fcc 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1251,6 +1251,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'negative number systems requested' case (145_pInt) msg = 'too many systems requested' + case (146_pInt) + msg = 'number of values does not match' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 36f0244ef..8d3e9c816 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -21,14 +21,9 @@ #include "source_damage_isoDuctile.f90" #include "source_damage_anisoBrittle.f90" #include "source_damage_anisoDuctile.f90" -#include "source_vacancy_phenoplasticity.f90" -#include "source_vacancy_irradiation.f90" -#include "source_vacancy_thermalfluc.f90" #include "kinematics_cleavage_opening.f90" #include "kinematics_slipplane_opening.f90" #include "kinematics_thermal_expansion.f90" -#include "kinematics_vacancy_strain.f90" -#include "kinematics_hydrogen_strain.f90" #include "plastic_none.f90" #include "plastic_isotropic.f90" #include "plastic_phenopowerlaw.f90" @@ -47,12 +42,5 @@ #include "damage_none.f90" #include "damage_local.f90" #include "damage_nonlocal.f90" -#include "vacancyflux_isoconc.f90" -#include "vacancyflux_isochempot.f90" -#include "vacancyflux_cahnhilliard.f90" -#include "porosity_none.f90" -#include "porosity_phasefield.f90" -#include "hydrogenflux_isoconc.f90" -#include "hydrogenflux_cahnhilliard.f90" #include "homogenization.f90" #include "CPFEM.f90" diff --git a/src/config.f90 b/src/config.f90 index 441dd953c..c7fd95b43 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -550,7 +550,7 @@ end function getString !> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all !! values from the last occurrence. If key is not found exits with error unless default is given. !-------------------------------------------------------------------------------------------------- -function getFloats(this,key,defaultVal,requiredShape) +function getFloats(this,key,defaultVal,requiredShape,requiredSize) use IO, only: & IO_error, & IO_stringValue, & @@ -561,7 +561,8 @@ function getFloats(this,key,defaultVal,requiredShape) class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key real(pReal), dimension(:), intent(in), optional :: defaultVal - integer(pInt), dimension(:), intent(in), optional :: requiredShape + integer(pInt), dimension(:), intent(in), optional :: requiredShape ! not useful (is always 1D array) + integer(pInt), intent(in), optional :: requiredSize type(tPartitionedStringList), pointer :: item integer(pInt) :: i logical :: found, & @@ -588,6 +589,9 @@ function getFloats(this,key,defaultVal,requiredShape) if (.not. found) then if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif endif + if (present(requiredSize)) then + if(requiredSize /= size(getFloats)) call IO_error(146,ext_msg=key) + endif end function getFloats @@ -597,7 +601,7 @@ end function getFloats !> @details for cumulative keys, "()", values from all occurrences are return. Otherwise only all !! values from the last occurrence. If key is not found exits with error unless default is given. !-------------------------------------------------------------------------------------------------- -function getInts(this,key,defaultVal,requiredShape) +function getInts(this,key,defaultVal,requiredShape,requiredSize) use IO, only: & IO_error, & IO_stringValue, & @@ -608,7 +612,8 @@ function getInts(this,key,defaultVal,requiredShape) class(tPartitionedStringList), target, intent(in) :: this character(len=*), intent(in) :: key integer(pInt), dimension(:), intent(in), optional :: defaultVal, & - requiredShape + requiredShape ! not useful (is always 1D array) + integer(pInt), intent(in), optional :: requiredSize type(tPartitionedStringList), pointer :: item integer(pInt) :: i logical :: found, & @@ -635,6 +640,9 @@ function getInts(this,key,defaultVal,requiredShape) if (.not. found) then if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif endif + if (present(requiredSize)) then + if(requiredSize /= size(getInts)) call IO_error(146,ext_msg=key) + endif end function getInts diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 8294047e7..b20d6e65d 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -88,14 +88,9 @@ subroutine constitutive_init() SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & SOURCE_damage_anisoDuctile_ID, & - SOURCE_vacancy_phenoplasticity_ID, & - SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID, & KINEMATICS_cleavage_opening_ID, & KINEMATICS_slipplane_opening_ID, & KINEMATICS_thermal_expansion_ID, & - KINEMATICS_vacancy_strain_ID, & - KINEMATICS_hydrogen_strain_ID, & ELASTICITY_HOOKE_label, & PLASTICITY_NONE_label, & PLASTICITY_ISOTROPIC_label, & @@ -110,9 +105,6 @@ subroutine constitutive_init() SOURCE_damage_isoDuctile_label, & SOURCE_damage_anisoBrittle_label, & SOURCE_damage_anisoDuctile_label, & - SOURCE_vacancy_phenoplasticity_label, & - SOURCE_vacancy_irradiation_label, & - SOURCE_vacancy_thermalfluc_label, & plasticState, & sourceState @@ -129,14 +121,9 @@ subroutine constitutive_init() use source_damage_isoDuctile use source_damage_anisoBrittle use source_damage_anisoDuctile - use source_vacancy_phenoplasticity - use source_vacancy_irradiation - use source_vacancy_thermalfluc use kinematics_cleavage_opening use kinematics_slipplane_opening use kinematics_thermal_expansion - use kinematics_vacancy_strain - use kinematics_hydrogen_strain implicit none integer(pInt), parameter :: FILEUNIT = 204_pInt @@ -179,9 +166,6 @@ subroutine constitutive_init() if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init(FILEUNIT) if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init(FILEUNIT) if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init(FILEUNIT) - if (any(phase_source == SOURCE_vacancy_phenoplasticity_ID)) call source_vacancy_phenoplasticity_init(FILEUNIT) - if (any(phase_source == SOURCE_vacancy_irradiation_ID)) call source_vacancy_irradiation_init(FILEUNIT) - if (any(phase_source == SOURCE_vacancy_thermalfluc_ID)) call source_vacancy_thermalfluc_init(FILEUNIT) !-------------------------------------------------------------------------------------------------- ! parse kinematic mechanisms from config file @@ -189,8 +173,6 @@ subroutine constitutive_init() if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT) if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT) - if (any(phase_kinematics == KINEMATICS_vacancy_strain_ID)) call kinematics_vacancy_strain_init(FILEUNIT) - if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID)) call kinematics_hydrogen_strain_init(FILEUNIT) close(FILEUNIT) call config_deallocate('material.config/phase') @@ -283,21 +265,6 @@ subroutine constitutive_init() outputName = SOURCE_damage_anisoDuctile_label thisOutput => source_damage_anisoDuctile_output thisSize => source_damage_anisoDuctile_sizePostResult - case (SOURCE_vacancy_phenoplasticity_ID) sourceType - ins = source_vacancy_phenoplasticity_instance(ph) - outputName = SOURCE_vacancy_phenoplasticity_label - thisOutput => source_vacancy_phenoplasticity_output - thisSize => source_vacancy_phenoplasticity_sizePostResult - case (SOURCE_vacancy_irradiation_ID) sourceType - ins = source_vacancy_irradiation_instance(ph) - outputName = SOURCE_vacancy_irradiation_label - thisOutput => source_vacancy_irradiation_output - thisSize => source_vacancy_irradiation_sizePostResult - case (SOURCE_vacancy_thermalfluc_ID) sourceType - ins = source_vacancy_thermalfluc_instance(ph) - outputName = SOURCE_vacancy_thermalfluc_label - thisOutput => source_vacancy_thermalfluc_output - thisSize => source_vacancy_thermalfluc_sizePostResult case default sourceType knownSource = .false. end select sourceType @@ -512,8 +479,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e dLp_dMp = 0.0_pReal case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) - dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = phasememberAt(ipc,ip,el) @@ -560,6 +528,7 @@ end subroutine constitutive_LpAndItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: MD: S is Mi? !-------------------------------------------------------------------------------------------------- subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el) use prec, only: & @@ -568,8 +537,12 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e math_I3, & math_inv33, & math_det33, & - math_mul33x33 + math_mul33x33, & + math_Mandel6to33 use material, only: & + phasememberAt, & + phase_plasticity, & + phase_plasticityInstance, & phase_plasticity, & material_phase, & phase_kinematics, & @@ -577,9 +550,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e PLASTICITY_isotropic_ID, & KINEMATICS_cleavage_opening_ID, & KINEMATICS_slipplane_opening_ID, & - KINEMATICS_thermal_expansion_ID, & - KINEMATICS_vacancy_strain_ID, & - KINEMATICS_hydrogen_strain_ID + KINEMATICS_thermal_expansion_ID use plastic_isotropic, only: & plastic_isotropic_LiAndItsTangent use kinematics_cleavage_opening, only: & @@ -588,10 +559,6 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e kinematics_slipplane_opening_LiAndItsTangent use kinematics_thermal_expansion, only: & kinematics_thermal_expansion_LiAndItsTangent - use kinematics_vacancy_strain, only: & - kinematics_vacancy_strain_LiAndItsTangent - use kinematics_hydrogen_strain, only: & - kinematics_hydrogen_strain_LiAndItsTangent implicit none integer(pInt), intent(in) :: & @@ -607,19 +574,18 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e real(pReal), intent(out), dimension(3,3,3,3) :: & dLi_dS, & !< derivative of Li with respect to S dLi_dFi + real(pReal), dimension(3,3) :: & - my_Li !< intermediate velocity gradient - real(pReal), dimension(3,3,3,3) :: & - my_dLi_dS - real(pReal), dimension(3,3) :: & + my_Li, & !< intermediate velocity gradient FiInv, & temp_33 + real(pReal), dimension(3,3,3,3) :: & + my_dLi_dS real(pReal) :: & detFi integer(pInt) :: & - k !< counter in kinematics loop - integer(pInt) :: & - i, j + k, i, j, & + instance, of Li = 0.0_pReal dLi_dS = 0.0_pReal @@ -627,7 +593,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_isotropic_ID) plasticityType - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_Mandel6to33(S6),instance,of) case default plasticityType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal @@ -644,10 +612,6 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) - case (KINEMATICS_vacancy_strain_ID) kinematicsType - call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) - case (KINEMATICS_hydrogen_strain_ID) kinematicsType - call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) case default kinematicsType my_Li = 0.0_pReal my_dLi_dS = 0.0_pReal @@ -684,15 +648,9 @@ pure function constitutive_initialFi(ipc, ip, el) phase_kinematics, & phase_Nkinematics, & material_phase, & - KINEMATICS_thermal_expansion_ID, & - KINEMATICS_vacancy_strain_ID, & - KINEMATICS_hydrogen_strain_ID + KINEMATICS_thermal_expansion_ID use kinematics_thermal_expansion, only: & kinematics_thermal_expansion_initialStrain - use kinematics_vacancy_strain, only: & - kinematics_vacancy_strain_initialStrain - use kinematics_hydrogen_strain, only: & - kinematics_hydrogen_strain_initialStrain implicit none integer(pInt), intent(in) :: & @@ -711,12 +669,6 @@ pure function constitutive_initialFi(ipc, ip, el) case (KINEMATICS_thermal_expansion_ID) kinematicsType constitutive_initialFi = & constitutive_initialFi + kinematics_thermal_expansion_initialStrain(ipc, ip, el) - case (KINEMATICS_vacancy_strain_ID) kinematicsType - constitutive_initialFi = & - constitutive_initialFi + kinematics_vacancy_strain_initialStrain(ipc, ip, el) - case (KINEMATICS_hydrogen_strain_ID) kinematicsType - constitutive_initialFi = & - constitutive_initialFi + kinematics_hydrogen_strain_initialStrain(ipc, ip, el) end select kinematicsType enddo KinematicsLoop @@ -771,10 +723,7 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip phase_stiffnessDegradation, & damage, & damageMapping, & - porosity, & - porosityMapping, & - STIFFNESS_DEGRADATION_damage_ID, & - STIFFNESS_DEGRADATION_porosity_ID + STIFFNESS_DEGRADATION_damage_ID implicit none integer(pInt), intent(in) :: & @@ -804,8 +753,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip degradationType: select case(phase_stiffnessDegradation(d,material_phase(ipc,ip,el))) case (STIFFNESS_DEGRADATION_damage_ID) degradationType C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2_pInt - case (STIFFNESS_DEGRADATION_porosity_ID) degradationType - C = C * porosity(ho)%p(porosityMapping(ho)%p(ip,el))**2_pInt end select degradationType enddo DegradationLoop @@ -916,7 +863,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (math_Mandel33to6(Mp),ipc,ip,el) + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_isotropic_dotState (Mp,instance,of) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = phasememberAt(ipc,ip,el) @@ -986,19 +935,13 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) material_phase, & PLASTICITY_KINEHARDENING_ID, & PLASTICITY_NONLOCAL_ID, & - SOURCE_damage_isoBrittle_ID, & - SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID + SOURCE_damage_isoBrittle_ID use plastic_kinehardening, only: & plastic_kinehardening_deltaState use plastic_nonlocal, only: & plastic_nonlocal_deltaState use source_damage_isoBrittle, only: & source_damage_isoBrittle_deltaState - use source_vacancy_irradiation, only: & - source_vacancy_irradiation_deltaState - use source_vacancy_thermalfluc, only: & - source_vacancy_thermalfluc_deltaState implicit none integer(pInt), intent(in) :: & @@ -1035,12 +978,6 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) - case (SOURCE_vacancy_irradiation_ID) sourceType - call source_vacancy_irradiation_deltaState(ipc, ip, el) - - case (SOURCE_vacancy_thermalfluc_ID) sourceType - call source_vacancy_thermalfluc_deltaState(ipc, ip, el) - end select sourceType enddo SourceLoop @@ -1140,8 +1077,10 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) constitutive_postResults(startPos:endPos) = & - plastic_isotropic_postResults(S6,ipc,ip,el) + plastic_isotropic_postResults(Mp,instance,of) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType of = phasememberAt(ipc,ip,el) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 82a97dc53..4663caa9d 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -25,10 +25,7 @@ module homogenization materialpoint_sizeResults, & homogenization_maxSizePostResults, & thermal_maxSizePostResults, & - damage_maxSizePostResults, & - vacancyflux_maxSizePostResults, & - porosity_maxSizePostResults, & - hydrogenflux_maxSizePostResults + damage_maxSizePostResults real(pReal), dimension(:,:,:,:), allocatable, private :: & materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment @@ -100,13 +97,6 @@ subroutine homogenization_init use damage_none use damage_local use damage_nonlocal - use vacancyflux_isoconc - use vacancyflux_isochempot - use vacancyflux_cahnhilliard - use porosity_none - use porosity_phasefield - use hydrogenflux_isoconc - use hydrogenflux_cahnhilliard use IO use numerics, only: & worldrank @@ -155,33 +145,6 @@ subroutine homogenization_init if (any(damage_type == DAMAGE_nonlocal_ID)) & call damage_nonlocal_init(FILEUNIT) -!-------------------------------------------------------------------------------------------------- -! parse vacancy transport from config file - call IO_checkAndRewind(FILEUNIT) - if (any(vacancyflux_type == VACANCYFLUX_isoconc_ID)) & - call vacancyflux_isoconc_init() - if (any(vacancyflux_type == VACANCYFLUX_isochempot_ID)) & - call vacancyflux_isochempot_init(FILEUNIT) - if (any(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID)) & - call vacancyflux_cahnhilliard_init(FILEUNIT) - -!-------------------------------------------------------------------------------------------------- -! parse porosity from config file - call IO_checkAndRewind(FILEUNIT) - if (any(porosity_type == POROSITY_none_ID)) & - call porosity_none_init() - if (any(porosity_type == POROSITY_phasefield_ID)) & - call porosity_phasefield_init(FILEUNIT) - -!-------------------------------------------------------------------------------------------------- -! parse hydrogen transport from config file - call IO_checkAndRewind(FILEUNIT) - if (any(hydrogenflux_type == HYDROGENFLUX_isoconc_ID)) & - call hydrogenflux_isoconc_init() - if (any(hydrogenflux_type == HYDROGENFLUX_cahnhilliard_ID)) & - call hydrogenflux_cahnhilliard_init(FILEUNIT) - close(FILEUNIT) - !-------------------------------------------------------------------------------------------------- ! write description file for homogenization output mainProcess2: if (worldrank == 0) then @@ -277,83 +240,6 @@ subroutine homogenization_init enddo endif endif - i = vacancyflux_typeInstance(p) ! which instance of this vacancy flux type - valid = .true. ! assume valid - select case(vacancyflux_type(p)) ! split per vacancy flux type - case (VACANCYFLUX_isoconc_ID) - outputName = VACANCYFLUX_isoconc_label - thisNoutput => null() - thisOutput => null() - thisSize => null() - case (VACANCYFLUX_isochempot_ID) - outputName = VACANCYFLUX_isochempot_label - thisNoutput => vacancyflux_isochempot_Noutput - thisOutput => vacancyflux_isochempot_output - thisSize => vacancyflux_isochempot_sizePostResult - case (VACANCYFLUX_cahnhilliard_ID) - outputName = VACANCYFLUX_cahnhilliard_label - thisNoutput => vacancyflux_cahnhilliard_Noutput - thisOutput => vacancyflux_cahnhilliard_output - thisSize => vacancyflux_cahnhilliard_sizePostResult - case default - valid = .false. - end select - if (valid) then - write(FILEUNIT,'(a)') '(vacancyflux)'//char(9)//trim(outputName) - if (vacancyflux_type(p) /= VACANCYFLUX_isoconc_ID) then - do e = 1,thisNoutput(i) - write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - endif - i = porosity_typeInstance(p) ! which instance of this porosity type - valid = .true. ! assume valid - select case(porosity_type(p)) ! split per porosity type - case (POROSITY_none_ID) - outputName = POROSITY_none_label - thisNoutput => null() - thisOutput => null() - thisSize => null() - case (POROSITY_phasefield_ID) - outputName = POROSITY_phasefield_label - thisNoutput => porosity_phasefield_Noutput - thisOutput => porosity_phasefield_output - thisSize => porosity_phasefield_sizePostResult - case default - valid = .false. - end select - if (valid) then - write(FILEUNIT,'(a)') '(porosity)'//char(9)//trim(outputName) - if (porosity_type(p) /= POROSITY_none_ID) then - do e = 1,thisNoutput(i) - write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - endif - i = hydrogenflux_typeInstance(p) ! which instance of this hydrogen flux type - valid = .true. ! assume valid - select case(hydrogenflux_type(p)) ! split per hydrogen flux type - case (HYDROGENFLUX_isoconc_ID) - outputName = HYDROGENFLUX_isoconc_label - thisNoutput => null() - thisOutput => null() - thisSize => null() - case (HYDROGENFLUX_cahnhilliard_ID) - outputName = HYDROGENFLUX_cahnhilliard_label - thisNoutput => hydrogenflux_cahnhilliard_Noutput - thisOutput => hydrogenflux_cahnhilliard_output - thisSize => hydrogenflux_cahnhilliard_sizePostResult - case default - valid = .false. - end select - if (valid) then - write(FILEUNIT,'(a)') '(hydrogenflux)'//char(9)//trim(outputName) - if (hydrogenflux_type(p) /= HYDROGENFLUX_isoconc_ID) then - do e = 1,thisNoutput(i) - write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) - enddo - endif - endif endif enddo close(FILEUNIT) @@ -383,25 +269,16 @@ subroutine homogenization_init homogenization_maxSizePostResults = 0_pInt thermal_maxSizePostResults = 0_pInt damage_maxSizePostResults = 0_pInt - vacancyflux_maxSizePostResults = 0_pInt - porosity_maxSizePostResults = 0_pInt - hydrogenflux_maxSizePostResults = 0_pInt do p = 1,size(config_homogenization) homogenization_maxSizePostResults = max(homogenization_maxSizePostResults,homogState (p)%sizePostResults) thermal_maxSizePostResults = max(thermal_maxSizePostResults, thermalState (p)%sizePostResults) damage_maxSizePostResults = max(damage_maxSizePostResults ,damageState (p)%sizePostResults) - vacancyflux_maxSizePostResults = max(vacancyflux_maxSizePostResults ,vacancyfluxState (p)%sizePostResults) - porosity_maxSizePostResults = max(porosity_maxSizePostResults ,porosityState (p)%sizePostResults) - hydrogenflux_maxSizePostResults = max(hydrogenflux_maxSizePostResults ,hydrogenfluxState(p)%sizePostResults) enddo materialpoint_sizeResults = 1 & ! grain count + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + thermal_maxSizePostResults & + damage_maxSizePostResults & - + vacancyflux_maxSizePostResults & - + porosity_maxSizePostResults & - + hydrogenflux_maxSizePostResults & + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + 1 + constitutive_plasticity_maxSizePostResults & ! constitutive size & constitutive results + constitutive_source_maxSizePostResults) @@ -460,9 +337,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) homogState, & thermalState, & damageState, & - vacancyfluxState, & - porosityState, & - hydrogenfluxState, & phase_Nsources, & mappingHomogenization, & phaseAt, phasememberAt, & @@ -569,18 +443,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) damageState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & damageState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & damageState(mappingHomogenization(2,i,e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal damage state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - vacancyfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - vacancyfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & - vacancyfluxState(mappingHomogenization(2,i,e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal vacancy transport state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - porosityState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - porosityState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & - porosityState(mappingHomogenization(2,i,e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal porosity state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & - hydrogenfluxState(mappingHomogenization(2,i,e))%State0( :,mappingHomogenization(1,i,e)) ! ...internal hydrogen transport state enddo NiterationHomog = 0_pInt @@ -654,18 +516,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) damageState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & damageState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & damageState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) ! ...internal damage state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - vacancyfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - vacancyfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & - vacancyfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal vacancy transport state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - porosityState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - porosityState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & - porosityState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal porosity state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & - hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad endif steppingNeeded @@ -729,18 +579,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) damageState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & damageState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) = & damageState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) ! ...internal damage state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - vacancyfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - vacancyfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) = & - vacancyfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal vacancy transport state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - porosityState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - porosityState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) = & - porosityState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal porosity state - forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), & - hydrogenfluxState(mappingHomogenization(2,i,e))%sizeState > 0_pInt) & - hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e)) = & - hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e))! ...internal hydrogen transport state endif endif converged @@ -846,9 +684,6 @@ subroutine materialpoint_postResults homogState, & thermalState, & damageState, & - vacancyfluxState, & - porosityState, & - hydrogenfluxState, & plasticState, & sourceState, & material_phase, & @@ -877,10 +712,7 @@ subroutine materialpoint_postResults theSize = homogState (mappingHomogenization(2,i,e))%sizePostResults & + thermalState (mappingHomogenization(2,i,e))%sizePostResults & - + damageState (mappingHomogenization(2,i,e))%sizePostResults & - + vacancyfluxState (mappingHomogenization(2,i,e))%sizePostResults & - + porosityState (mappingHomogenization(2,i,e))%sizePostResults & - + hydrogenfluxState(mappingHomogenization(2,i,e))%sizePostResults + + damageState (mappingHomogenization(2,i,e))%sizePostResults materialpoint_results(thePos+1,i,e) = real(theSize,pReal) ! tell size of homogenization results thePos = thePos + 1_pInt @@ -964,12 +796,10 @@ function homogenization_updateState(ip,el) homogenization_type, & thermal_type, & damage_type, & - vacancyflux_type, & homogenization_maxNgrains, & HOMOGENIZATION_RGC_ID, & THERMAL_adiabatic_ID, & - DAMAGE_local_ID, & - VACANCYFLUX_isochempot_ID + DAMAGE_local_ID use crystallite, only: & crystallite_P, & crystallite_dPdF, & @@ -981,8 +811,6 @@ function homogenization_updateState(ip,el) thermal_adiabatic_updateState use damage_local, only: & damage_local_updateState - use vacancyflux_isochempot, only: & - vacancyflux_isochempot_updateState implicit none integer(pInt), intent(in) :: & @@ -1023,15 +851,6 @@ function homogenization_updateState(ip,el) el) end select chosenDamage - chosenVacancyflux: select case (vacancyflux_type(mesh_element(3,el))) - case (VACANCYFLUX_isochempot_ID) chosenVacancyflux - homogenization_updateState = & - homogenization_updateState .and. & - vacancyflux_isochempot_updateState(materialpoint_subdt(ip,el), & - ip, & - el) - end select chosenVacancyflux - end function homogenization_updateState @@ -1095,15 +914,9 @@ function homogenization_postResults(ip,el) homogState, & thermalState, & damageState, & - vacancyfluxState, & - porosityState, & - hydrogenfluxState, & homogenization_type, & thermal_type, & damage_type, & - vacancyflux_type, & - porosity_type, & - hydrogenflux_type, & HOMOGENIZATION_NONE_ID, & HOMOGENIZATION_ISOSTRAIN_ID, & HOMOGENIZATION_RGC_ID, & @@ -1112,14 +925,7 @@ function homogenization_postResults(ip,el) THERMAL_conduction_ID, & DAMAGE_none_ID, & DAMAGE_local_ID, & - DAMAGE_nonlocal_ID, & - VACANCYFLUX_isoconc_ID, & - VACANCYFLUX_isochempot_ID, & - VACANCYFLUX_cahnhilliard_ID, & - POROSITY_none_ID, & - POROSITY_phasefield_ID, & - HYDROGENFLUX_isoconc_ID, & - HYDROGENFLUX_cahnhilliard_ID + DAMAGE_nonlocal_ID use homogenization_isostrain, only: & homogenization_isostrain_postResults use homogenization_RGC, only: & @@ -1132,14 +938,6 @@ function homogenization_postResults(ip,el) damage_local_postResults use damage_nonlocal, only: & damage_nonlocal_postResults - use vacancyflux_isochempot, only: & - vacancyflux_isochempot_postResults - use vacancyflux_cahnhilliard, only: & - vacancyflux_cahnhilliard_postResults - use porosity_phasefield, only: & - porosity_phasefield_postResults - use hydrogenflux_cahnhilliard, only: & - hydrogenflux_cahnhilliard_postResults implicit none integer(pInt), intent(in) :: & @@ -1147,10 +945,7 @@ function homogenization_postResults(ip,el) el !< element number real(pReal), dimension( homogState (mappingHomogenization(2,ip,el))%sizePostResults & + thermalState (mappingHomogenization(2,ip,el))%sizePostResults & - + damageState (mappingHomogenization(2,ip,el))%sizePostResults & - + vacancyfluxState (mappingHomogenization(2,ip,el))%sizePostResults & - + porosityState (mappingHomogenization(2,ip,el))%sizePostResults & - + hydrogenfluxState(mappingHomogenization(2,ip,el))%sizePostResults) :: & + + damageState (mappingHomogenization(2,ip,el))%sizePostResults) :: & homogenization_postResults integer(pInt) :: & startPos, endPos @@ -1205,39 +1000,6 @@ function homogenization_postResults(ip,el) damage_nonlocal_postResults(ip, el) end select chosenDamage - startPos = endPos + 1_pInt - endPos = endPos + vacancyfluxState(mappingHomogenization(2,ip,el))%sizePostResults - chosenVacancyflux: select case (vacancyflux_type(mesh_element(3,el))) - case (VACANCYFLUX_isoconc_ID) chosenVacancyflux - - case (VACANCYFLUX_isochempot_ID) chosenVacancyflux - homogenization_postResults(startPos:endPos) = & - vacancyflux_isochempot_postResults(ip, el) - case (VACANCYFLUX_cahnhilliard_ID) chosenVacancyflux - homogenization_postResults(startPos:endPos) = & - vacancyflux_cahnhilliard_postResults(ip, el) - end select chosenVacancyflux - - startPos = endPos + 1_pInt - endPos = endPos + porosityState(mappingHomogenization(2,ip,el))%sizePostResults - chosenPorosity: select case (porosity_type(mesh_element(3,el))) - case (POROSITY_none_ID) chosenPorosity - - case (POROSITY_phasefield_ID) chosenPorosity - homogenization_postResults(startPos:endPos) = & - porosity_phasefield_postResults(ip, el) - end select chosenPorosity - - startPos = endPos + 1_pInt - endPos = endPos + hydrogenfluxState(mappingHomogenization(2,ip,el))%sizePostResults - chosenHydrogenflux: select case (hydrogenflux_type(mesh_element(3,el))) - case (HYDROGENFLUX_isoconc_ID) chosenHydrogenflux - - case (HYDROGENFLUX_cahnhilliard_ID) chosenHydrogenflux - homogenization_postResults(startPos:endPos) = & - hydrogenflux_cahnhilliard_postResults(ip, el) - end select chosenHydrogenflux - end function homogenization_postResults end module homogenization diff --git a/src/hydrogenflux_cahnhilliard.f90 b/src/hydrogenflux_cahnhilliard.f90 deleted file mode 100644 index 3a42a49e1..000000000 --- a/src/hydrogenflux_cahnhilliard.f90 +++ /dev/null @@ -1,508 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for conservative transport of solute hydrogen -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module hydrogenflux_cahnhilliard - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - hydrogenflux_cahnhilliard_sizePostResults !< cumulative size of post results - - integer(pInt), dimension(:,:), allocatable, target, public :: & - hydrogenflux_cahnhilliard_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - hydrogenflux_cahnhilliard_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - hydrogenflux_cahnhilliard_Noutput !< number of outputs per instance of this damage - - real(pReal), parameter, private :: & - kB = 1.3806488e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: undefined_ID, & - hydrogenConc_ID - end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - hydrogenflux_cahnhilliard_outputID !< ID of each post result output - - - public :: & - hydrogenflux_cahnhilliard_init, & - hydrogenflux_cahnhilliard_getMobility33, & - hydrogenflux_cahnhilliard_getDiffusion33, & - hydrogenflux_cahnhilliard_getFormationEnergy, & - hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent, & - hydrogenflux_cahnhilliard_getChemPotAndItsTangent, & - hydrogenflux_cahnhilliard_putHydrogenConcAndItsRate, & - hydrogenflux_cahnhilliard_postResults - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine hydrogenflux_cahnhilliard_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - hydrogenflux_type, & - hydrogenflux_typeInstance, & - homogenization_Noutput, & - HYDROGENFLUX_cahnhilliard_label, & - HYDROGENFLUX_cahnhilliard_ID, & - material_homog, & - mappingHomogenization, & - hydrogenfluxState, & - hydrogenfluxMapping, & - hydrogenConc, & - hydrogenConcRate, & - hydrogenflux_initialCh - use config, only: & - material_partHomogenization, & - material_partPhase - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o - integer(pInt) :: sizeState - integer(pInt) :: NofMyHomog - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- hydrogenflux_'//HYDROGENFLUX_cahnhilliard_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(hydrogenflux_type == HYDROGENFLUX_cahnhilliard_ID),pInt) - if (maxNinstance == 0_pInt) return - - allocate(hydrogenflux_cahnhilliard_sizePostResults(maxNinstance), source=0_pInt) - allocate(hydrogenflux_cahnhilliard_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) - allocate(hydrogenflux_cahnhilliard_output (maxval(homogenization_Noutput),maxNinstance)) - hydrogenflux_cahnhilliard_output = '' - allocate(hydrogenflux_cahnhilliard_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) - allocate(hydrogenflux_cahnhilliard_Noutput (maxNinstance), source=0_pInt) - - rewind(fileUnit) - section = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to - line = IO_read(fileUnit) - enddo - - parsingHomog: do while (trim(line) /= IO_EOF) ! read through sections of homog part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next homog section - section = section + 1_pInt ! advance homog section counter - cycle ! skip to next line - endif - - if (section > 0_pInt ) then; if (hydrogenflux_type(section) == HYDROGENFLUX_cahnhilliard_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = hydrogenflux_typeInstance(section) ! which instance of my hydrogenflux is present homog - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('hydrogenconc') - hydrogenflux_cahnhilliard_Noutput(instance) = hydrogenflux_cahnhilliard_Noutput(instance) + 1_pInt - hydrogenflux_cahnhilliard_outputID(hydrogenflux_cahnhilliard_Noutput(instance),instance) = hydrogenConc_ID - hydrogenflux_cahnhilliard_output(hydrogenflux_cahnhilliard_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - end select - - end select - endif; endif - enddo parsingHomog - - rewind(fileUnit) - section = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - initializeInstances: do section = 1_pInt, size(hydrogenflux_type) - if (hydrogenflux_type(section) == HYDROGENFLUX_cahnhilliard_ID) then - NofMyHomog=count(material_homog==section) - instance = hydrogenflux_typeInstance(section) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,hydrogenflux_cahnhilliard_Noutput(instance) - select case(hydrogenflux_cahnhilliard_outputID(o,instance)) - case(hydrogenConc_ID) - mySize = 1_pInt - end select - - if (mySize > 0_pInt) then ! any meaningful output found - hydrogenflux_cahnhilliard_sizePostResult(o,instance) = mySize - hydrogenflux_cahnhilliard_sizePostResults(instance) = hydrogenflux_cahnhilliard_sizePostResults(instance) + mySize - endif - enddo outputsLoop - -! allocate state arrays - sizeState = 0_pInt - hydrogenfluxState(section)%sizeState = sizeState - hydrogenfluxState(section)%sizePostResults = hydrogenflux_cahnhilliard_sizePostResults(instance) - allocate(hydrogenfluxState(section)%state0 (sizeState,NofMyHomog)) - allocate(hydrogenfluxState(section)%subState0(sizeState,NofMyHomog)) - allocate(hydrogenfluxState(section)%state (sizeState,NofMyHomog)) - - nullify(hydrogenfluxMapping(section)%p) - hydrogenfluxMapping(section)%p => mappingHomogenization(1,:,:) - deallocate(hydrogenConc (section)%p) - deallocate(hydrogenConcRate(section)%p) - allocate (hydrogenConc (section)%p(NofMyHomog), source=hydrogenflux_initialCh(section)) - allocate (hydrogenConcRate(section)%p(NofMyHomog), source=0.0_pReal) - - endif - - enddo initializeInstances - -end subroutine hydrogenflux_cahnhilliard_init - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized solute mobility tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function hydrogenflux_cahnhilliard_getMobility33(ip,el) - use lattice, only: & - lattice_hydrogenfluxMobility33 - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_push33ToRef - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - hydrogenflux_cahnhilliard_getMobility33 - integer(pInt) :: & - grain - - hydrogenflux_cahnhilliard_getMobility33 = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - hydrogenflux_cahnhilliard_getMobility33 = hydrogenflux_cahnhilliard_getMobility33 + & - crystallite_push33ToRef(grain,ip,el,lattice_hydrogenfluxMobility33(:,:,material_phase(grain,ip,el))) - enddo - - hydrogenflux_cahnhilliard_getMobility33 = & - hydrogenflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function hydrogenflux_cahnhilliard_getMobility33 - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized solute nonlocal diffusion tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function hydrogenflux_cahnhilliard_getDiffusion33(ip,el) - use lattice, only: & - lattice_hydrogenfluxDiffusion33 - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_push33ToRef - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - hydrogenflux_cahnhilliard_getDiffusion33 - integer(pInt) :: & - grain - - hydrogenflux_cahnhilliard_getDiffusion33 = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - hydrogenflux_cahnhilliard_getDiffusion33 = hydrogenflux_cahnhilliard_getDiffusion33 + & - crystallite_push33ToRef(grain,ip,el,lattice_hydrogenfluxDiffusion33(:,:,material_phase(grain,ip,el))) - enddo - - hydrogenflux_cahnhilliard_getDiffusion33 = & - hydrogenflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function hydrogenflux_cahnhilliard_getDiffusion33 - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized solution energy -!-------------------------------------------------------------------------------------------------- -function hydrogenflux_cahnhilliard_getFormationEnergy(ip,el) - use lattice, only: & - lattice_hydrogenFormationEnergy, & - lattice_hydrogenVol, & - lattice_hydrogenSurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - hydrogenflux_cahnhilliard_getFormationEnergy - integer(pInt) :: & - grain - - hydrogenflux_cahnhilliard_getFormationEnergy = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - hydrogenflux_cahnhilliard_getFormationEnergy = hydrogenflux_cahnhilliard_getFormationEnergy + & - lattice_hydrogenFormationEnergy(material_phase(grain,ip,el))/ & - lattice_hydrogenVol(material_phase(grain,ip,el))/ & - lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el)) - enddo - - hydrogenflux_cahnhilliard_getFormationEnergy = & - hydrogenflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function hydrogenflux_cahnhilliard_getFormationEnergy - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized hydrogen entropy coefficient -!-------------------------------------------------------------------------------------------------- -function hydrogenflux_cahnhilliard_getEntropicCoeff(ip,el) - use lattice, only: & - lattice_hydrogenVol, & - lattice_hydrogenSurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_homog, & - material_phase, & - temperature, & - thermalMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - hydrogenflux_cahnhilliard_getEntropicCoeff - integer(pInt) :: & - grain - - hydrogenflux_cahnhilliard_getEntropicCoeff = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homog(ip,el)) - hydrogenflux_cahnhilliard_getEntropicCoeff = hydrogenflux_cahnhilliard_getEntropicCoeff + & - kB/ & - lattice_hydrogenVol(material_phase(grain,ip,el))/ & - lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el)) - enddo - - hydrogenflux_cahnhilliard_getEntropicCoeff = hydrogenflux_cahnhilliard_getEntropicCoeff* & - temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ & - real(homogenization_Ngrains(material_homog(ip,el)),pReal) - -end function hydrogenflux_cahnhilliard_getEntropicCoeff - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized kinematic contribution to chemical potential -!-------------------------------------------------------------------------------------------------- -subroutine hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCh, Ch, ip, el) - use lattice, only: & - lattice_hydrogenSurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_homog, & - phase_kinematics, & - phase_Nkinematics, & - material_phase, & - KINEMATICS_hydrogen_strain_ID - use crystallite, only: & - crystallite_Tstar_v, & - crystallite_Fi0, & - crystallite_Fi - use kinematics_hydrogen_strain, only: & - kinematics_hydrogen_strain_ChemPotAndItsTangent - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Ch - real(pReal), intent(out) :: & - KPot, dKPot_dCh - real(pReal) :: & - my_KPot, my_dKPot_dCh - integer(pInt) :: & - grain, kinematics - - KPot = 0.0_pReal - dKPot_dCh = 0.0_pReal - do grain = 1_pInt,homogenization_Ngrains(material_homog(ip,el)) - do kinematics = 1_pInt, phase_Nkinematics(material_phase(grain,ip,el)) - select case (phase_kinematics(kinematics,material_phase(grain,ip,el))) - case (KINEMATICS_hydrogen_strain_ID) - call kinematics_hydrogen_strain_ChemPotAndItsTangent(my_KPot, my_dKPot_dCh, & - crystallite_Tstar_v(1:6,grain,ip,el), & - crystallite_Fi0(1:3,1:3,grain,ip,el), & - crystallite_Fi (1:3,1:3,grain,ip,el), & - grain,ip, el) - - case default - my_KPot = 0.0_pReal - my_dKPot_dCh = 0.0_pReal - - end select - KPot = KPot + my_KPot/lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el)) - dKPot_dCh = dKPot_dCh + my_dKPot_dCh/lattice_hydrogenSurfaceEnergy(material_phase(grain,ip,el)) - enddo - enddo - - KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal) - dKPot_dCh = dKPot_dCh/real(homogenization_Ngrains(material_homog(ip,el)),pReal) - -end subroutine hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized chemical potential -!-------------------------------------------------------------------------------------------------- -subroutine hydrogenflux_cahnhilliard_getChemPotAndItsTangent(ChemPot,dChemPot_dCh,Ch,ip,el) - use numerics, only: & - hydrogenBoundPenalty, & - hydrogenPolyOrder - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Ch - real(pReal), intent(out) :: & - ChemPot, & - dChemPot_dCh - real(pReal) :: & - kBT, KPot, dKPot_dCh - integer(pInt) :: & - o - - ChemPot = hydrogenflux_cahnhilliard_getFormationEnergy(ip,el) - dChemPot_dCh = 0.0_pReal - kBT = hydrogenflux_cahnhilliard_getEntropicCoeff(ip,el) - do o = 1_pInt, hydrogenPolyOrder - ChemPot = ChemPot + kBT*((2.0_pReal*Ch - 1.0_pReal)**real(2_pInt*o-1_pInt,pReal))/ & - real(2_pInt*o-1_pInt,pReal) - dChemPot_dCh = dChemPot_dCh + 2.0_pReal*kBT*(2.0_pReal*Ch - 1.0_pReal)**real(2_pInt*o-2_pInt,pReal) - enddo - - call hydrogenflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCh, Ch, ip, el) - ChemPot = ChemPot + KPot - dChemPot_dCh = dChemPot_dCh + dKPot_dCh - - if (Ch < 0.0_pReal) then - ChemPot = ChemPot - 3.0_pReal*hydrogenBoundPenalty*Ch*Ch - dChemPot_dCh = dChemPot_dCh - 6.0_pReal*hydrogenBoundPenalty*Ch - elseif (Ch > 1.0_pReal) then - ChemPot = ChemPot + 3.0_pReal*hydrogenBoundPenalty*(1.0_pReal - Ch)*(1.0_pReal - Ch) - dChemPot_dCh = dChemPot_dCh - 6.0_pReal*hydrogenBoundPenalty*(1.0_pReal - Ch) - endif - -end subroutine hydrogenflux_cahnhilliard_getChemPotAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief updates hydrogen concentration with solution from Cahn-Hilliard PDE for solute transport -!-------------------------------------------------------------------------------------------------- -subroutine hydrogenflux_cahnhilliard_putHydrogenConcAndItsRate(Ch,Chdot,ip,el) - use material, only: & - mappingHomogenization, & - hydrogenConc, & - hydrogenConcRate, & - hydrogenfluxMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Ch, & - Chdot - integer(pInt) :: & - homog, & - offset - - homog = mappingHomogenization(2,ip,el) - offset = hydrogenfluxMapping(homog)%p(ip,el) - hydrogenConc (homog)%p(offset) = Ch - hydrogenConcRate(homog)%p(offset) = Chdot - -end subroutine hydrogenflux_cahnhilliard_putHydrogenConcAndItsRate - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of hydrogen transport results -!-------------------------------------------------------------------------------------------------- -function hydrogenflux_cahnhilliard_postResults(ip,el) - use material, only: & - mappingHomogenization, & - hydrogenflux_typeInstance, & - hydrogenConc, & - hydrogenfluxMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(hydrogenflux_cahnhilliard_sizePostResults(hydrogenflux_typeInstance(mappingHomogenization(2,ip,el)))) :: & - hydrogenflux_cahnhilliard_postResults - - integer(pInt) :: & - instance, homog, offset, o, c - - homog = mappingHomogenization(2,ip,el) - offset = hydrogenfluxMapping(homog)%p(ip,el) - instance = hydrogenflux_typeInstance(homog) - - c = 0_pInt - hydrogenflux_cahnhilliard_postResults = 0.0_pReal - - do o = 1_pInt,hydrogenflux_cahnhilliard_Noutput(instance) - select case(hydrogenflux_cahnhilliard_outputID(o,instance)) - - case (hydrogenConc_ID) - hydrogenflux_cahnhilliard_postResults(c+1_pInt) = hydrogenConc(homog)%p(offset) - c = c + 1 - end select - enddo -end function hydrogenflux_cahnhilliard_postResults - -end module hydrogenflux_cahnhilliard diff --git a/src/hydrogenflux_isoconc.f90 b/src/hydrogenflux_isoconc.f90 deleted file mode 100644 index 836d29198..000000000 --- a/src/hydrogenflux_isoconc.f90 +++ /dev/null @@ -1,62 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for constant hydrogen concentration -!-------------------------------------------------------------------------------------------------- -module hydrogenflux_isoconc - - implicit none - private - - public :: & - hydrogenflux_isoconc_init - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields, reads information from material configuration file -!-------------------------------------------------------------------------------------------------- -subroutine hydrogenflux_isoconc_init() -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use prec, only: & - pReal, & - pInt - use IO, only: & - IO_timeStamp - use material - use config - - implicit none - integer(pInt) :: & - homog, & - NofMyHomog - - write(6,'(/,a)') ' <<<+- hydrogenflux_'//HYDROGENFLUX_isoconc_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - initializeInstances: do homog = 1_pInt, material_Nhomogenization - - myhomog: if (hydrogenflux_type(homog) == HYDROGENFLUX_isoconc_ID) then - NofMyHomog = count(material_homog == homog) - hydrogenfluxState(homog)%sizeState = 0_pInt - hydrogenfluxState(homog)%sizePostResults = 0_pInt - allocate(hydrogenfluxState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal) - allocate(hydrogenfluxState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal) - allocate(hydrogenfluxState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal) - - deallocate(hydrogenConc (homog)%p) - deallocate(hydrogenConcRate(homog)%p) - allocate (hydrogenConc (homog)%p(1), source=hydrogenflux_initialCh(homog)) - allocate (hydrogenConcRate(homog)%p(1), source=0.0_pReal) - - endif myhomog - enddo initializeInstances - - -end subroutine hydrogenflux_isoconc_init - -end module hydrogenflux_isoconc diff --git a/src/kinematics_hydrogen_strain.f90 b/src/kinematics_hydrogen_strain.f90 deleted file mode 100644 index 516ca286f..000000000 --- a/src/kinematics_hydrogen_strain.f90 +++ /dev/null @@ -1,263 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incorporating kinematics resulting from interstitial hydrogen -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module kinematics_hydrogen_strain - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - kinematics_hydrogen_strain_sizePostResults, & !< cumulative size of post results - kinematics_hydrogen_strain_offset, & !< which kinematics is my current damage mechanism? - kinematics_hydrogen_strain_instance !< instance of damage kinematics mechanism - - integer(pInt), dimension(:,:), allocatable, target, public :: & - kinematics_hydrogen_strain_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - kinematics_hydrogen_strain_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - kinematics_hydrogen_strain_Noutput !< number of outputs per instance of this damage - - real(pReal), dimension(:), allocatable, private :: & - kinematics_hydrogen_strain_coeff - - public :: & - kinematics_hydrogen_strain_init, & - kinematics_hydrogen_strain_initialStrain, & - kinematics_hydrogen_strain_LiAndItsTangent, & - kinematics_hydrogen_strain_ChemPotAndItsTangent - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine kinematics_hydrogen_strain_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - phase_kinematics, & - phase_Nkinematics, & - phase_Noutput, & - KINEMATICS_hydrogen_strain_label, & - KINEMATICS_hydrogen_strain_ID - use config, only: & - material_Nphase, & - MATERIAL_partPhase - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt) - if (maxNinstance == 0_pInt) return - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(kinematics_hydrogen_strain_offset(material_Nphase), source=0_pInt) - allocate(kinematics_hydrogen_strain_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - kinematics_hydrogen_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_hydrogen_strain_ID) - do kinematics = 1, phase_Nkinematics(phase) - if (phase_kinematics(kinematics,phase) == kinematics_hydrogen_strain_ID) & - kinematics_hydrogen_strain_offset(phase) = kinematics - enddo - enddo - - allocate(kinematics_hydrogen_strain_sizePostResults(maxNinstance), source=0_pInt) - allocate(kinematics_hydrogen_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(kinematics_hydrogen_strain_output(maxval(phase_Noutput),maxNinstance)) - kinematics_hydrogen_strain_output = '' - allocate(kinematics_hydrogen_strain_Noutput(maxNinstance), source=0_pInt) - allocate(kinematics_hydrogen_strain_coeff(maxNinstance), source=0.0_pReal) - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_hydrogen_strain_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = kinematics_hydrogen_strain_instance(phase) ! which instance of my damage is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('hydrogen_strain_coeff') - kinematics_hydrogen_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt) - - end select - endif; endif - enddo parsingFile - -end subroutine kinematics_hydrogen_strain_init - -!-------------------------------------------------------------------------------------------------- -!> @brief report initial hydrogen strain based on current hydrogen conc deviation from -!> equillibrium (0) -!-------------------------------------------------------------------------------------------------- -pure function kinematics_hydrogen_strain_initialStrain(ipc, ip, el) - use math, only: & - math_I3 - use material, only: & - material_phase, & - material_homog, & - hydrogenConc, & - hydrogenfluxMapping - use lattice, only: & - lattice_equilibriumHydrogenConcentration - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - kinematics_hydrogen_strain_initialStrain !< initial thermal strain (should be small strain, though) - integer(pInt) :: & - phase, & - homog, offset, instance - - phase = material_phase(ipc,ip,el) - instance = kinematics_hydrogen_strain_instance(phase) - homog = material_homog(ip,el) - offset = hydrogenfluxMapping(homog)%p(ip,el) - - kinematics_hydrogen_strain_initialStrain = & - (hydrogenConc(homog)%p(offset) - lattice_equilibriumHydrogenConcentration(phase)) * & - kinematics_hydrogen_strain_coeff(instance)* math_I3 - -end function kinematics_hydrogen_strain_initialStrain - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -!-------------------------------------------------------------------------------------------------- -subroutine kinematics_hydrogen_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el) - use material, only: & - material_phase, & - material_homog, & - hydrogenConc, & - hydrogenConcRate, & - hydrogenfluxMapping - use math, only: & - math_I3 - use lattice, only: & - lattice_equilibriumHydrogenConcentration - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(out), dimension(3,3) :: & - Li !< thermal velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor) - integer(pInt) :: & - phase, & - instance, & - homog, offset - real(pReal) :: & - Ch, ChEq, ChDot - - phase = material_phase(ipc,ip,el) - instance = kinematics_hydrogen_strain_instance(phase) - homog = material_homog(ip,el) - offset = hydrogenfluxMapping(homog)%p(ip,el) - Ch = hydrogenConc(homog)%p(offset) - ChDot = hydrogenConcRate(homog)%p(offset) - ChEq = lattice_equilibriumHydrogenConcentration(phase) - - Li = ChDot*math_I3* & - kinematics_hydrogen_strain_coeff(instance)/ & - (1.0_pReal + kinematics_hydrogen_strain_coeff(instance)*(Ch - ChEq)) - dLi_dTstar3333 = 0.0_pReal - -end subroutine kinematics_hydrogen_strain_LiAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the kinematic contribution to hydrogen chemical potential -!-------------------------------------------------------------------------------------------------- -subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCh, Tstar_v, Fi0, Fi, ipc, ip, el) - use material, only: & - material_phase - use math, only: & - math_inv33, & - math_mul33x33, & - math_Mandel6to33, & - math_transpose33 - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(6) :: & - Tstar_v - real(pReal), intent(in), dimension(3,3) :: & - Fi0, Fi - real(pReal), intent(out) :: & - ChemPot, dChemPot_dCh - integer(pInt) :: & - phase, & - instance - - phase = material_phase(ipc,ip,el) - instance = kinematics_hydrogen_strain_instance(phase) - - ChemPot = -kinematics_hydrogen_strain_coeff(instance)* & - sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* & - math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi)) - dChemPot_dCh = 0.0_pReal - -end subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent - -end module kinematics_hydrogen_strain diff --git a/src/kinematics_vacancy_strain.f90 b/src/kinematics_vacancy_strain.f90 deleted file mode 100644 index 7ecc7fe6e..000000000 --- a/src/kinematics_vacancy_strain.f90 +++ /dev/null @@ -1,264 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incorporating kinematics resulting from vacancy point defects -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module kinematics_vacancy_strain - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - kinematics_vacancy_strain_sizePostResults, & !< cumulative size of post results - kinematics_vacancy_strain_offset, & !< which kinematics is my current damage mechanism? - kinematics_vacancy_strain_instance !< instance of damage kinematics mechanism - - integer(pInt), dimension(:,:), allocatable, target, public :: & - kinematics_vacancy_strain_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - kinematics_vacancy_strain_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - kinematics_vacancy_strain_Noutput !< number of outputs per instance of this damage - - real(pReal), dimension(:), allocatable, private :: & - kinematics_vacancy_strain_coeff - - public :: & - kinematics_vacancy_strain_init, & - kinematics_vacancy_strain_initialStrain, & - kinematics_vacancy_strain_LiAndItsTangent, & - kinematics_vacancy_strain_ChemPotAndItsTangent - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine kinematics_vacancy_strain_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - phase_kinematics, & - phase_Nkinematics, & - phase_Noutput, & - KINEMATICS_vacancy_strain_label, & - KINEMATICS_vacancy_strain_ID - use config, only: & - material_Nphase, & - MATERIAL_partPhase - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,kinematics - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_vacancy_strain_LABEL//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_kinematics == KINEMATICS_vacancy_strain_ID),pInt) - if (maxNinstance == 0_pInt) return - - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(kinematics_vacancy_strain_offset(material_Nphase), source=0_pInt) - allocate(kinematics_vacancy_strain_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - kinematics_vacancy_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_vacancy_strain_ID) - do kinematics = 1, phase_Nkinematics(phase) - if (phase_kinematics(kinematics,phase) == kinematics_vacancy_strain_ID) & - kinematics_vacancy_strain_offset(phase) = kinematics - enddo - enddo - - allocate(kinematics_vacancy_strain_sizePostResults(maxNinstance), source=0_pInt) - allocate(kinematics_vacancy_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(kinematics_vacancy_strain_output(maxval(phase_Noutput),maxNinstance)) - kinematics_vacancy_strain_output = '' - allocate(kinematics_vacancy_strain_Noutput(maxNinstance), source=0_pInt) - allocate(kinematics_vacancy_strain_coeff(maxNinstance), source=0.0_pReal) - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_vacancy_strain_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - instance = kinematics_vacancy_strain_instance(phase) ! which instance of my damage is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('vacancy_strain_coeff') - kinematics_vacancy_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt) - - end select - endif; endif - enddo parsingFile - -end subroutine kinematics_vacancy_strain_init - -!-------------------------------------------------------------------------------------------------- -!> @brief report initial vacancy strain based on current vacancy conc deviation from equillibrium -!-------------------------------------------------------------------------------------------------- -pure function kinematics_vacancy_strain_initialStrain(ipc, ip, el) - use math, only: & - math_I3 - use material, only: & - material_phase, & - material_homog, & - vacancyConc, & - vacancyfluxMapping - use lattice, only: & - lattice_equilibriumVacancyConcentration - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - kinematics_vacancy_strain_initialStrain !< initial thermal strain (should be small strain, though) - integer(pInt) :: & - phase, & - homog, offset, instance - - phase = material_phase(ipc,ip,el) - instance = kinematics_vacancy_strain_instance(phase) - homog = material_homog(ip,el) - offset = vacancyfluxMapping(homog)%p(ip,el) - - kinematics_vacancy_strain_initialStrain = & - (vacancyConc(homog)%p(offset) - lattice_equilibriumVacancyConcentration(phase)) * & - kinematics_vacancy_strain_coeff(instance)* math_I3 - -end function kinematics_vacancy_strain_initialStrain - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -!-------------------------------------------------------------------------------------------------- -subroutine kinematics_vacancy_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el) - use material, only: & - material_phase, & - material_homog, & - vacancyConc, & - vacancyConcRate, & - vacancyfluxMapping - use math, only: & - math_I3 - use lattice, only: & - lattice_equilibriumVacancyConcentration - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(out), dimension(3,3) :: & - Li !< thermal velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dTstar3333 !< derivative of Li with respect to Tstar (4th-order tensor) - integer(pInt) :: & - phase, & - instance, & - homog, offset - real(pReal) :: & - Cv, CvEq, CvDot - - phase = material_phase(ipc,ip,el) - instance = kinematics_vacancy_strain_instance(phase) - homog = material_homog(ip,el) - offset = vacancyfluxMapping(homog)%p(ip,el) - - Cv = vacancyConc(homog)%p(offset) - CvDot = vacancyConcRate(homog)%p(offset) - CvEq = lattice_equilibriumvacancyConcentration(phase) - - Li = CvDot*math_I3* & - kinematics_vacancy_strain_coeff(instance)/ & - (1.0_pReal + kinematics_vacancy_strain_coeff(instance)*(Cv - CvEq)) - - dLi_dTstar3333 = 0.0_pReal - -end subroutine kinematics_vacancy_strain_LiAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the kinematic contribution to vacancy chemical potential -!-------------------------------------------------------------------------------------------------- -subroutine kinematics_vacancy_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCv, Tstar_v, Fi0, Fi, ipc, ip, el) - use material, only: & - material_phase - use math, only: & - math_inv33, & - math_mul33x33, & - math_Mandel6to33, & - math_transpose33 - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(in), dimension(6) :: & - Tstar_v - real(pReal), intent(in), dimension(3,3) :: & - Fi0, Fi - real(pReal), intent(out) :: & - ChemPot, dChemPot_dCv - integer(pInt) :: & - phase, & - instance - - phase = material_phase(ipc,ip,el) - instance = kinematics_vacancy_strain_instance(phase) - - ChemPot = -kinematics_vacancy_strain_coeff(instance)* & - sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* & - math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi)) - dChemPot_dCv = 0.0_pReal - -end subroutine kinematics_vacancy_strain_ChemPotAndItsTangent - -end module kinematics_vacancy_strain diff --git a/src/material.f90 b/src/material.f90 index c9dd28079..8356f43c7 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -36,29 +36,16 @@ module material SOURCE_damage_isoDuctile_label = 'damage_isoductile', & SOURCE_damage_anisoBrittle_label = 'damage_anisobrittle', & SOURCE_damage_anisoDuctile_label = 'damage_anisoductile', & - SOURCE_vacancy_phenoplasticity_label = 'vacancy_phenoplasticity', & - SOURCE_vacancy_irradiation_label = 'vacancy_irradiation', & - SOURCE_vacancy_thermalfluc_label = 'vacancy_thermalfluctuation', & KINEMATICS_thermal_expansion_label = 'thermal_expansion', & KINEMATICS_cleavage_opening_label = 'cleavage_opening', & KINEMATICS_slipplane_opening_label = 'slipplane_opening', & - KINEMATICS_vacancy_strain_label = 'vacancy_strain', & - KINEMATICS_hydrogen_strain_label = 'hydrogen_strain', & STIFFNESS_DEGRADATION_damage_label = 'damage', & - STIFFNESS_DEGRADATION_porosity_label = 'porosity', & THERMAL_isothermal_label = 'isothermal', & THERMAL_adiabatic_label = 'adiabatic', & THERMAL_conduction_label = 'conduction', & DAMAGE_none_label = 'none', & DAMAGE_local_label = 'local', & DAMAGE_nonlocal_label = 'nonlocal', & - VACANCYFLUX_isoconc_label = 'isoconcentration', & - VACANCYFLUX_isochempot_label = 'isochemicalpotential', & - VACANCYFLUX_cahnhilliard_label = 'cahnhilliard', & - POROSITY_none_label = 'none', & - POROSITY_phasefield_label = 'phasefield', & - HYDROGENFLUX_isoconc_label = 'isoconcentration', & - HYDROGENFLUX_cahnhilliard_label = 'cahnhilliard', & HOMOGENIZATION_none_label = 'none', & HOMOGENIZATION_isostrain_label = 'isostrain', & HOMOGENIZATION_rgc_label = 'rgc' @@ -87,25 +74,19 @@ module material SOURCE_damage_isoBrittle_ID, & SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & - SOURCE_damage_anisoDuctile_ID, & - SOURCE_vacancy_phenoplasticity_ID, & - SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID + SOURCE_damage_anisoDuctile_ID end enum enum, bind(c) enumerator :: KINEMATICS_undefined_ID, & KINEMATICS_cleavage_opening_ID, & KINEMATICS_slipplane_opening_ID, & - KINEMATICS_thermal_expansion_ID, & - KINEMATICS_vacancy_strain_ID, & - KINEMATICS_hydrogen_strain_ID + KINEMATICS_thermal_expansion_ID end enum enum, bind(c) enumerator :: STIFFNESS_DEGRADATION_undefined_ID, & - STIFFNESS_DEGRADATION_damage_ID, & - STIFFNESS_DEGRADATION_porosity_ID + STIFFNESS_DEGRADATION_damage_ID end enum enum, bind(c) @@ -120,21 +101,6 @@ module material DAMAGE_nonlocal_ID end enum - enum, bind(c) - enumerator :: VACANCYFLUX_isoconc_ID, & - VACANCYFLUX_isochempot_ID, & - VACANCYFLUX_cahnhilliard_ID - end enum - - enum, bind(c) - enumerator :: POROSITY_none_ID, & - POROSITY_phasefield_ID - end enum - enum, bind(c) - enumerator :: HYDROGENFLUX_isoconc_ID, & - HYDROGENFLUX_cahnhilliard_ID - end enum - enum, bind(c) enumerator :: HOMOGENIZATION_undefined_ID, & HOMOGENIZATION_none_ID, & @@ -150,12 +116,6 @@ module material thermal_type !< thermal transport model integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & damage_type !< nonlocal damage model - integer(kind(VACANCYFLUX_isoconc_ID)), dimension(:), allocatable, public, protected :: & - vacancyflux_type !< vacancy transport model - integer(kind(POROSITY_none_ID)), dimension(:), allocatable, public, protected :: & - porosity_type !< porosity evolution model - integer(kind(HYDROGENFLUX_isoconc_ID)), dimension(:), allocatable, public, protected :: & - hydrogenflux_type !< hydrogen transport model integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable, public, protected :: & phase_source, & !< active sources mechanisms of each phase @@ -181,17 +141,11 @@ module material homogenization_typeInstance, & !< instance of particular type of each homogenization thermal_typeInstance, & !< instance of particular type of each thermal transport damage_typeInstance, & !< instance of particular type of each nonlocal damage - vacancyflux_typeInstance, & !< instance of particular type of each vacancy flux - porosity_typeInstance, & !< instance of particular type of each porosity model - hydrogenflux_typeInstance, & !< instance of particular type of each hydrogen flux microstructure_crystallite !< crystallite setting ID of each microstructure ! DEPRECATED !!!! real(pReal), dimension(:), allocatable, public, protected :: & thermal_initialT, & !< initial temperature per each homogenization - damage_initialPhi, & !< initial damage per each homogenization - vacancyflux_initialCv, & !< initial vacancy concentration per each homogenization - porosity_initialPhi, & !< initial posority per each homogenization - hydrogenflux_initialCh !< initial hydrogen concentration per each homogenization + damage_initialPhi !< initial damage per each homogenization ! NEW MAPPINGS integer(pInt), dimension(:), allocatable, public, protected :: & @@ -221,10 +175,7 @@ module material type(tState), allocatable, dimension(:), public :: & homogState, & thermalState, & - damageState, & - vacancyfluxState, & - porosityState, & - hydrogenfluxState + damageState integer(pInt), dimension(:,:,:), allocatable, public, protected :: & material_texture !< texture (index) of each grain,IP,element @@ -274,20 +225,12 @@ module material type(tHomogMapping), allocatable, dimension(:), public :: & thermalMapping, & !< mapping for thermal state/fields - damageMapping, & !< mapping for damage state/fields - vacancyfluxMapping, & !< mapping for vacancy conc state/fields - porosityMapping, & !< mapping for porosity state/fields - hydrogenfluxMapping !< mapping for hydrogen conc state/fields + damageMapping !< mapping for damage state/fields type(group_float), allocatable, dimension(:), public :: & temperature, & !< temperature field damage, & !< damage field - vacancyConc, & !< vacancy conc field - porosity, & !< porosity field - hydrogenConc, & !< hydrogen conc field - temperatureRate, & !< temperature change rate field - vacancyConcRate, & !< vacancy conc change field - hydrogenConcRate !< hydrogen conc change field + temperatureRate !< temperature change rate field public :: & material_init, & @@ -306,29 +249,16 @@ module material SOURCE_damage_isoDuctile_ID, & SOURCE_damage_anisoBrittle_ID, & SOURCE_damage_anisoDuctile_ID, & - SOURCE_vacancy_phenoplasticity_ID, & - SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID, & KINEMATICS_cleavage_opening_ID, & KINEMATICS_slipplane_opening_ID, & KINEMATICS_thermal_expansion_ID, & - KINEMATICS_vacancy_strain_ID, & - KINEMATICS_hydrogen_strain_ID, & STIFFNESS_DEGRADATION_damage_ID, & - STIFFNESS_DEGRADATION_porosity_ID, & THERMAL_isothermal_ID, & THERMAL_adiabatic_ID, & THERMAL_conduction_ID, & DAMAGE_none_ID, & DAMAGE_local_ID, & DAMAGE_nonlocal_ID, & - VACANCYFLUX_isoconc_ID, & - VACANCYFLUX_isochempot_ID, & - VACANCYFLUX_cahnhilliard_ID, & - POROSITY_none_ID, & - POROSITY_phasefield_ID, & - HYDROGENFLUX_isoconc_ID, & - HYDROGENFLUX_cahnhilliard_ID, & HOMOGENIZATION_none_ID, & HOMOGENIZATION_isostrain_ID, & HOMOGENIZATION_RGC_ID @@ -420,25 +350,14 @@ subroutine material_init() allocate(homogState (size(config_homogenization))) allocate(thermalState (size(config_homogenization))) allocate(damageState (size(config_homogenization))) - allocate(vacancyfluxState (size(config_homogenization))) - allocate(porosityState (size(config_homogenization))) - allocate(hydrogenfluxState (size(config_homogenization))) allocate(thermalMapping (size(config_homogenization))) allocate(damageMapping (size(config_homogenization))) - allocate(vacancyfluxMapping (size(config_homogenization))) - allocate(porosityMapping (size(config_homogenization))) - allocate(hydrogenfluxMapping(size(config_homogenization))) allocate(temperature (size(config_homogenization))) allocate(damage (size(config_homogenization))) - allocate(vacancyConc (size(config_homogenization))) - allocate(porosity (size(config_homogenization))) - allocate(hydrogenConc (size(config_homogenization))) allocate(temperatureRate (size(config_homogenization))) - allocate(vacancyConcRate (size(config_homogenization))) - allocate(hydrogenConcRate (size(config_homogenization))) do m = 1_pInt,size(config_microstructure) if(microstructure_crystallite(m) < 1_pInt .or. & @@ -511,17 +430,9 @@ subroutine material_init() do myHomog = 1,size(config_homogenization) thermalMapping (myHomog)%p => mappingHomogenizationConst damageMapping (myHomog)%p => mappingHomogenizationConst - vacancyfluxMapping (myHomog)%p => mappingHomogenizationConst - porosityMapping (myHomog)%p => mappingHomogenizationConst - hydrogenfluxMapping(myHomog)%p => mappingHomogenizationConst allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog)) allocate(damage (myHomog)%p(1), source=damage_initialPhi(myHomog)) - allocate(vacancyConc (myHomog)%p(1), source=vacancyflux_initialCv(myHomog)) - allocate(porosity (myHomog)%p(1), source=porosity_initialPhi(myHomog)) - allocate(hydrogenConc (myHomog)%p(1), source=hydrogenflux_initialCh(myHomog)) allocate(temperatureRate (myHomog)%p(1), source=0.0_pReal) - allocate(vacancyConcRate (myHomog)%p(1), source=0.0_pReal) - allocate(hydrogenConcRate(myHomog)%p(1), source=0.0_pReal) enddo end subroutine material_init @@ -545,23 +456,14 @@ subroutine material_parseHomogenization allocate(homogenization_type(size(config_homogenization)), source=HOMOGENIZATION_undefined_ID) allocate(thermal_type(size(config_homogenization)), source=THERMAL_isothermal_ID) allocate(damage_type (size(config_homogenization)), source=DAMAGE_none_ID) - allocate(vacancyflux_type(size(config_homogenization)), source=VACANCYFLUX_isoconc_ID) - allocate(porosity_type (size(config_homogenization)), source=POROSITY_none_ID) - allocate(hydrogenflux_type(size(config_homogenization)), source=HYDROGENFLUX_isoconc_ID) allocate(homogenization_typeInstance(size(config_homogenization)), source=0_pInt) allocate(thermal_typeInstance(size(config_homogenization)), source=0_pInt) allocate(damage_typeInstance(size(config_homogenization)), source=0_pInt) - allocate(vacancyflux_typeInstance(size(config_homogenization)), source=0_pInt) - allocate(porosity_typeInstance(size(config_homogenization)), source=0_pInt) - allocate(hydrogenflux_typeInstance(size(config_homogenization)), source=0_pInt) allocate(homogenization_Ngrains(size(config_homogenization)), source=0_pInt) allocate(homogenization_Noutput(size(config_homogenization)), source=0_pInt) allocate(homogenization_active(size(config_homogenization)), source=.false.) !!!!!!!!!!!!!!! allocate(thermal_initialT(size(config_homogenization)), source=300.0_pReal) allocate(damage_initialPhi(size(config_homogenization)), source=1.0_pReal) - allocate(vacancyflux_initialCv(size(config_homogenization)), source=0.0_pReal) - allocate(porosity_initialPhi(size(config_homogenization)), source=1.0_pReal) - allocate(hydrogenflux_initialCh(size(config_homogenization)), source=0.0_pReal) forall (h = 1_pInt:size(config_homogenization)) & homogenization_active(h) = any(mesh_homogenizationAt == h) @@ -620,53 +522,6 @@ subroutine material_parseHomogenization end select endif - - if (config_homogenization(h)%keyExists('vacancyflux')) then - vacancyflux_initialCv(h) = config_homogenization(h)%getFloat('cv0',defaultVal=0.0_pReal) - - tag = config_homogenization(h)%getString('vacancyflux') - select case (trim(tag)) - case(VACANCYFLUX_isoconc_label) - vacancyflux_type(h) = VACANCYFLUX_isoconc_ID - case(VACANCYFLUX_isochempot_label) - vacancyflux_type(h) = VACANCYFLUX_isochempot_ID - case(VACANCYFLUX_cahnhilliard_label) - vacancyflux_type(h) = VACANCYFLUX_cahnhilliard_ID - case default - call IO_error(500_pInt,ext_msg=trim(tag)) - end select - - endif - - if (config_homogenization(h)%keyExists('porosity')) then - !ToDo? - - tag = config_homogenization(h)%getString('porosity') - select case (trim(tag)) - case(POROSITY_NONE_label) - porosity_type(h) = POROSITY_none_ID - case(POROSITY_phasefield_label) - porosity_type(h) = POROSITY_phasefield_ID - case default - call IO_error(500_pInt,ext_msg=trim(tag)) - end select - - endif - - if (config_homogenization(h)%keyExists('hydrogenflux')) then - hydrogenflux_initialCh(h) = config_homogenization(h)%getFloat('ch0',defaultVal=0.0_pReal) - - tag = config_homogenization(h)%getString('hydrogenflux') - select case (trim(tag)) - case(HYDROGENFLUX_isoconc_label) - hydrogenflux_type(h) = HYDROGENFLUX_isoconc_ID - case(HYDROGENFLUX_cahnhilliard_label) - hydrogenflux_type(h) = HYDROGENFLUX_cahnhilliard_ID - case default - call IO_error(500_pInt,ext_msg=trim(tag)) - end select - - endif enddo @@ -674,9 +529,6 @@ subroutine material_parseHomogenization homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) - vacancyflux_typeInstance(h) = count(vacancyflux_type (1:h) == vacancyflux_type (h)) - porosity_typeInstance(h) = count(porosity_type (1:h) == porosity_type (h)) - hydrogenflux_typeInstance(h) = count(hydrogenflux_type (1:h) == hydrogenflux_type (h)) enddo homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active) @@ -866,12 +718,6 @@ subroutine material_parsePhase phase_source(sourceCtr,p) = SOURCE_damage_anisoBrittle_ID case (SOURCE_damage_anisoDuctile_label) phase_source(sourceCtr,p) = SOURCE_damage_anisoDuctile_ID - case (SOURCE_vacancy_phenoplasticity_label) - phase_source(sourceCtr,p) = SOURCE_vacancy_phenoplasticity_ID - case (SOURCE_vacancy_irradiation_label) - phase_source(sourceCtr,p) = SOURCE_vacancy_irradiation_ID - case (SOURCE_vacancy_thermalfluc_label) - phase_source(sourceCtr,p) = SOURCE_vacancy_thermalfluc_ID end select enddo @@ -890,10 +736,6 @@ subroutine material_parsePhase phase_kinematics(kinematicsCtr,p) = KINEMATICS_slipplane_opening_ID case (KINEMATICS_thermal_expansion_label) phase_kinematics(kinematicsCtr,p) = KINEMATICS_thermal_expansion_ID - case (KINEMATICS_vacancy_strain_label) - phase_kinematics(kinematicsCtr,p) = KINEMATICS_vacancy_strain_ID - case (KINEMATICS_hydrogen_strain_label) - phase_kinematics(kinematicsCtr,p) = KINEMATICS_hydrogen_strain_ID end select enddo #if defined(__GFORTRAN__) @@ -907,8 +749,6 @@ subroutine material_parsePhase select case (trim(str(stiffDegradationCtr))) case (STIFFNESS_DEGRADATION_damage_label) phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID - case (STIFFNESS_DEGRADATION_porosity_label) - phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_porosity_ID end select enddo enddo diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index d65fe583f..a44681676 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -1,8 +1,9 @@ !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for isotropic (ISOTROPIC) plasticity -!> @details Isotropic (ISOTROPIC) Plasticity which resembles the phenopowerlaw plasticity without +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine for isotropic plasticity +!> @details Isotropic Plasticity which resembles the phenopowerlaw plasticity without !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- @@ -10,57 +11,58 @@ module plastic_isotropic use prec, only: & pReal,& pInt - + implicit none private integer(pInt), dimension(:,:), allocatable, target, public :: & - plastic_isotropic_sizePostResult !< size of each post result output + plastic_isotropic_sizePostResult !< size of each post result output character(len=64), dimension(:,:), allocatable, target, public :: & - plastic_isotropic_output !< name of each post result output - integer(pInt), dimension(:), allocatable, target, public :: & - plastic_isotropic_Noutput !< number of outputs per instance - - enum, bind(c) - enumerator :: undefined_ID, & - flowstress_ID, & - strainrate_ID + plastic_isotropic_output !< name of each post result output + + enum, bind(c) + enumerator :: & + undefined_ID, & + flowstress_ID, & + strainrate_ID end enum - type, private :: tParameters !< container type for internal constitutive parameters - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID - real(pReal) :: & - fTaylor, & - tau0, & - gdot0, & - n, & + type, private :: tParameters + real(pReal) :: & + fTaylor, & !< Taylor factor + tau0, & !< initial critical stress + gdot0, & !< reference strain rate + n, & !< stress exponent h0, & h0_slopeLnRate, & - tausat, & + tausat, & !< maximum critical stress a, & - aTolFlowstress, & - aTolShear, & tausat_SinhFitA, & tausat_SinhFitB, & tausat_SinhFitC, & - tausat_SinhFitD - logical :: & + tausat_SinhFitD, & + aTolFlowstress, & + aTolShear + integer(pInt) :: & + of_debug = 0_pInt + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID + logical :: & dilatation end type - type(tParameters), dimension(:), allocatable, target, private :: param !< containers of constitutive parameters (len Ninstance) - - type, private :: tIsotropicState !< internal state aliases - real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance + type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + + type, private :: tIsotropicState + real(pReal), pointer, dimension(:) :: & flowstress, & accumulatedShear end type - type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance - state, & - dotState + type(tIsotropicState), allocatable, dimension(:), private :: & + dotState, & + state - public :: & + public :: & plastic_isotropic_init, & plastic_isotropic_LpAndItsTangent, & plastic_isotropic_LiAndItsTangent, & @@ -80,20 +82,29 @@ subroutine plastic_isotropic_init() compiler_version, & compiler_options #endif -use IO + use prec, only: & + pStringLen use debug, only: & +#ifdef DEBUG + debug_e, & + debug_i, & + debug_g, & + debug_levelExtensive, & +#endif debug_level, & - debug_constitutive, & + debug_constitutive,& debug_levelBasic - use numerics, only: & - numerics_integrator - use math, only: & - math_Mandel3333to66, & - math_Voigt66to3333 + use IO, only: & + IO_error, & + IO_timeStamp use material, only: & +#ifdef DEBUG + phasememberAt, & +#endif phase_plasticity, & phase_plasticityInstance, & phase_Noutput, & + material_allocatePlasticState, & PLASTICITY_ISOTROPIC_label, & PLASTICITY_ISOTROPIC_ID, & material_phase, & @@ -101,148 +112,140 @@ use IO use config, only: & MATERIAL_partPhase, & config_phase - - use lattice + use lattice implicit none - - type(tParameters), pointer :: prm - integer(pInt) :: & - phase, & - instance, & - maxNinstance, & - sizeDotState, & - sizeState, & - sizeDeltaState - character(len=65536) :: & - extmsg = '' - integer(pInt) :: NipcMyPhase,i - character(len=65536), dimension(:), allocatable :: outputs + Ninstance, & + p, i, & + NipcMyPhase, & + sizeState, sizeDotState + + character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] + + integer(kind(undefined_ID)) :: & + outputID - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' + character(len=pStringLen) :: & + extmsg = '' + character(len=65536), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_ISOTROPIC_label//' init -+>>>' + write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia, 145:37-40, 2018' + write(6,'(/,a)') ' https://doi.org/10.1016/j.scriptamat.2017.09.047' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - - maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt) + + Ninstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance -! public variables - allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt) - allocate(plastic_isotropic_output(maxval(phase_Noutput), maxNinstance)) + allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), Ninstance),source=0_pInt) + allocate(plastic_isotropic_output(maxval(phase_Noutput), Ninstance)) plastic_isotropic_output = '' - allocate(plastic_isotropic_Noutput(maxNinstance), source=0_pInt) - allocate(param(maxNinstance)) ! one container of parameters per instance - allocate(state(maxNinstance)) ! internal state aliases - allocate(dotState(maxNinstance)) + allocate(param(Ninstance)) + allocate(state(Ninstance)) + allocate(dotState(Ninstance)) - do phase = 1_pInt, size(phase_plasticityInstance) - if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then - instance = phase_plasticityInstance(phase) - prm => param(instance) ! shorthand pointer to parameter object of my constitutive law - prm%tau0 = config_phase(phase)%getFloat('tau0') - prm%tausat = config_phase(phase)%getFloat('tausat') - prm%gdot0 = config_phase(phase)%getFloat('gdot0') - prm%n = config_phase(phase)%getFloat('n') - prm%h0 = config_phase(phase)%getFloat('h0') - prm%fTaylor = config_phase(phase)%getFloat('m') - prm%h0_slopeLnRate = config_phase(phase)%getFloat('h0_slopelnrate', defaultVal=0.0_pReal) - prm%tausat_SinhFitA = config_phase(phase)%getFloat('tausat_sinhfita',defaultVal=0.0_pReal) - prm%tausat_SinhFitB = config_phase(phase)%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal) - prm%tausat_SinhFitC = config_phase(phase)%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal) - prm%tausat_SinhFitD = config_phase(phase)%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal) - prm%a = config_phase(phase)%getFloat('a') - prm%aTolFlowStress = config_phase(phase)%getFloat('atol_flowstress',defaultVal=1.0_pReal) - prm%aTolShear = config_phase(phase)%getFloat('atol_shear',defaultVal=1.0e-6_pReal) - - prm%dilatation = config_phase(phase)%keyExists('/dilatation/') - -#if defined(__GFORTRAN__) - outputs = ['GfortranBug86277'] - outputs = config_phase(phase)%getStrings('(output)',defaultVal=outputs) - if (outputs(1) == 'GfortranBug86277') outputs = [character(len=65536)::] -#else - outputs = config_phase(phase)%getStrings('(output)',defaultVal=[character(len=65536)::]) + do p = 1_pInt, size(phase_plasticityInstance) + if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)), & + config => config_phase(p)) + +#ifdef DEBUG + if (p==material_phase(debug_g,debug_i,debug_e)) then + prm%of_debug = phasememberAt(debug_g,debug_i,debug_e) + endif #endif - allocate(prm%outputID(0)) - do i=1_pInt, size(outputs) - select case(outputs(i)) - case ('flowstress') - plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt - plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputs(i) - plasticState(phase)%sizePostResults = plasticState(phase)%sizePostResults + 1_pInt - plastic_isotropic_sizePostResult(i,instance) = 1_pInt - prm%outputID = [prm%outputID,flowstress_ID] - case ('strainrate') - plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt - plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputs(i) - plasticState(phase)%sizePostResults = & - plasticState(phase)%sizePostResults + 1_pInt - plastic_isotropic_sizePostResult(i,instance) = 1_pInt - prm%outputID = [prm%outputID,strainrate_ID] - end select - enddo + + prm%tau0 = config%getFloat('tau0') + prm%tausat = config%getFloat('tausat') + prm%gdot0 = config%getFloat('gdot0') + prm%n = config%getFloat('n') + prm%h0 = config%getFloat('h0') + prm%fTaylor = config%getFloat('m') + prm%h0_slopeLnRate = config%getFloat('h0_slopelnrate', defaultVal=0.0_pReal) + prm%tausat_SinhFitA = config%getFloat('tausat_sinhfita',defaultVal=0.0_pReal) + prm%tausat_SinhFitB = config%getFloat('tausat_sinhfitb',defaultVal=0.0_pReal) + prm%tausat_SinhFitC = config%getFloat('tausat_sinhfitc',defaultVal=0.0_pReal) + prm%tausat_SinhFitD = config%getFloat('tausat_sinhfitd',defaultVal=0.0_pReal) + prm%a = config%getFloat('a') + prm%aTolFlowStress = config%getFloat('atol_flowstress',defaultVal=1.0_pReal) + prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + + prm%dilatation = config%keyExists('/dilatation/') !-------------------------------------------------------------------------------------------------- ! sanity checks - extmsg = '' - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//"'aTolShear' " - if (prm%tau0 < 0.0_pReal) extmsg = trim(extmsg)//"'tau0' " - if (prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//"'gdot0' " - if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//"'n' " - if (prm%tausat <= prm%tau0) extmsg = trim(extmsg)//"'tausat' " - if (prm%a <= 0.0_pReal) extmsg = trim(extmsg)//"'a' " - if (prm%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//"'m' " - if (prm%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//"'atol_flowstress' " - if (extmsg /= '') call IO_error(211_pInt,ip=instance,& - ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')') + extmsg = '' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'aTolShear ' + if (prm%tau0 < 0.0_pReal) extmsg = trim(extmsg)//'tau0 ' + if (prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0 ' + if (prm%n <= 0.0_pReal) extmsg = trim(extmsg)//'n ' + if (prm%tausat <= prm%tau0) extmsg = trim(extmsg)//'tausat ' + if (prm%a <= 0.0_pReal) extmsg = trim(extmsg)//'a ' + if (prm%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//'m ' + if (prm%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//'atol_flowstress ' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'atol_shear ' + +!-------------------------------------------------------------------------------------------------- +! exit if any parameter is out of range + if (extmsg /= '') & + call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')') + +!-------------------------------------------------------------------------------------------------- +! output pararameters + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1_pInt, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + case ('flowstress') + outputID = flowstress_ID + case ('strainrate') + outputID = strainrate_ID + end select + + if (outputID /= undefined_ID) then + plastic_isotropic_output(i,phase_plasticityInstance(p)) = outputs(i) + plastic_isotropic_sizePostResult(i,phase_plasticityInstance(p)) = 1_pInt + prm%outputID = [prm%outputID , outputID] + endif + + end do !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) + NipcMyPhase = count(material_phase == p) + sizeState = size(["flowstress ","accumulated_shear"]) + sizeDotState = sizeState - sizeDotState = size(["flowstress ","accumulated_shear"]) - sizeDeltaState = 0_pInt ! no sudden jumps in state - sizeState = sizeDotState + sizeDeltaState - plasticState(phase)%sizeState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%nSlip = 1 - allocate(plasticState(phase)%aTolState ( sizeState)) - allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%state ( sizeState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & + 1_pInt,0_pInt,0_pInt) + plasticState(p)%sizePostResults = sum(plastic_isotropic_sizePostResult(:,phase_plasticityInstance(p))) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState + stt%flowstress => plasticState(p)%state (1,1:NipcMyPhase) + stt%flowstress = prm%tau0 + dot%flowstress => plasticState(p)%dotState (1,1:NipcMyPhase) + plasticState(p)%aTolState(1) = prm%aTolFlowstress - state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase) - dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase) - plasticState(phase)%state0(1,1:NipcMyPhase) = prm%tau0 - plasticState(phase)%aTolState(1) = prm%aTolFlowstress + stt%accumulatedShear => plasticState(p)%state (2,1:NipcMyPhase) + dot%accumulatedShear => plasticState(p)%dotState (2,1:NipcMyPhase) + plasticState(p)%aTolState(2) = prm%aTolShear + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(2:2,1:NipcMyPhase) + plasticState(p)%accumulatedSlip => plasticState(p)%state (2:2,1:NipcMyPhase) + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate - state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase) - dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase) - plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal - plasticState(phase)%aTolState(2) = prm%aTolShear - ! global alias - plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase) - plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase) - -endif enddo end subroutine plastic_isotropic_init @@ -251,308 +254,234 @@ end subroutine plastic_isotropic_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) +subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +#ifdef DEBUG use debug, only: & debug_level, & - debug_constitutive, & - debug_levelBasic, & + debug_constitutive,& debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_g + debug_levelSelective +#endif use math, only: & - math_mul6x6, & - math_Mandel6to33, & - math_Plain3333to99, & math_deviatoric33, & math_mul33xx33 - use material, only: & - phasememberAt, & - material_phase, & - phase_plasticityInstance implicit none - real(pReal), dimension(3,3), intent(out) :: & + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & - dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + instance, & + of - type(tParameters), pointer :: prm - real(pReal), dimension(3,3) :: & - Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor - real(pReal), dimension(3,3,3,3) :: & - dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor + Mp_dev !< deviatoric part of the Mandel stress real(pReal) :: & gamma_dot, & !< strainrate - norm_Tstar_dev, & !< euclidean norm of Tstar_dev - squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev + norm_Mp_dev, & !< norm of the deviatoric part of the Mandel stress + squarenorm_Mp_dev !< square of the norm of the deviatoric part of the Mandel stress integer(pInt) :: & - instance, of, & k, l, m, n - of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - prm => param(instance) + associate(prm => param(instance), stt => state(instance)) - Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress - squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) - norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) + Mp_dev = math_deviatoric33(Mp) + squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev) + norm_Mp_dev = sqrt(squarenorm_Mp_dev) - if (norm_Tstar_dev <= 0.0_pReal) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero - Lp = 0.0_pReal - dLp_dTstar99 = 0.0_pReal - else - gamma_dot = prm%gdot0 & - * ( sqrt(1.5_pReal) * norm_Tstar_dev / prm%fTaylor / state(instance)%flowstress(of) ) & - **prm%n - - Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/prm%fTaylor + if (norm_Mp_dev > 0.0_pReal) then + gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%fTaylor*stt%flowstress(of))) **prm%n + Lp = Mp_dev/norm_Mp_dev * gamma_dot/prm%fTaylor +#ifdef DEBUG if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then - write(6,'(a,i8,1x,i2,1x,i3)') '<< CONST isotropic >> at el ip g ',el,ip,ipc + .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Tstar (dev) / MPa', & - transpose(Tstar_dev_33(1:3,1:3))*1.0e-6_pReal - write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Tstar_dev*1.0e-6_pReal + transpose(Mp_dev)*1.0e-6_pReal + write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Mp_dev*1.0e-6_pReal write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', gamma_dot end if -!-------------------------------------------------------------------------------------------------- -! Calculation of the tangent of Lp +#endif forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar_3333(k,l,m,n) = (prm%n-1.0_pReal) * & - Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev + dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal + dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) & - dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal - dLp_dTstar99 = math_Plain3333to99(gamma_dot / prm%fTaylor * & - dLp_dTstar_3333 / norm_Tstar_dev) + dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal + dLp_dMp = gamma_dot / prm%fTaylor * dLp_dMp / norm_Mp_dev + else + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal end if + + end associate + end subroutine plastic_isotropic_LpAndItsTangent + !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent +! ToDo: Rename Tstar to Mi? !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,el) +subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of) use math, only: & - math_mul6x6, & - math_Mandel6to33, & - math_Plain3333to99, & math_spherical33, & math_mul33xx33 - use material, only: & - phasememberAt, & - material_phase, & - phase_plasticityInstance implicit none real(pReal), dimension(3,3), intent(out) :: & - Li !< plastic velocity gradient + Li !< inleastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & - dLi_dTstar_3333 !< derivative of Li with respect to Tstar as 4th order tensor - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - - type(tParameters), pointer :: prm + dLi_dTstar !< derivative of Li with respect to the Mandel stress + real(pReal), dimension(3,3), intent(in) :: & + Tstar !< Mandel stress ToDo: Mi? + integer(pInt), intent(in) :: & + instance, & + of + real(pReal), dimension(3,3) :: & - Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor + Tstar_sph !< sphiatoric part of the Mandel stress real(pReal) :: & gamma_dot, & !< strainrate norm_Tstar_sph, & !< euclidean norm of Tstar_sph squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph integer(pInt) :: & - instance, of, & k, l, m, n - of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - prm => param(instance) + associate(prm => param(instance), stt => state(instance)) - Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress - squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33) + Tstar_sph = math_spherical33(Tstar) + squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph,Tstar_sph) norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) - if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero - gamma_dot = prm%gdot0 & - * (sqrt(1.5_pReal) * norm_Tstar_sph / prm%fTaylor / state(instance)%flowstress(of) ) & - **prm%n + if (prm%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! no stress or J2 plastitiy --> Li and its derivative are zero + gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Tstar_sph /(prm%fTaylor*stt%flowstress(of))) **prm%n - Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/prm%fTaylor - - !-------------------------------------------------------------------------------------------------- - ! Calculation of the tangent of Li + Li = Tstar_sph/norm_Tstar_sph * gamma_dot/prm%fTaylor forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLi_dTstar_3333(k,l,m,n) = (prm%n-1.0_pReal) * & - Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph + dLi_dTstar(k,l,m,n) = (prm%n-1.0_pReal) * Tstar_sph(k,l)*Tstar_sph(m,n) / squarenorm_Tstar_sph forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal + dLi_dTstar(k,l,k,l) = dLi_dTstar(k,l,k,l) + 1.0_pReal - dLi_dTstar_3333 = gamma_dot / prm%fTaylor * & - dLi_dTstar_3333 / norm_Tstar_sph + dLi_dTstar = gamma_dot / prm%fTaylor * dLi_dTstar / norm_Tstar_sph else - Li = 0.0_pReal - dLi_dTstar_3333 = 0.0_pReal + Li = 0.0_pReal + dLi_dTstar = 0.0_pReal endif + + end associate + end subroutine plastic_isotropic_LiAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el) +subroutine plastic_isotropic_dotState(Mp,instance,of) use prec, only: & dEq0 use math, only: & - math_mul6x6 - use material, only: & - phasememberAt, & - material_phase, & - phase_plasticityInstance - + math_mul33xx33, & + math_deviatoric33 + implicit none - real(pReal), dimension(6), intent(in):: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - type(tParameters), pointer :: prm - real(pReal), dimension(6) :: & - Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer(pInt), intent(in) :: & + instance, & + of + real(pReal) :: & gamma_dot, & !< strainrate hardening, & !< hardening coefficient saturation, & !< saturation flowstress - norm_Tstar_v !< euclidean norm of Tstar_dev - integer(pInt) :: & - instance, & !< instance of my instance (unique number of my constitutive model) - of !< shortcut notation for offset position in state array + norm_Mp !< norm of the (deviatoric) Mandel stress - of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - prm => param(instance) + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) -!-------------------------------------------------------------------------------------------------- -! norm of (deviatoric) 2nd Piola-Kirchhoff stress if (prm%dilatation) then - norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) + norm_Mp = sqrt(math_mul33xx33(Mp,Mp)) else - Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal - Tstar_dev_v(4:6) = Tstar_v(4:6) - norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) - end if -!-------------------------------------------------------------------------------------------------- -! strain rate - gamma_dot = prm%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & - / &!----------------------------------------------------------------------------------- - (prm%fTaylor*state(instance)%flowstress(of) ))**prm%n + norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp))) + endif + + gamma_dot = prm%gdot0 * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor*stt%flowstress(of))) **prm%n -!-------------------------------------------------------------------------------------------------- -! hardening coefficient if (abs(gamma_dot) > 1e-12_pReal) then if (dEq0(prm%tausat_SinhFitA)) then saturation = prm%tausat else saturation = prm%tausat & - + asinh( (gamma_dot / prm%tausat_SinhFitA& - )**(1.0_pReal / prm%tausat_SinhFitD)& + + asinh( (gamma_dot / prm%tausat_SinhFitA)**(1.0_pReal / prm%tausat_SinhFitD) & )**(1.0_pReal / prm%tausat_SinhFitC) & - / ( prm%tausat_SinhFitB & - * (gamma_dot / prm%gdot0)**(1.0_pReal / prm%n) & - ) + / prm%tausat_SinhFitB * (gamma_dot / prm%gdot0)**(1.0_pReal / prm%n) endif hardening = ( prm%h0 + prm%h0_slopeLnRate * log(gamma_dot) ) & - * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**prm%a & - * sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation) + * abs( 1.0_pReal - stt%flowstress(of)/saturation )**prm%a & + * sign(1.0_pReal, 1.0_pReal - stt%flowstress(of)/saturation) else hardening = 0.0_pReal endif - dotState(instance)%flowstress (of) = hardening * gamma_dot - dotState(instance)%accumulatedShear(of) = gamma_dot + dot%flowstress (of) = hardening * gamma_dot + dot%accumulatedShear(of) = gamma_dot + + end associate end subroutine plastic_isotropic_dotState + !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function plastic_isotropic_postResults(Tstar_v,ipc,ip,el) +function plastic_isotropic_postResults(Mp,instance,of) result(postResults) use math, only: & - math_mul6x6 - use material, only: & - plasticState, & - material_phase, & - phasememberAt, & - phase_plasticityInstance + math_mul33xx33, & + math_deviatoric33 implicit none - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer(pInt), intent(in) :: & + instance, & + of - type(tParameters), pointer :: prm - - real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: & - plastic_isotropic_postResults + real(pReal), dimension(sum(plastic_isotropic_sizePostResult(:,instance))) :: & + postResults - real(pReal), dimension(6) :: & - Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal) :: & - norm_Tstar_v ! euclidean norm of Tstar_dev + norm_Mp !< norm of the Mandel stress integer(pInt) :: & - instance, & !< instance of my instance (unique number of my constitutive model) - of, & !< shortcut notation for offset position in state array - c, & - o + o,c - of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember - instance = phase_plasticityInstance(material_phase(ipc,ip,el)) - prm => param(instance) + associate(prm => param(instance), stt => state(instance)) -!-------------------------------------------------------------------------------------------------- -! norm of (deviatoric) 2nd Piola-Kirchhoff stress if (prm%dilatation) then - norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) + norm_Mp = sqrt(math_mul33xx33(Mp,Mp)) else - Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal - Tstar_dev_v(4:6) = Tstar_v(4:6) - norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) - end if + norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp))) + endif c = 0_pInt - plastic_isotropic_postResults = 0.0_pReal - outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) + outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) case (flowstress_ID) - plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of) + postResults(c+1_pInt) = stt%flowstress(of) c = c + 1_pInt case (strainrate_ID) - plastic_isotropic_postResults(c+1_pInt) = & - prm%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & - / &!---------------------------------------------------------------------------------- - (prm%fTaylor * state(instance)%flowstress(of)) ) ** prm%n + postResults(c+1_pInt) = prm%gdot0 & + * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor * stt%flowstress(of)))**prm%n c = c + 1_pInt end select enddo outputsLoop + + end associate end function plastic_isotropic_postResults diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 5470c4a43..5b5bb49d1 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -48,7 +48,7 @@ subroutine plastic_none_init phase, & NofMyPhase - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 053fe958b..1e42876f9 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -82,8 +82,7 @@ module plastic_phenopowerlaw xi_slip, & xi_twin, & gamma_slip, & - gamma_twin, & - whole + gamma_twin end type type(tPhenopowerlawState), allocatable, dimension(:), private :: & @@ -95,6 +94,9 @@ module plastic_phenopowerlaw plastic_phenopowerlaw_LpAndItsTangent, & plastic_phenopowerlaw_dotState, & plastic_phenopowerlaw_postResults + private :: & + kinetics_slip, & + kinetics_twin contains @@ -110,8 +112,7 @@ subroutine plastic_phenopowerlaw_init compiler_options #endif use prec, only: & - pStringLen, & - dEq0 + pStringLen use debug, only: & debug_level, & debug_constitutive,& @@ -119,7 +120,6 @@ subroutine plastic_phenopowerlaw_init use math, only: & math_expand use IO, only: & - IO_warning, & IO_error, & IO_timeStamp use material, only: & @@ -149,7 +149,7 @@ subroutine plastic_phenopowerlaw_init character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] integer(kind(undefined_ID)) :: & - outputID !< ID of each post result output + outputID character(len=pStringLen) :: & structure = '',& @@ -157,7 +157,7 @@ subroutine plastic_phenopowerlaw_init character(len=65536), dimension(:), allocatable :: & outputs - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -177,20 +177,21 @@ subroutine plastic_phenopowerlaw_init if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p))) + stt => state(phase_plasticityInstance(p)), & + config => config_phase(p)) - structure = config_phase(p)%getString('lattice_structure') + structure = config%getString('lattice_structure') !-------------------------------------------------------------------------------------------------- ! optional parameters that need to be defined - prm%twinB = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal) - prm%twinC = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal) - prm%twinD = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal) - prm%twinE = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal) + prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) + prm%twinC = config%getFloat('twin_c',defaultVal=0.0_pReal) + prm%twinD = config%getFloat('twin_d',defaultVal=0.0_pReal) + prm%twinE = config%getFloat('twin_e',defaultVal=0.0_pReal) - prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config_phase(p)%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - prm%aTolTwinfrac = config_phase(p)%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) ! sanity checks if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//'aTolresistance ' @@ -199,14 +200,14 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) slipActive: if (prm%totalNslip > 0_pInt) then prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& - config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) + config%getFloat('c/a',defaultVal=0.0_pReal)) if(structure=='bcc') then - prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray) + prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& + defaultVal = emptyRealArray) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) else @@ -214,18 +215,18 @@ subroutine plastic_phenopowerlaw_init prm%nonSchmid_neg = prm%Schmid_slip endif prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & - config_phase(p)%getFloats('interaction_slipslip'), & + config%getFloats('interaction_slipslip'), & structure(1:3)) - prm%xi_slip_0 = config_phase(p)%getFloats('tau0_slip', requiredShape=shape(prm%Nslip)) - prm%xi_slip_sat = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip)) - prm%H_int = config_phase(p)%getFloats('h_int', requiredShape=shape(prm%Nslip), & - defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) + prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) + prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) + prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & + defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) - prm%gdot0_slip = config_phase(p)%getFloat('gdot0_slip') - prm%n_slip = config_phase(p)%getFloat('n_slip') - prm%a_slip = config_phase(p)%getFloat('a_slip') - prm%h0_SlipSlip = config_phase(p)%getFloat('h0_slipslip') + prm%gdot0_slip = config%getFloat('gdot0_slip') + prm%n_slip = config%getFloat('n_slip') + prm%a_slip = config%getFloat('a_slip') + prm%h0_SlipSlip = config%getFloat('h0_slipslip') ! expand: family => system prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) @@ -233,9 +234,9 @@ subroutine plastic_phenopowerlaw_init prm%H_int = math_expand(prm%H_int, prm%Nslip) ! sanity checks - if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_slip ' - if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//'a_slip ' ! ToDo: negative values ok? - if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//'n_slip ' ! ToDo: negative values ok? + if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_slip ' + if (prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//'a_slip ' + if (prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//'n_slip ' if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//'xi_slip_0 ' if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat ' else slipActive @@ -245,30 +246,30 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! twin related parameters - prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray) + prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) prm%totalNtwin = sum(prm%Ntwin) twinActive: if (prm%totalNtwin > 0_pInt) then prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),& - config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) + config%getFloat('c/a',defaultVal=0.0_pReal)) prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,& - config_phase(p)%getFloats('interaction_twintwin'), & + config%getFloats('interaction_twintwin'), & structure(1:3)) prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& - config_phase(p)%getFloat('c/a')) + config%getFloat('c/a')) - prm%xi_twin_0 = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin)) + prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) - prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') - prm%n_twin = config_phase(p)%getFloat('n_twin') - prm%spr = config_phase(p)%getFloat('s_pr') - prm%h0_TwinTwin = config_phase(p)%getFloat('h0_twintwin') + prm%gdot0_twin = config%getFloat('gdot0_twin') + prm%n_twin = config%getFloat('n_twin') + prm%spr = config%getFloat('s_pr') + prm%h0_TwinTwin = config%getFloat('h0_twintwin') ! expand: family => system prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) ! sanity checks if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_twin ' - if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//'n_twin ' ! ToDo: negative values ok? + if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//'n_twin ' else twinActive allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%xi_twin_0(0)) @@ -279,10 +280,10 @@ subroutine plastic_phenopowerlaw_init ! slip-twin related parameters slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& - config_phase(p)%getFloats('interaction_sliptwin'), & + config%getFloats('interaction_sliptwin'), & structure(1:3)) prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& - config_phase(p)%getFloats('interaction_twinslip'), & + config%getFloats('interaction_twinslip'), & structure(1:3)) else slipAndTwinActive allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 @@ -297,7 +298,7 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! output pararameters - outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray) + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) do i=1_pInt, size(outputs) outputID = undefined_ID @@ -340,7 +341,7 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase + NipcMyPhase = count(material_phase == p) sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin sizeDotState = sizeState @@ -349,7 +350,6 @@ subroutine plastic_phenopowerlaw_init prm%totalNslip,prm%totalNtwin,0_pInt) plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,phase_plasticityInstance(p))) - !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1_pInt @@ -381,8 +381,7 @@ subroutine plastic_phenopowerlaw_init dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - dot%whole => plasticState(p)%dotState + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally end associate enddo @@ -398,7 +397,7 @@ end subroutine plastic_phenopowerlaw_init subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) implicit none - real(pReal), dimension(3,3), intent(out) :: & + real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & dLp_dMp !< derivative of Lp with respect to the Mandel stress @@ -420,9 +419,9 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) Lp = 0.0_pReal dLp_dMp = 0.0_pReal - associate(prm => param(instance), stt => state(instance)) + associate(prm => param(instance)) - call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1_pInt, prm%totalNslip Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & @@ -431,7 +430,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin) + call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) twinSystems: do i = 1_pInt, prm%totalNtwin Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & @@ -452,7 +451,7 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) implicit none real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & instance, & of @@ -468,7 +467,6 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - dot%whole(:,of) = 0.0_pReal sumGamma = sum(stt%gamma_slip(:,of)) sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) @@ -487,9 +485,9 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) !-------------------------------------------------------------------------------------------------- ! shear rates - call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of)) + call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) !-------------------------------------------------------------------------------------------------- ! hardening @@ -510,40 +508,108 @@ end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details Shear rates are calculated only optionally. +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) + use math, only: & + math_mul33xx33 + + implicit none + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer(pInt), intent(in) :: & + instance, & + of + + real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: & + postResults + + integer(pInt) :: & + o,c,i + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos,gdot_slip_neg + + c = 0_pInt + + associate( prm => param(instance), stt => state(instance)) + + outputsLoop: do o = 1_pInt,size(prm%outputID) + select case(prm%outputID(o)) + + case (resistance_slip_ID) + postResults(c+1_pInt:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of) + c = c + prm%totalNslip + case (accumulatedshear_slip_ID) + postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of) + c = c + prm%totalNslip + case (shearrate_slip_ID) + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) + postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg + c = c + prm%totalNslip + case (resolvedstress_slip_ID) + do i = 1_pInt, prm%totalNslip + postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) + enddo + c = c + prm%totalNslip + + case (resistance_twin_ID) + postResults(c+1_pInt:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (accumulatedshear_twin_ID) + postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of) + c = c + prm%totalNtwin + case (shearrate_twin_ID) + call kinetics_twin(Mp,instance,of,postResults(c+1_pInt:c+prm%totalNtwin)) + c = c + prm%totalNtwin + case (resolvedstress_twin_ID) + do i = 1_pInt, prm%totalNtwin + postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) + enddo + c = c + prm%totalNtwin + + end select + enddo outputsLoop + + end associate + +end function plastic_phenopowerlaw_postResults + + +!-------------------------------------------------------------------------------------------------- +!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress +!> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & - dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) +pure subroutine kinetics_slip(Mp,instance,of, & + gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) use prec, only: & dNeq0 use math, only: & math_mul33xx33 implicit none - type(tParameters), intent(in) :: & - prm - type(tPhenopowerlawState), intent(in) :: & - stt - integer(pInt), intent(in) :: & + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer(pInt), intent(in) :: & + instance, & of - real(pReal), dimension(prm%totalNslip), intent(out) :: & + + real(pReal), dimension(param(instance)%totalNslip), intent(out) :: & gdot_slip_pos, & gdot_slip_neg - real(pReal), dimension(prm%totalNslip), optional, intent(out) :: & + real(pReal), dimension(param(instance)%totalNslip), intent(out), optional :: & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg - real(pReal), dimension(3,3), intent(in) :: & - Mp - real(pReal), dimension(prm%totalNslip) :: & + real(pReal), dimension(param(instance)%totalNslip) :: & tau_slip_pos, & tau_slip_neg integer(pInt) :: i logical :: nonSchmidActive + associate(prm => param(instance), stt => state(instance)) + nonSchmidActive = size(prm%nonSchmidCoeff) > 0_pInt do i = 1_pInt, prm%totalNslip @@ -580,40 +646,42 @@ pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & dgdot_dtau_slip_neg = 0.0_pReal end where endif + end associate end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress. +!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress. ! twinning is assumed to take place only in untwinned volume. !> @details Derivates are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) +pure subroutine kinetics_twin(Mp,instance,of,& + gdot_twin,dgdot_dtau_twin) use prec, only: & dNeq0 use math, only: & math_mul33xx33 implicit none - type(tParameters), intent(in) :: & - prm - type(tPhenopowerlawState), intent(in) :: & - stt - integer(pInt), intent(in) :: & + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer(pInt), intent(in) :: & + instance, & of - real(pReal), dimension(3,3), intent(in) :: & - Mp - real(pReal), dimension(prm%totalNtwin), intent(out) :: & + + real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & gdot_twin - real(pReal), dimension(prm%totalNtwin), optional, intent(out) :: & + real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & dgdot_dtau_twin - real(pReal), dimension(prm%totalNtwin) :: & + real(pReal), dimension(param(instance)%totalNtwin) :: & tau_twin integer(pInt) :: i + + associate(prm => param(instance), stt => state(instance)) do i = 1_pInt, prm%totalNtwin tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) @@ -633,76 +701,9 @@ pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) dgdot_dtau_twin = 0.0_pReal end where endif - -end subroutine kinetics_twin - - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of constitutive results -!-------------------------------------------------------------------------------------------------- -function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) - use math, only: & - math_mul33xx33 - - implicit none - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer(pInt), intent(in) :: & - instance, & - of - - real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: & - postResults - - integer(pInt) :: & - o,c,i - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos,gdot_slip_neg - - postResults = 0.0_pReal - c = 0_pInt - - associate( prm => param(instance), stt => state(instance)) - - outputsLoop: do o = 1_pInt,size(prm%outputID) - select case(prm%outputID(o)) - - case (resistance_slip_ID) - postResults(c+1_pInt:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of) - c = c + prm%totalNslip - case (accumulatedshear_slip_ID) - postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of) - c = c + prm%totalNslip - case (shearrate_slip_ID) - call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) - postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg - c = c + prm%totalNslip - case (resolvedstress_slip_ID) - do i = 1_pInt, prm%totalNslip - postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) - enddo - c = c + prm%totalNslip - - case (resistance_twin_ID) - postResults(c+1_pInt:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of) - c = c + prm%totalNtwin - case (accumulatedshear_twin_ID) - postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of) - c = c + prm%totalNtwin - case (shearrate_twin_ID) - call kinetics_twin(prm,stt,of,Mp,postResults(c+1_pInt:c+prm%totalNtwin)) - c = c + prm%totalNtwin - case (resolvedstress_twin_ID) - do i = 1_pInt, prm%totalNtwin - postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) - enddo - c = c + prm%totalNtwin - - end select - enddo outputsLoop end associate -end function plastic_phenopowerlaw_postResults +end subroutine kinetics_twin end module plastic_phenopowerlaw diff --git a/src/porosity_none.f90 b/src/porosity_none.f90 deleted file mode 100644 index d8175cd9e..000000000 --- a/src/porosity_none.f90 +++ /dev/null @@ -1,60 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for constant porosity -!-------------------------------------------------------------------------------------------------- -module porosity_none - - implicit none - private - - public :: & - porosity_none_init - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields, reads information from material configuration file -!-------------------------------------------------------------------------------------------------- -subroutine porosity_none_init() -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use prec, only: & - pReal, & - pInt - use IO, only: & - IO_timeStamp - use material - use config - - implicit none - integer(pInt) :: & - homog, & - NofMyHomog - - write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_none_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - initializeInstances: do homog = 1_pInt, material_Nhomogenization - - myhomog: if (porosity_type(homog) == POROSITY_none_ID) then - NofMyHomog = count(material_homog == homog) - porosityState(homog)%sizeState = 0_pInt - porosityState(homog)%sizePostResults = 0_pInt - allocate(porosityState(homog)%state0 (0_pInt,NofMyHomog), source=0.0_pReal) - allocate(porosityState(homog)%subState0(0_pInt,NofMyHomog), source=0.0_pReal) - allocate(porosityState(homog)%state (0_pInt,NofMyHomog), source=0.0_pReal) - - deallocate(porosity(homog)%p) - allocate (porosity(homog)%p(1), source=porosity_initialPhi(homog)) - - endif myhomog - enddo initializeInstances - - -end subroutine porosity_none_init - -end module porosity_none diff --git a/src/porosity_phasefield.f90 b/src/porosity_phasefield.f90 deleted file mode 100644 index 1975ba64c..000000000 --- a/src/porosity_phasefield.f90 +++ /dev/null @@ -1,448 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for phase field modelling of pore nucleation and growth -!> @details phase field model for pore nucleation and growth based on vacancy clustering -!-------------------------------------------------------------------------------------------------- -module porosity_phasefield - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - porosity_phasefield_sizePostResults !< cumulative size of post results - - integer(pInt), dimension(:,:), allocatable, target, public :: & - porosity_phasefield_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - porosity_phasefield_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - porosity_phasefield_Noutput !< number of outputs per instance of this porosity - - enum, bind(c) - enumerator :: undefined_ID, & - porosity_ID - end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - porosity_phasefield_outputID !< ID of each post result output - - - public :: & - porosity_phasefield_init, & - porosity_phasefield_getFormationEnergy, & - porosity_phasefield_getSurfaceEnergy, & - porosity_phasefield_getSourceAndItsTangent, & - porosity_phasefield_getDiffusion33, & - porosity_phasefield_getMobility, & - porosity_phasefield_putPorosity, & - porosity_phasefield_postResults - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine porosity_phasefield_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - porosity_type, & - porosity_typeInstance, & - homogenization_Noutput, & - POROSITY_phasefield_label, & - POROSITY_phasefield_ID, & - material_homog, & - mappingHomogenization, & - porosityState, & - porosityMapping, & - porosity, & - porosity_initialPhi - use config, only: & - material_partHomogenization, & - material_partPhase - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o - integer(pInt) :: sizeState - integer(pInt) :: NofMyHomog - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- porosity_'//POROSITY_phasefield_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(porosity_type == POROSITY_phasefield_ID),pInt) - if (maxNinstance == 0_pInt) return - - allocate(porosity_phasefield_sizePostResults(maxNinstance), source=0_pInt) - allocate(porosity_phasefield_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) - allocate(porosity_phasefield_output (maxval(homogenization_Noutput),maxNinstance)) - porosity_phasefield_output = '' - allocate(porosity_phasefield_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) - allocate(porosity_phasefield_Noutput (maxNinstance), source=0_pInt) - - rewind(fileUnit) - section = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to - line = IO_read(fileUnit) - enddo - - parsingHomog: do while (trim(line) /= IO_EOF) ! read through sections of homog part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next homog section - section = section + 1_pInt ! advance homog section counter - cycle ! skip to next line - endif - - if (section > 0_pInt ) then; if (porosity_type(section) == POROSITY_phasefield_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = porosity_typeInstance(section) ! which instance of my porosity is present homog - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('porosity') - porosity_phasefield_Noutput(instance) = porosity_phasefield_Noutput(instance) + 1_pInt - porosity_phasefield_outputID(porosity_phasefield_Noutput(instance),instance) = porosity_ID - porosity_phasefield_output(porosity_phasefield_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - end select - - end select - endif; endif - enddo parsingHomog - - initializeInstances: do section = 1_pInt, size(porosity_type) - if (porosity_type(section) == POROSITY_phasefield_ID) then - NofMyHomog=count(material_homog==section) - instance = porosity_typeInstance(section) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,porosity_phasefield_Noutput(instance) - select case(porosity_phasefield_outputID(o,instance)) - case(porosity_ID) - mySize = 1_pInt - end select - - if (mySize > 0_pInt) then ! any meaningful output found - porosity_phasefield_sizePostResult(o,instance) = mySize - porosity_phasefield_sizePostResults(instance) = porosity_phasefield_sizePostResults(instance) + mySize - endif - enddo outputsLoop - -! allocate state arrays - sizeState = 0_pInt - porosityState(section)%sizeState = sizeState - porosityState(section)%sizePostResults = porosity_phasefield_sizePostResults(instance) - allocate(porosityState(section)%state0 (sizeState,NofMyHomog)) - allocate(porosityState(section)%subState0(sizeState,NofMyHomog)) - allocate(porosityState(section)%state (sizeState,NofMyHomog)) - - nullify(porosityMapping(section)%p) - porosityMapping(section)%p => mappingHomogenization(1,:,:) - deallocate(porosity(section)%p) - allocate(porosity(section)%p(NofMyHomog), source=porosity_initialPhi(section)) - - endif - - enddo initializeInstances -end subroutine porosity_phasefield_init - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized vacancy formation energy -!-------------------------------------------------------------------------------------------------- -function porosity_phasefield_getFormationEnergy(ip,el) - use lattice, only: & - lattice_vacancyFormationEnergy, & - lattice_vacancyVol - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - porosity_phasefield_getFormationEnergy - integer(pInt) :: & - grain - - porosity_phasefield_getFormationEnergy = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - porosity_phasefield_getFormationEnergy = porosity_phasefield_getFormationEnergy + & - lattice_vacancyFormationEnergy(material_phase(grain,ip,el))/ & - lattice_vacancyVol(material_phase(grain,ip,el)) - enddo - - porosity_phasefield_getFormationEnergy = & - porosity_phasefield_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function porosity_phasefield_getFormationEnergy - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized pore surface energy (normalized by characteristic length) -!-------------------------------------------------------------------------------------------------- -function porosity_phasefield_getSurfaceEnergy(ip,el) - use lattice, only: & - lattice_vacancySurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal) :: & - porosity_phasefield_getSurfaceEnergy - integer(pInt) :: & - grain - - porosity_phasefield_getSurfaceEnergy = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - porosity_phasefield_getSurfaceEnergy = porosity_phasefield_getSurfaceEnergy + & - lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) - enddo - - porosity_phasefield_getSurfaceEnergy = & - porosity_phasefield_getSurfaceEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function porosity_phasefield_getSurfaceEnergy - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized local driving force for pore nucleation and growth -!-------------------------------------------------------------------------------------------------- -subroutine porosity_phasefield_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - use math, only : & - math_mul33x33, & - math_mul66x6, & - math_Mandel33to6, & - math_transpose33, & - math_I3 - use material, only: & - homogenization_Ngrains, & - material_homog, & - material_phase, & - phase_NstiffnessDegradations, & - phase_stiffnessDegradation, & - vacancyConc, & - vacancyfluxMapping, & - damage, & - damageMapping, & - STIFFNESS_DEGRADATION_damage_ID - use crystallite, only: & - crystallite_Fe - use constitutive, only: & - constitutive_homogenizedC - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - integer(pInt) :: & - phase, & - grain, & - homog, & - mech - real(pReal) :: & - phiDot, dPhiDot_dPhi, Cv, W_e, strain(6), C(6,6) - - homog = material_homog(ip,el) - Cv = vacancyConc(homog)%p(vacancyfluxMapping(homog)%p(ip,el)) - - W_e = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - phase = material_phase(grain,ip,el) - strain = math_Mandel33to6(math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,grain,ip,el)), & - crystallite_Fe(1:3,1:3,grain,ip,el)) - math_I3)/2.0_pReal - C = constitutive_homogenizedC(grain,ip,el) - do mech = 1_pInt, phase_NstiffnessDegradations(phase) - select case(phase_stiffnessDegradation(mech,phase)) - case (STIFFNESS_DEGRADATION_damage_ID) - C = damage(homog)%p(damageMapping(homog)%p(ip,el))* & - damage(homog)%p(damageMapping(homog)%p(ip,el))* & - C - - end select - enddo - W_e = W_e + sum(abs(strain*math_mul66x6(C,strain))) - enddo - W_e = W_e/real(homogenization_Ngrains(homog),pReal) - - phiDot = 2.0_pReal*(1.0_pReal - phi)*(1.0_pReal - Cv)*(1.0_pReal - Cv) - & - 2.0_pReal*phi*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ & - porosity_phasefield_getSurfaceEnergy (ip,el) - dPhiDot_dPhi = - 2.0_pReal*(1.0_pReal - Cv)*(1.0_pReal - Cv) & - - 2.0_pReal*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ & - porosity_phasefield_getSurfaceEnergy (ip,el) - -end subroutine porosity_phasefield_getSourceAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized nonlocal diffusion tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function porosity_phasefield_getDiffusion33(ip,el) - use lattice, only: & - lattice_PorosityDiffusion33 - use material, only: & - homogenization_Ngrains, & - material_phase, & - mappingHomogenization - use crystallite, only: & - crystallite_push33ToRef - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - porosity_phasefield_getDiffusion33 - integer(pInt) :: & - homog, & - grain - - homog = mappingHomogenization(2,ip,el) - porosity_phasefield_getDiffusion33 = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - porosity_phasefield_getDiffusion33 = porosity_phasefield_getDiffusion33 + & - crystallite_push33ToRef(grain,ip,el,lattice_PorosityDiffusion33(1:3,1:3,material_phase(grain,ip,el))) - enddo - - porosity_phasefield_getDiffusion33 = & - porosity_phasefield_getDiffusion33/real(homogenization_Ngrains(homog),pReal) - -end function porosity_phasefield_getDiffusion33 - -!-------------------------------------------------------------------------------------------------- -!> @brief Returns homogenized phase field mobility -!-------------------------------------------------------------------------------------------------- -real(pReal) function porosity_phasefield_getMobility(ip,el) - use mesh, only: & - mesh_element - use lattice, only: & - lattice_PorosityMobility - use material, only: & - material_phase, & - homogenization_Ngrains - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - integer(pInt) :: & - ipc - - porosity_phasefield_getMobility = 0.0_pReal - - do ipc = 1, homogenization_Ngrains(mesh_element(3,el)) - porosity_phasefield_getMobility = porosity_phasefield_getMobility & - + lattice_PorosityMobility(material_phase(ipc,ip,el)) - enddo - - porosity_phasefield_getMobility = & - porosity_phasefield_getMobility/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function porosity_phasefield_getMobility - -!-------------------------------------------------------------------------------------------------- -!> @brief updates porosity with solution from phasefield PDE -!-------------------------------------------------------------------------------------------------- -subroutine porosity_phasefield_putPorosity(phi,ip,el) - use material, only: & - material_homog, & - porosityMapping, & - porosity - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - integer(pInt) :: & - homog, & - offset - - homog = material_homog(ip,el) - offset = porosityMapping(homog)%p(ip,el) - porosity(homog)%p(offset) = phi - -end subroutine porosity_phasefield_putPorosity - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of porosity results -!-------------------------------------------------------------------------------------------------- -function porosity_phasefield_postResults(ip,el) - use material, only: & - mappingHomogenization, & - porosity_typeInstance, & - porosity - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(porosity_phasefield_sizePostResults(porosity_typeInstance(mappingHomogenization(2,ip,el)))) :: & - porosity_phasefield_postResults - - integer(pInt) :: & - instance, homog, offset, o, c - - homog = mappingHomogenization(2,ip,el) - offset = mappingHomogenization(1,ip,el) - instance = porosity_typeInstance(homog) - - c = 0_pInt - porosity_phasefield_postResults = 0.0_pReal - - do o = 1_pInt,porosity_phasefield_Noutput(instance) - select case(porosity_phasefield_outputID(o,instance)) - - case (porosity_ID) - porosity_phasefield_postResults(c+1_pInt) = porosity(homog)%p(offset) - c = c + 1 - end select - enddo -end function porosity_phasefield_postResults - -end module porosity_phasefield diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 041761afe..5aa3648f3 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -246,10 +246,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) sourceState, & material_homog, & phase_NstiffnessDegradations, & - phase_stiffnessDegradation, & - porosity, & - porosityMapping, & - STIFFNESS_DEGRADATION_porosity_ID + phase_stiffnessDegradation use math, only : & math_mul33x33, & math_mul66x6, & @@ -279,15 +276,7 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) instance = source_damage_isoBrittle_instance(phase) !< instance of damage_isoBrittle source sourceOffset = source_damage_isoBrittle_offset(phase) - stiffness = C - do mech = 1_pInt, phase_NstiffnessDegradations(phase) - select case(phase_stiffnessDegradation(mech,phase)) - case (STIFFNESS_DEGRADATION_porosity_ID) - stiffness = porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* & - porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* & - stiffness - end select - enddo + stiffness = C strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) strainenergy = 2.0_pReal*sum(strain*math_mul66x6(stiffness,strain))/ & diff --git a/src/source_vacancy_irradiation.f90 b/src/source_vacancy_irradiation.f90 deleted file mode 100644 index 67b4cabcf..000000000 --- a/src/source_vacancy_irradiation.f90 +++ /dev/null @@ -1,248 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for vacancy generation due to irradiation -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module source_vacancy_irradiation - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - source_vacancy_irradiation_sizePostResults, & !< cumulative size of post results - source_vacancy_irradiation_offset, & !< which source is my current damage mechanism? - source_vacancy_irradiation_instance !< instance of damage source mechanism - - integer(pInt), dimension(:,:), allocatable, target, public :: & - source_vacancy_irradiation_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - source_vacancy_irradiation_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - source_vacancy_irradiation_Noutput !< number of outputs per instance of this damage - - real(pReal), dimension(:), allocatable, private :: & - source_vacancy_irradiation_cascadeProb, & - source_vacancy_irradiation_cascadeVolume - - public :: & - source_vacancy_irradiation_init, & - source_vacancy_irradiation_deltaState, & - source_vacancy_irradiation_getRateAndItsTangent - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_irradiation_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_vacancy_irradiation_label, & - SOURCE_vacancy_irradiation_ID, & - material_phase, & - sourceState - use config, only: & - material_Nphase, & - MATERIAL_partPhase - use numerics,only: & - numerics_integrator - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_irradiation_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_source == SOURCE_vacancy_irradiation_ID),pInt) - if (maxNinstance == 0_pInt) return - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(source_vacancy_irradiation_offset(material_Nphase), source=0_pInt) - allocate(source_vacancy_irradiation_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - source_vacancy_irradiation_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_irradiation_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_vacancy_irradiation_ID) & - source_vacancy_irradiation_offset(phase) = source - enddo - enddo - - allocate(source_vacancy_irradiation_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_vacancy_irradiation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_vacancy_irradiation_output(maxval(phase_Noutput),maxNinstance)) - source_vacancy_irradiation_output = '' - allocate(source_vacancy_irradiation_Noutput(maxNinstance), source=0_pInt) - allocate(source_vacancy_irradiation_cascadeProb(maxNinstance), source=0.0_pReal) - allocate(source_vacancy_irradiation_cascadeVolume(maxNinstance), source=0.0_pReal) - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - - if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = source_vacancy_irradiation_instance(phase) ! which instance of my vacancy is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('irradiation_cascadeprobability') - source_vacancy_irradiation_cascadeProb(instance) = IO_floatValue(line,chunkPos,2_pInt) - - case ('irradiation_cascadevolume') - source_vacancy_irradiation_cascadeVolume(instance) = IO_floatValue(line,chunkPos,2_pInt) - - end select - endif; endif - enddo parsingFile - - initializeInstances: do phase = 1_pInt, material_Nphase - if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then - NofMyPhase=count(material_phase==phase) - instance = source_vacancy_irradiation_instance(phase) - sourceOffset = source_vacancy_irradiation_offset(phase) - - sizeDotState = 2_pInt - sizeDeltaState = 2_pInt - sizeState = 2_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState - sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState - sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_irradiation_sizePostResults(instance) - allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.1_pReal) - allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal) - - allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) - - endif - - enddo initializeInstances -end subroutine source_vacancy_irradiation_init - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_irradiation_deltaState(ipc, ip, el) - use material, only: & - phaseAt, phasememberAt, & - sourceState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - integer(pInt) :: & - phase, constituent, sourceOffset - real(pReal) :: & - randNo - - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) - sourceOffset = source_vacancy_irradiation_offset(phase) - - call random_number(randNo) - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - randNo - sourceState(phase)%p(sourceOffset)%state(1,constituent) - call random_number(randNo) - sourceState(phase)%p(sourceOffset)%deltaState(2,constituent) = & - randNo - sourceState(phase)%p(sourceOffset)%state(2,constituent) - -end subroutine source_vacancy_irradiation_deltaState - -!-------------------------------------------------------------------------------------------------- -!> @brief returns local vacancy generation rate -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_irradiation_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el) - use material, only: & - phaseAt, phasememberAt, & - sourceState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(out) :: & - CvDot, dCvDot_dCv - integer(pInt) :: & - instance, phase, constituent, sourceOffset - - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) - instance = source_vacancy_irradiation_instance(phase) - sourceOffset = source_vacancy_irradiation_offset(phase) - - CvDot = 0.0_pReal - dCvDot_dCv = 0.0_pReal - if (sourceState(phase)%p(sourceOffset)%state0(1,constituent) < source_vacancy_irradiation_cascadeProb(instance)) & - CvDot = sourceState(phase)%p(sourceOffset)%state0(2,constituent)*source_vacancy_irradiation_cascadeVolume(instance) - -end subroutine source_vacancy_irradiation_getRateAndItsTangent - -end module source_vacancy_irradiation diff --git a/src/source_vacancy_phenoplasticity.f90 b/src/source_vacancy_phenoplasticity.f90 deleted file mode 100644 index e20d8ec06..000000000 --- a/src/source_vacancy_phenoplasticity.f90 +++ /dev/null @@ -1,210 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for vacancy generation due to plasticity -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module source_vacancy_phenoplasticity - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - source_vacancy_phenoplasticity_sizePostResults, & !< cumulative size of post results - source_vacancy_phenoplasticity_offset, & !< which source is my current damage mechanism? - source_vacancy_phenoplasticity_instance !< instance of damage source mechanism - - integer(pInt), dimension(:,:), allocatable, target, public :: & - source_vacancy_phenoplasticity_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - source_vacancy_phenoplasticity_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - source_vacancy_phenoplasticity_Noutput !< number of outputs per instance of this damage - - real(pReal), dimension(:), allocatable, private :: & - source_vacancy_phenoplasticity_rateCoeff - - public :: & - source_vacancy_phenoplasticity_init, & - source_vacancy_phenoplasticity_getRateAndItsTangent - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_phenoplasticity_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_vacancy_phenoplasticity_label, & - SOURCE_vacancy_phenoplasticity_ID, & - material_phase, & - sourceState - use config, only: & - material_Nphase, & - MATERIAL_partPhase - use numerics,only: & - numerics_integrator - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_phenoplasticity_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_source == SOURCE_vacancy_phenoplasticity_ID),pInt) - if (maxNinstance == 0_pInt) return - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(source_vacancy_phenoplasticity_offset(material_Nphase), source=0_pInt) - allocate(source_vacancy_phenoplasticity_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - source_vacancy_phenoplasticity_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_phenoplasticity_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_vacancy_phenoplasticity_ID) & - source_vacancy_phenoplasticity_offset(phase) = source - enddo - enddo - - allocate(source_vacancy_phenoplasticity_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_vacancy_phenoplasticity_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_vacancy_phenoplasticity_output(maxval(phase_Noutput),maxNinstance)) - source_vacancy_phenoplasticity_output = '' - allocate(source_vacancy_phenoplasticity_Noutput(maxNinstance), source=0_pInt) - allocate(source_vacancy_phenoplasticity_rateCoeff(maxNinstance), source=0.0_pReal) - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - - if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = source_vacancy_phenoplasticity_instance(phase) ! which instance of my vacancy is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('phenoplasticity_ratecoeff') - source_vacancy_phenoplasticity_rateCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt) - - end select - endif; endif - enddo parsingFile - - initializeInstances: do phase = 1_pInt, material_Nphase - if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then - NofMyPhase=count(material_phase==phase) - instance = source_vacancy_phenoplasticity_instance(phase) - sourceOffset = source_vacancy_phenoplasticity_offset(phase) - - sizeDotState = 0_pInt - sizeDeltaState = 0_pInt - sizeState = 0_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState - sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState - sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_phenoplasticity_sizePostResults(instance) - allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal) - - allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) - - endif - - enddo initializeInstances -end subroutine source_vacancy_phenoplasticity_init - -!-------------------------------------------------------------------------------------------------- -!> @brief returns local vacancy generation rate -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_phenoplasticity_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el) - use material, only: & - phaseAt, phasememberAt, & - plasticState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(out) :: & - CvDot, dCvDot_dCv - integer(pInt) :: & - instance, phase, constituent - - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) - instance = source_vacancy_phenoplasticity_instance(phase) - - CvDot = & - source_vacancy_phenoplasticity_rateCoeff(instance)* & - sum(plasticState(phase)%slipRate(:,constituent)) - dCvDot_dCv = 0.0_pReal - -end subroutine source_vacancy_phenoplasticity_getRateAndItsTangent - -end module source_vacancy_phenoplasticity diff --git a/src/source_vacancy_thermalfluc.f90 b/src/source_vacancy_thermalfluc.f90 deleted file mode 100644 index cea52aa75..000000000 --- a/src/source_vacancy_thermalfluc.f90 +++ /dev/null @@ -1,250 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for vacancy generation due to thermal fluctuations -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module source_vacancy_thermalfluc - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - source_vacancy_thermalfluc_sizePostResults, & !< cumulative size of post results - source_vacancy_thermalfluc_offset, & !< which source is my current damage mechanism? - source_vacancy_thermalfluc_instance !< instance of damage source mechanism - - integer(pInt), dimension(:,:), allocatable, target, public :: & - source_vacancy_thermalfluc_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - source_vacancy_thermalfluc_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - source_vacancy_thermalfluc_Noutput !< number of outputs per instance of this damage - - real(pReal), dimension(:), allocatable, private :: & - source_vacancy_thermalfluc_amplitude, & - source_vacancy_thermalfluc_normVacancyEnergy - - public :: & - source_vacancy_thermalfluc_init, & - source_vacancy_thermalfluc_deltaState, & - source_vacancy_thermalfluc_getRateAndItsTangent - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_thermalfluc_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use debug, only: & - debug_level,& - debug_constitutive,& - debug_levelBasic - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use lattice, only: & - lattice_vacancyFormationEnergy - use material, only: & - phase_source, & - phase_Nsources, & - phase_Noutput, & - SOURCE_vacancy_thermalfluc_label, & - SOURCE_vacancy_thermalfluc_ID, & - material_phase, & - sourceState - use config, only: & - material_Nphase, & - MATERIAL_partPhase - use numerics,only: & - numerics_integrator - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset - integer(pInt) :: sizeState, sizeDotState, sizeDeltaState - integer(pInt) :: NofMyPhase - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- source_'//SOURCE_vacancy_thermalfluc_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(phase_source == SOURCE_vacancy_thermalfluc_ID),pInt) - if (maxNinstance == 0_pInt) return - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance - - allocate(source_vacancy_thermalfluc_offset(material_Nphase), source=0_pInt) - allocate(source_vacancy_thermalfluc_instance(material_Nphase), source=0_pInt) - do phase = 1, material_Nphase - source_vacancy_thermalfluc_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_thermalfluc_ID) - do source = 1, phase_Nsources(phase) - if (phase_source(source,phase) == source_vacancy_thermalfluc_ID) & - source_vacancy_thermalfluc_offset(phase) = source - enddo - enddo - - allocate(source_vacancy_thermalfluc_sizePostResults(maxNinstance), source=0_pInt) - allocate(source_vacancy_thermalfluc_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) - allocate(source_vacancy_thermalfluc_output(maxval(phase_Noutput),maxNinstance)) - source_vacancy_thermalfluc_output = '' - allocate(source_vacancy_thermalfluc_Noutput(maxNinstance), source=0_pInt) - allocate(source_vacancy_thermalfluc_amplitude(maxNinstance), source=0.0_pReal) - allocate(source_vacancy_thermalfluc_normVacancyEnergy(maxNinstance), source=0.0_pReal) - - rewind(fileUnit) - phase = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next phase section - phase = phase + 1_pInt ! advance phase section counter - cycle ! skip to next line - endif - - if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = source_vacancy_thermalfluc_instance(phase) ! which instance of my vacancy is present phase - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('thermalfluctuation_amplitude') - source_vacancy_thermalfluc_amplitude(instance) = IO_floatValue(line,chunkPos,2_pInt) - - end select - endif; endif - enddo parsingFile - - initializeInstances: do phase = 1_pInt, material_Nphase - if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then - NofMyPhase=count(material_phase==phase) - instance = source_vacancy_thermalfluc_instance(phase) - source_vacancy_thermalfluc_normVacancyEnergy(instance) = & - lattice_vacancyFormationEnergy(phase)/1.3806488e-23_pReal - sourceOffset = source_vacancy_thermalfluc_offset(phase) - - sizeDotState = 1_pInt - sizeDeltaState = 1_pInt - sizeState = 1_pInt - sourceState(phase)%p(sourceOffset)%sizeState = sizeState - sourceState(phase)%p(sourceOffset)%sizeDotState = sizeDotState - sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState - sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance) - allocate(sourceState(phase)%p(sourceOffset)%aTolState (sizeState), source=0.1_pReal) - allocate(sourceState(phase)%p(sourceOffset)%state0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%state (sizeState,NofMyPhase), source=0.0_pReal) - - allocate(sourceState(phase)%p(sourceOffset)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%deltaState (sizeDeltaState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 1_pInt)) then - allocate(sourceState(phase)%p(sourceOffset)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(sourceState(phase)%p(sourceOffset)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) - if (any(numerics_integrator == 5_pInt)) & - allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) - - endif - - enddo initializeInstances -end subroutine source_vacancy_thermalfluc_init - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_thermalfluc_deltaState(ipc, ip, el) - use material, only: & - phaseAt, phasememberAt, & - sourceState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - integer(pInt) :: & - phase, constituent, sourceOffset - real(pReal) :: & - randNo - - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) - sourceOffset = source_vacancy_thermalfluc_offset(phase) - - call random_number(randNo) - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - randNo - 0.5_pReal - sourceState(phase)%p(sourceOffset)%state(1,constituent) - -end subroutine source_vacancy_thermalfluc_deltaState - -!-------------------------------------------------------------------------------------------------- -!> @brief returns local vacancy generation rate -!-------------------------------------------------------------------------------------------------- -subroutine source_vacancy_thermalfluc_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el) - use material, only: & - phaseAt, phasememberAt, & - material_homog, & - temperature, & - thermalMapping, & - sourceState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), intent(out) :: & - CvDot, dCvDot_dCv - integer(pInt) :: & - instance, phase, constituent, sourceOffset - - phase = phaseAt(ipc,ip,el) - constituent = phasememberAt(ipc,ip,el) - instance = source_vacancy_thermalfluc_instance(phase) - sourceOffset = source_vacancy_thermalfluc_offset(phase) - - CvDot = source_vacancy_thermalfluc_amplitude(instance)* & - sourceState(phase)%p(sourceOffset)%state0(2,constituent)* & - exp(-source_vacancy_thermalfluc_normVacancyEnergy(instance)/ & - temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))) - dCvDot_dCv = 0.0_pReal - -end subroutine source_vacancy_thermalfluc_getRateAndItsTangent - -end module source_vacancy_thermalfluc diff --git a/src/vacancyflux_cahnhilliard.f90 b/src/vacancyflux_cahnhilliard.f90 deleted file mode 100644 index ae5bd1cbc..000000000 --- a/src/vacancyflux_cahnhilliard.f90 +++ /dev/null @@ -1,602 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for conservative transport of vacancy concentration field -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module vacancyflux_cahnhilliard - use prec, only: & - pReal, & - pInt, & - group_float - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - vacancyflux_cahnhilliard_sizePostResults !< cumulative size of post results - - integer(pInt), dimension(:,:), allocatable, target, public :: & - vacancyflux_cahnhilliard_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - vacancyflux_cahnhilliard_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - vacancyflux_cahnhilliard_Noutput !< number of outputs per instance of this damage - - real(pReal), dimension(:), allocatable, private :: & - vacancyflux_cahnhilliard_flucAmplitude - - type(group_float), dimension(:), allocatable, private :: & - vacancyflux_cahnhilliard_thermalFluc - - real(pReal), parameter, private :: & - kB = 1.3806488e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: undefined_ID, & - vacancyConc_ID - end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - vacancyflux_cahnhilliard_outputID !< ID of each post result output - - - public :: & - vacancyflux_cahnhilliard_init, & - vacancyflux_cahnhilliard_getSourceAndItsTangent, & - vacancyflux_cahnhilliard_getMobility33, & - vacancyflux_cahnhilliard_getDiffusion33, & - vacancyflux_cahnhilliard_getChemPotAndItsTangent, & - vacancyflux_cahnhilliard_putVacancyConcAndItsRate, & - vacancyflux_cahnhilliard_postResults - private :: & - vacancyflux_cahnhilliard_getFormationEnergy, & - vacancyflux_cahnhilliard_getEntropicCoeff, & - vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_cahnhilliard_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - vacancyflux_type, & - vacancyflux_typeInstance, & - homogenization_Noutput, & - VACANCYFLUX_cahnhilliard_label, & - VACANCYFLUX_cahnhilliard_ID, & - material_homog, & - mappingHomogenization, & - vacancyfluxState, & - vacancyfluxMapping, & - vacancyConc, & - vacancyConcRate, & - vacancyflux_initialCv - use config, only: & - material_partPhase, & - material_partHomogenization - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o,offset - integer(pInt) :: sizeState - integer(pInt) :: NofMyHomog - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_cahnhilliard_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_cahnhilliard_ID),pInt) - if (maxNinstance == 0_pInt) return - - allocate(vacancyflux_cahnhilliard_sizePostResults(maxNinstance), source=0_pInt) - allocate(vacancyflux_cahnhilliard_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) - allocate(vacancyflux_cahnhilliard_output (maxval(homogenization_Noutput),maxNinstance)) - vacancyflux_cahnhilliard_output = '' - allocate(vacancyflux_cahnhilliard_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) - allocate(vacancyflux_cahnhilliard_Noutput (maxNinstance), source=0_pInt) - - allocate(vacancyflux_cahnhilliard_flucAmplitude (maxNinstance)) - allocate(vacancyflux_cahnhilliard_thermalFluc (maxNinstance)) - - rewind(fileUnit) - section = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to - line = IO_read(fileUnit) - enddo - - parsingHomog: do while (trim(line) /= IO_EOF) ! read through sections of homog part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next homog section - section = section + 1_pInt ! advance homog section counter - cycle ! skip to next line - endif - - if (section > 0_pInt ) then; if (vacancyflux_type(section) == VACANCYFLUX_cahnhilliard_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = vacancyflux_typeInstance(section) ! which instance of my vacancyflux is present homog - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('vacancyconc') - vacancyflux_cahnhilliard_Noutput(instance) = vacancyflux_cahnhilliard_Noutput(instance) + 1_pInt - vacancyflux_cahnhilliard_outputID(vacancyflux_cahnhilliard_Noutput(instance),instance) = vacancyConc_ID - vacancyflux_cahnhilliard_output(vacancyflux_cahnhilliard_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - end select - - case ('vacancyflux_flucamplitude') - vacancyflux_cahnhilliard_flucAmplitude(instance) = IO_floatValue(line,chunkPos,2_pInt) - - end select - endif; endif - enddo parsingHomog - - initializeInstances: do section = 1_pInt, size(vacancyflux_type) - if (vacancyflux_type(section) == VACANCYFLUX_cahnhilliard_ID) then - NofMyHomog=count(material_homog==section) - instance = vacancyflux_typeInstance(section) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,vacancyflux_cahnhilliard_Noutput(instance) - select case(vacancyflux_cahnhilliard_outputID(o,instance)) - case(vacancyConc_ID) - mySize = 1_pInt - end select - - if (mySize > 0_pInt) then ! any meaningful output found - vacancyflux_cahnhilliard_sizePostResult(o,instance) = mySize - vacancyflux_cahnhilliard_sizePostResults(instance) = vacancyflux_cahnhilliard_sizePostResults(instance) + mySize - endif - enddo outputsLoop - -! allocate state arrays - sizeState = 0_pInt - vacancyfluxState(section)%sizeState = sizeState - vacancyfluxState(section)%sizePostResults = vacancyflux_cahnhilliard_sizePostResults(instance) - allocate(vacancyfluxState(section)%state0 (sizeState,NofMyHomog)) - allocate(vacancyfluxState(section)%subState0(sizeState,NofMyHomog)) - allocate(vacancyfluxState(section)%state (sizeState,NofMyHomog)) - - allocate(vacancyflux_cahnhilliard_thermalFluc(instance)%p(NofMyHomog)) - do offset = 1_pInt, NofMyHomog - call random_number(vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset)) - vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset) = & - 1.0_pReal - & - vacancyflux_cahnhilliard_flucAmplitude(instance)* & - (vacancyflux_cahnhilliard_thermalFluc(instance)%p(offset) - 0.5_pReal) - enddo - - nullify(vacancyfluxMapping(section)%p) - vacancyfluxMapping(section)%p => mappingHomogenization(1,:,:) - deallocate(vacancyConc (section)%p) - allocate (vacancyConc (section)%p(NofMyHomog), source=vacancyflux_initialCv(section)) - deallocate(vacancyConcRate(section)%p) - allocate (vacancyConcRate(section)%p(NofMyHomog), source=0.0_pReal) - - endif - - enddo initializeInstances - -end subroutine vacancyflux_cahnhilliard_init - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized vacancy driving forces -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el) - use material, only: & - homogenization_Ngrains, & - mappingHomogenization, & - phaseAt, & - phase_source, & - phase_Nsources, & - SOURCE_vacancy_phenoplasticity_ID, & - SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID - use source_vacancy_phenoplasticity, only: & - source_vacancy_phenoplasticity_getRateAndItsTangent - use source_vacancy_irradiation, only: & - source_vacancy_irradiation_getRateAndItsTangent - use source_vacancy_thermalfluc, only: & - source_vacancy_thermalfluc_getRateAndItsTangent - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Cv - integer(pInt) :: & - phase, & - grain, & - source - real(pReal) :: & - CvDot, dCvDot_dCv, localCvDot, dLocalCvDot_dCv - - CvDot = 0.0_pReal - dCvDot_dCv = 0.0_pReal - do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el)) - phase = phaseAt(grain,ip,el) - do source = 1_pInt, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_vacancy_phenoplasticity_ID) - call source_vacancy_phenoplasticity_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) - - case (SOURCE_vacancy_irradiation_ID) - call source_vacancy_irradiation_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) - - case (SOURCE_vacancy_thermalfluc_ID) - call source_vacancy_thermalfluc_getRateAndItsTangent(localCvDot, dLocalCvDot_dCv, grain, ip, el) - - end select - CvDot = CvDot + localCvDot - dCvDot_dCv = dCvDot_dCv + dLocalCvDot_dCv - enddo - enddo - - CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal) - dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal) - -end subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized vacancy mobility tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function vacancyflux_cahnhilliard_getMobility33(ip,el) - use lattice, only: & - lattice_vacancyfluxMobility33 - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_push33ToRef - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - vacancyflux_cahnhilliard_getMobility33 - integer(pInt) :: & - grain - - vacancyflux_cahnhilliard_getMobility33 = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - vacancyflux_cahnhilliard_getMobility33 = vacancyflux_cahnhilliard_getMobility33 + & - crystallite_push33ToRef(grain,ip,el,lattice_vacancyfluxMobility33(:,:,material_phase(grain,ip,el))) - enddo - - vacancyflux_cahnhilliard_getMobility33 = & - vacancyflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function vacancyflux_cahnhilliard_getMobility33 - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized vacancy diffusion tensor in reference configuration -!-------------------------------------------------------------------------------------------------- -function vacancyflux_cahnhilliard_getDiffusion33(ip,el) - use lattice, only: & - lattice_vacancyfluxDiffusion33 - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - use crystallite, only: & - crystallite_push33ToRef - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - vacancyflux_cahnhilliard_getDiffusion33 - integer(pInt) :: & - grain - - vacancyflux_cahnhilliard_getDiffusion33 = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - vacancyflux_cahnhilliard_getDiffusion33 = vacancyflux_cahnhilliard_getDiffusion33 + & - crystallite_push33ToRef(grain,ip,el,lattice_vacancyfluxDiffusion33(:,:,material_phase(grain,ip,el))) - enddo - - vacancyflux_cahnhilliard_getDiffusion33 = & - vacancyflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function vacancyflux_cahnhilliard_getDiffusion33 - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized vacancy formation energy -!-------------------------------------------------------------------------------------------------- -real(pReal) function vacancyflux_cahnhilliard_getFormationEnergy(ip,el) - use lattice, only: & - lattice_vacancyFormationEnergy, & - lattice_vacancyVol, & - lattice_vacancySurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_phase - use mesh, only: & - mesh_element - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - integer(pInt) :: & - grain - - vacancyflux_cahnhilliard_getFormationEnergy = 0.0_pReal - do grain = 1, homogenization_Ngrains(mesh_element(3,el)) - vacancyflux_cahnhilliard_getFormationEnergy = vacancyflux_cahnhilliard_getFormationEnergy + & - lattice_vacancyFormationEnergy(material_phase(grain,ip,el))/ & - lattice_vacancyVol(material_phase(grain,ip,el))/ & - lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) - enddo - - vacancyflux_cahnhilliard_getFormationEnergy = & - vacancyflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal) - -end function vacancyflux_cahnhilliard_getFormationEnergy - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized vacancy entropy coefficient -!-------------------------------------------------------------------------------------------------- -real(pReal) function vacancyflux_cahnhilliard_getEntropicCoeff(ip,el) - use lattice, only: & - lattice_vacancyVol, & - lattice_vacancySurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_homog, & - material_phase, & - temperature, & - thermalMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - integer(pInt) :: & - grain - - vacancyflux_cahnhilliard_getEntropicCoeff = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homog(ip,el)) - vacancyflux_cahnhilliard_getEntropicCoeff = vacancyflux_cahnhilliard_getEntropicCoeff + & - kB/ & - lattice_vacancyVol(material_phase(grain,ip,el))/ & - lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) - enddo - - vacancyflux_cahnhilliard_getEntropicCoeff = & - vacancyflux_cahnhilliard_getEntropicCoeff* & - temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ & - real(homogenization_Ngrains(material_homog(ip,el)),pReal) - -end function vacancyflux_cahnhilliard_getEntropicCoeff - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized kinematic contribution to chemical potential -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCv, Cv, ip, el) - use lattice, only: & - lattice_vacancySurfaceEnergy - use material, only: & - homogenization_Ngrains, & - material_homog, & - phase_kinematics, & - phase_Nkinematics, & - material_phase, & - KINEMATICS_vacancy_strain_ID - use crystallite, only: & - crystallite_Tstar_v, & - crystallite_Fi0, & - crystallite_Fi - use kinematics_vacancy_strain, only: & - kinematics_vacancy_strain_ChemPotAndItsTangent - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Cv - real(pReal), intent(out) :: & - KPot, dKPot_dCv - real(pReal) :: & - my_KPot, my_dKPot_dCv - integer(pInt) :: & - grain, kinematics - - KPot = 0.0_pReal - dKPot_dCv = 0.0_pReal - do grain = 1_pInt,homogenization_Ngrains(material_homog(ip,el)) - do kinematics = 1_pInt, phase_Nkinematics(material_phase(grain,ip,el)) - select case (phase_kinematics(kinematics,material_phase(grain,ip,el))) - case (KINEMATICS_vacancy_strain_ID) - call kinematics_vacancy_strain_ChemPotAndItsTangent(my_KPot, my_dKPot_dCv, & - crystallite_Tstar_v(1:6,grain,ip,el), & - crystallite_Fi0(1:3,1:3,grain,ip,el), & - crystallite_Fi (1:3,1:3,grain,ip,el), & - grain,ip, el) - - case default - my_KPot = 0.0_pReal - my_dKPot_dCv = 0.0_pReal - - end select - KPot = KPot + my_KPot/lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) - dKPot_dCv = dKPot_dCv + my_dKPot_dCv/lattice_vacancySurfaceEnergy(material_phase(grain,ip,el)) - enddo - enddo - - KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal) - dKPot_dCv = dKPot_dCv/real(homogenization_Ngrains(material_homog(ip,el)),pReal) - -end subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief returns homogenized chemical potential and its tangent -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_cahnhilliard_getChemPotAndItsTangent(ChemPot,dChemPot_dCv,Cv,ip,el) - use numerics, only: & - vacancyBoundPenalty, & - vacancyPolyOrder - use material, only: & - mappingHomogenization, & - vacancyflux_typeInstance, & - porosity, & - porosityMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Cv - real(pReal), intent(out) :: & - ChemPot, & - dChemPot_dCv - real(pReal) :: & - VoidPhaseFrac, kBT, KPot, dKPot_dCv - integer(pInt) :: & - homog, o - - homog = mappingHomogenization(2,ip,el) - VoidPhaseFrac = porosity(homog)%p(porosityMapping(homog)%p(ip,el)) - kBT = vacancyflux_cahnhilliard_getEntropicCoeff(ip,el) - - ChemPot = vacancyflux_cahnhilliard_getFormationEnergy(ip,el) - dChemPot_dCv = 0.0_pReal - do o = 1_pInt, vacancyPolyOrder - ChemPot = ChemPot + kBT*((2.0_pReal*Cv - 1.0_pReal)**real(2_pInt*o-1_pInt,pReal))/ & - real(2_pInt*o-1_pInt,pReal) - dChemPot_dCv = dChemPot_dCv + 2.0_pReal*kBT*(2.0_pReal*Cv - 1.0_pReal)**real(2_pInt*o-2_pInt,pReal) - enddo - - ChemPot = VoidPhaseFrac*VoidPhaseFrac*ChemPot & - - 2.0_pReal*(1.0_pReal - Cv)*(1.0_pReal - VoidPhaseFrac)*(1.0_pReal - VoidPhaseFrac) - - dChemPot_dCv = VoidPhaseFrac*VoidPhaseFrac*dChemPot_dCv & - + 2.0_pReal*(1.0_pReal - VoidPhaseFrac)*(1.0_pReal - VoidPhaseFrac) - - call vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dCv, Cv, ip, el) - ChemPot = ChemPot + KPot - dChemPot_dCv = dChemPot_dCv + dKPot_dCv - - if (Cv < 0.0_pReal) then - ChemPot = ChemPot - 3.0_pReal*vacancyBoundPenalty*Cv*Cv - dChemPot_dCv = dChemPot_dCv - 6.0_pReal*vacancyBoundPenalty*Cv - elseif (Cv > 1.0_pReal) then - ChemPot = ChemPot + 3.0_pReal*vacancyBoundPenalty*(1.0_pReal - Cv)*(1.0_pReal - Cv) - dChemPot_dCv = dChemPot_dCv - 6.0_pReal*vacancyBoundPenalty*(1.0_pReal - Cv) - endif - - ChemPot = ChemPot* & - vacancyflux_cahnhilliard_thermalFluc(vacancyflux_typeInstance(homog))%p(mappingHomogenization(1,ip,el)) - dChemPot_dCv = dChemPot_dCv* & - vacancyflux_cahnhilliard_thermalFluc(vacancyflux_typeInstance(homog))%p(mappingHomogenization(1,ip,el)) - -end subroutine vacancyflux_cahnhilliard_getChemPotAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief updated vacancy concentration and its rate with solution from transport PDE -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_cahnhilliard_putVacancyConcAndItsRate(Cv,Cvdot,ip,el) - use material, only: & - mappingHomogenization, & - vacancyConc, & - vacancyConcRate, & - vacancyfluxMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Cv, & - Cvdot - integer(pInt) :: & - homog, & - offset - - homog = mappingHomogenization(2,ip,el) - offset = vacancyfluxMapping(homog)%p(ip,el) - vacancyConc (homog)%p(offset) = Cv - vacancyConcRate(homog)%p(offset) = Cvdot - -end subroutine vacancyflux_cahnhilliard_putVacancyConcAndItsRate - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of vacancy transport results -!-------------------------------------------------------------------------------------------------- -function vacancyflux_cahnhilliard_postResults(ip,el) - use material, only: & - mappingHomogenization, & - vacancyflux_typeInstance, & - vacancyConc, & - vacancyfluxMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(vacancyflux_cahnhilliard_sizePostResults(vacancyflux_typeInstance(mappingHomogenization(2,ip,el)))) :: & - vacancyflux_cahnhilliard_postResults - - integer(pInt) :: & - instance, homog, offset, o, c - - homog = mappingHomogenization(2,ip,el) - offset = vacancyfluxMapping(homog)%p(ip,el) - instance = vacancyflux_typeInstance(homog) - - c = 0_pInt - vacancyflux_cahnhilliard_postResults = 0.0_pReal - - do o = 1_pInt,vacancyflux_cahnhilliard_Noutput(instance) - select case(vacancyflux_cahnhilliard_outputID(o,instance)) - - case (vacancyConc_ID) - vacancyflux_cahnhilliard_postResults(c+1_pInt) = vacancyConc(homog)%p(offset) - c = c + 1 - end select - enddo -end function vacancyflux_cahnhilliard_postResults - -end module vacancyflux_cahnhilliard diff --git a/src/vacancyflux_isochempot.f90 b/src/vacancyflux_isochempot.f90 deleted file mode 100644 index 761a0ba22..000000000 --- a/src/vacancyflux_isochempot.f90 +++ /dev/null @@ -1,328 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for locally evolving vacancy concentration -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module vacancyflux_isochempot - use prec, only: & - pReal, & - pInt - - implicit none - private - integer(pInt), dimension(:), allocatable, public, protected :: & - vacancyflux_isochempot_sizePostResults !< cumulative size of post results - - integer(pInt), dimension(:,:), allocatable, target, public :: & - vacancyflux_isochempot_sizePostResult !< size of each post result output - - character(len=64), dimension(:,:), allocatable, target, public :: & - vacancyflux_isochempot_output !< name of each post result output - - integer(pInt), dimension(:), allocatable, target, public :: & - vacancyflux_isochempot_Noutput !< number of outputs per instance of this damage - - enum, bind(c) - enumerator :: undefined_ID, & - vacancyconc_ID - end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & - vacancyflux_isochempot_outputID !< ID of each post result output - - - public :: & - vacancyflux_isochempot_init, & - vacancyflux_isochempot_updateState, & - vacancyflux_isochempot_getSourceAndItsTangent, & - vacancyflux_isochempot_postResults - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_isochempot_init(fileUnit) -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use IO, only: & - IO_read, & - IO_lc, & - IO_getTag, & - IO_isBlank, & - IO_stringPos, & - IO_stringValue, & - IO_floatValue, & - IO_intValue, & - IO_warning, & - IO_error, & - IO_timeStamp, & - IO_EOF - use material, only: & - vacancyflux_type, & - vacancyflux_typeInstance, & - homogenization_Noutput, & - VACANCYFLUX_isochempot_label, & - VACANCYFLUX_isochempot_ID, & - material_homog, & - mappingHomogenization, & - vacancyfluxState, & - vacancyfluxMapping, & - vacancyConc, & - vacancyConcRate, & - vacancyflux_initialCv - use config, only: & - material_partHomogenization - - implicit none - integer(pInt), intent(in) :: fileUnit - - integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: maxNinstance,mySize=0_pInt,section,instance,o - integer(pInt) :: sizeState - integer(pInt) :: NofMyHomog - character(len=65536) :: & - tag = '', & - line = '' - - write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isochempot_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - maxNinstance = int(count(vacancyflux_type == VACANCYFLUX_isochempot_ID),pInt) - if (maxNinstance == 0_pInt) return - - allocate(vacancyflux_isochempot_sizePostResults(maxNinstance), source=0_pInt) - allocate(vacancyflux_isochempot_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0_pInt) - allocate(vacancyflux_isochempot_output (maxval(homogenization_Noutput),maxNinstance)) - vacancyflux_isochempot_output = '' - allocate(vacancyflux_isochempot_outputID (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID) - allocate(vacancyflux_isochempot_Noutput (maxNinstance), source=0_pInt) - - rewind(fileUnit) - section = 0_pInt - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization)! wind forward to - line = IO_read(fileUnit) - enddo - - parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of homog part - line = IO_read(fileUnit) - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') then ! stop at next part - line = IO_read(fileUnit, .true.) ! reset IO_read - exit - endif - if (IO_getTag(line,'[',']') /= '') then ! next homog section - section = section + 1_pInt ! advance homog section counter - cycle ! skip to next line - endif - - if (section > 0_pInt ) then; if (vacancyflux_type(section) == VACANCYFLUX_isochempot_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - - instance = vacancyflux_typeInstance(section) ! which instance of my vacancyflux is present homog - chunkPos = IO_stringPos(line) - tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key - select case(tag) - case ('(output)') - select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) - case ('vacancyconc') - vacancyflux_isochempot_Noutput(instance) = vacancyflux_isochempot_Noutput(instance) + 1_pInt - vacancyflux_isochempot_outputID(vacancyflux_isochempot_Noutput(instance),instance) = vacancyconc_ID - vacancyflux_isochempot_output(vacancyflux_isochempot_Noutput(instance),instance) = & - IO_lc(IO_stringValue(line,chunkPos,2_pInt)) - end select - - end select - endif; endif - enddo parsingFile - - initializeInstances: do section = 1_pInt, size(vacancyflux_type) - if (vacancyflux_type(section) == VACANCYFLUX_isochempot_ID) then - NofMyHomog=count(material_homog==section) - instance = vacancyflux_typeInstance(section) - -!-------------------------------------------------------------------------------------------------- -! Determine size of postResults array - outputsLoop: do o = 1_pInt,vacancyflux_isochempot_Noutput(instance) - select case(vacancyflux_isochempot_outputID(o,instance)) - case(vacancyconc_ID) - mySize = 1_pInt - end select - - if (mySize > 0_pInt) then ! any meaningful output found - vacancyflux_isochempot_sizePostResult(o,instance) = mySize - vacancyflux_isochempot_sizePostResults(instance) = vacancyflux_isochempot_sizePostResults(instance) + mySize - endif - enddo outputsLoop - -! allocate state arrays - sizeState = 1_pInt - vacancyfluxState(section)%sizeState = sizeState - vacancyfluxState(section)%sizePostResults = vacancyflux_isochempot_sizePostResults(instance) - allocate(vacancyfluxState(section)%state0 (sizeState,NofMyHomog), source=vacancyflux_initialCv(section)) - allocate(vacancyfluxState(section)%subState0(sizeState,NofMyHomog), source=vacancyflux_initialCv(section)) - allocate(vacancyfluxState(section)%state (sizeState,NofMyHomog), source=vacancyflux_initialCv(section)) - - nullify(vacancyfluxMapping(section)%p) - vacancyfluxMapping(section)%p => mappingHomogenization(1,:,:) - deallocate(vacancyConc(section)%p) - vacancyConc(section)%p => vacancyfluxState(section)%state(1,:) - deallocate(vacancyConcRate(section)%p) - allocate(vacancyConcRate(section)%p(NofMyHomog), source=0.0_pReal) - - endif - - enddo initializeInstances -end subroutine vacancyflux_isochempot_init - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates change in vacancy concentration based on local vacancy generation model -!-------------------------------------------------------------------------------------------------- -function vacancyflux_isochempot_updateState(subdt, ip, el) - use numerics, only: & - err_vacancyflux_tolAbs, & - err_vacancyflux_tolRel - use material, only: & - mappingHomogenization, & - vacancyflux_typeInstance, & - vacancyfluxState, & - vacancyConc, & - vacancyConcRate, & - vacancyfluxMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - subdt - logical, dimension(2) :: & - vacancyflux_isochempot_updateState - integer(pInt) :: & - homog, & - offset, & - instance - real(pReal) :: & - Cv, Cvdot, dCvDot_dCv - - homog = mappingHomogenization(2,ip,el) - offset = mappingHomogenization(1,ip,el) - instance = vacancyflux_typeInstance(homog) - - Cv = vacancyfluxState(homog)%subState0(1,offset) - call vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el) - Cv = Cv + subdt*Cvdot - - vacancyflux_isochempot_updateState = [ abs(Cv - vacancyfluxState(homog)%state(1,offset)) & - <= err_vacancyflux_tolAbs & - .or. abs(Cv - vacancyfluxState(homog)%state(1,offset)) & - <= err_vacancyflux_tolRel*abs(vacancyfluxState(homog)%state(1,offset)), & - .true.] - - vacancyConc (homog)%p(vacancyfluxMapping(homog)%p(ip,el)) = Cv - vacancyConcRate(homog)%p(vacancyfluxMapping(homog)%p(ip,el)) = & - (vacancyfluxState(homog)%state(1,offset) - vacancyfluxState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal)) - -end function vacancyflux_isochempot_updateState - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized vacancy driving forces -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv, ip, el) - use material, only: & - homogenization_Ngrains, & - mappingHomogenization, & - phaseAt, & - phase_source, & - phase_Nsources, & - SOURCE_vacancy_phenoplasticity_ID, & - SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID - use source_vacancy_phenoplasticity, only: & - source_vacancy_phenoplasticity_getRateAndItsTangent - use source_vacancy_irradiation, only: & - source_vacancy_irradiation_getRateAndItsTangent - use source_vacancy_thermalfluc, only: & - source_vacancy_thermalfluc_getRateAndItsTangent - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - Cv - integer(pInt) :: & - phase, & - grain, & - source - real(pReal) :: & - CvDot, dCvDot_dCv, localCvDot, dLocalCvDot_dCv - - CvDot = 0.0_pReal - dCvDot_dCv = 0.0_pReal - do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el)) - phase = phaseAt(grain,ip,el) - do source = 1_pInt, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_vacancy_phenoplasticity_ID) - call source_vacancy_phenoplasticity_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) - - case (SOURCE_vacancy_irradiation_ID) - call source_vacancy_irradiation_getRateAndItsTangent (localCvDot, dLocalCvDot_dCv, grain, ip, el) - - case (SOURCE_vacancy_thermalfluc_ID) - call source_vacancy_thermalfluc_getRateAndItsTangent(localCvDot, dLocalCvDot_dCv, grain, ip, el) - - end select - CvDot = CvDot + localCvDot - dCvDot_dCv = dCvDot_dCv + dLocalCvDot_dCv - enddo - enddo - - CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal) - dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal) - -end subroutine vacancyflux_isochempot_getSourceAndItsTangent - -!-------------------------------------------------------------------------------------------------- -!> @brief return array of vacancy transport results -!-------------------------------------------------------------------------------------------------- -function vacancyflux_isochempot_postResults(ip,el) - use material, only: & - mappingHomogenization, & - vacancyflux_typeInstance, & - vacancyConc, & - vacancyfluxMapping - - implicit none - integer(pInt), intent(in) :: & - ip, & !< integration point - el !< element - real(pReal), dimension(vacancyflux_isochempot_sizePostResults(vacancyflux_typeInstance(mappingHomogenization(2,ip,el)))) :: & - vacancyflux_isochempot_postResults - - integer(pInt) :: & - instance, homog, offset, o, c - - homog = mappingHomogenization(2,ip,el) - offset = vacancyfluxMapping(homog)%p(ip,el) - instance = vacancyflux_typeInstance(homog) - - c = 0_pInt - vacancyflux_isochempot_postResults = 0.0_pReal - - do o = 1_pInt,vacancyflux_isochempot_Noutput(instance) - select case(vacancyflux_isochempot_outputID(o,instance)) - - case (vacancyconc_ID) - vacancyflux_isochempot_postResults(c+1_pInt) = vacancyConc(homog)%p(offset) - c = c + 1 - end select - enddo -end function vacancyflux_isochempot_postResults - -end module vacancyflux_isochempot diff --git a/src/vacancyflux_isoconc.f90 b/src/vacancyflux_isoconc.f90 deleted file mode 100644 index 135509aa1..000000000 --- a/src/vacancyflux_isoconc.f90 +++ /dev/null @@ -1,62 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for constant vacancy concentration -!-------------------------------------------------------------------------------------------------- -module vacancyflux_isoconc - - implicit none - private - - public :: & - vacancyflux_isoconc_init - -contains - -!-------------------------------------------------------------------------------------------------- -!> @brief allocates all neccessary fields, reads information from material configuration file -!-------------------------------------------------------------------------------------------------- -subroutine vacancyflux_isoconc_init() -#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 - use, intrinsic :: iso_fortran_env, only: & - compiler_version, & - compiler_options -#endif - use prec, only: & - pReal, & - pInt - use IO, only: & - IO_timeStamp - use material - use config - - implicit none - integer(pInt) :: & - homog, & - NofMyHomog - - write(6,'(/,a)') ' <<<+- vacancyflux_'//VACANCYFLUX_isoconc_label//' init -+>>>' - write(6,'(a15,a)') ' Current time: ',IO_timeStamp() -#include "compilation_info.f90" - - initializeInstances: do homog = 1_pInt, material_Nhomogenization - - myhomog: if (vacancyflux_type(homog) == VACANCYFLUX_isoconc_ID) then - NofMyHomog = count(material_homog == homog) - vacancyfluxState(homog)%sizeState = 0_pInt - vacancyfluxState(homog)%sizePostResults = 0_pInt - allocate(vacancyfluxState(homog)%state0 (0_pInt,NofMyHomog)) - allocate(vacancyfluxState(homog)%subState0(0_pInt,NofMyHomog)) - allocate(vacancyfluxState(homog)%state (0_pInt,NofMyHomog)) - - deallocate(vacancyConc (homog)%p) - allocate (vacancyConc (homog)%p(1), source=vacancyflux_initialCv(homog)) - deallocate(vacancyConcRate(homog)%p) - allocate (vacancyConcRate(homog)%p(1), source=0.0_pReal) - - endif myhomog - enddo initializeInstances - - -end subroutine vacancyflux_isoconc_init - -end module vacancyflux_isoconc