diff --git a/PRIVATE b/PRIVATE index 599d5accf..e9254133c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 599d5accf4f5249257972cef997e8078e857c5a1 +Subproject commit e9254133c1e9ea3855a4fd17078d4c6f7d8728b1 diff --git a/examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml b/examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml index 3567f55cf..e6cc8f509 100644 --- a/examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml +++ b/examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml @@ -1,4 +1,4 @@ type: externalheat -f_T: [1, 1, 0, 0] -t_n: [0, 500, 500.001, 1000] +f: [1, 1, 0, 0] +t: [0, 500, 500.001, 1000] diff --git a/src/IO.f90 b/src/IO.f90 index 1643cac71..8cc350b10 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -484,6 +484,8 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2) ! user errors case (602) msg = 'invalid selection for debug' + case (603) + msg = 'invalid data for table' !------------------------------------------------------------------------------------------------ ! errors related to YAML data diff --git a/src/Marc/DAMASK_Marc.f90 b/src/Marc/DAMASK_Marc.f90 index 33d1268ef..024fd410b 100644 --- a/src/Marc/DAMASK_Marc.f90 +++ b/src/Marc/DAMASK_Marc.f90 @@ -153,6 +153,7 @@ end module DAMASK_interface #include "../math.f90" #include "../rotations.f90" #include "../polynomials.f90" +#include "../tables.f90" #include "../lattice.f90" #include "element.f90" #include "../geometry_plastic_nonlocal.f90" diff --git a/src/Marc/materialpoint_Marc.f90 b/src/Marc/materialpoint_Marc.f90 index dc878f507..2b910217d 100644 --- a/src/Marc/materialpoint_Marc.f90 +++ b/src/Marc/materialpoint_Marc.f90 @@ -15,6 +15,7 @@ module materialpoint_Marc use math use rotations use polynomials + use tables use lattice use material use phase @@ -85,6 +86,7 @@ subroutine materialpoint_initAll() call math_init() call rotations_init() call polynomials_init() + call tables_init() call lattice_init() call discretization_Marc_init() call material_init(.false.) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 04ff4d5cc..b1f5aaf71 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -55,6 +55,7 @@ end subroutine YAML_parse_init !-------------------------------------------------------------------------------------------------- !> @brief Parse a YAML string with list as root into a a structure of nodes. +!> @details The string needs to end with a newline (unless using libfyaml). !-------------------------------------------------------------------------------------------------- function YAML_parse_str_asList(str) result(list) @@ -72,6 +73,7 @@ end function YAML_parse_str_asList !-------------------------------------------------------------------------------------------------- !> @brief Parse a YAML string with dict as root into a a structure of nodes. +!> @details The string needs to end with a newline (unless using libfyaml). !-------------------------------------------------------------------------------------------------- function YAML_parse_str_asDict(str) result(dict) @@ -815,7 +817,8 @@ end subroutine decide !-------------------------------------------------------------------------------------------------- -! @brief Convert all block style YAML parts to flow style. +!> @brief Convert all block style YAML parts to flow style. +!> @details The input needs to end with a newline. !-------------------------------------------------------------------------------------------------- function to_flow(blck) @@ -998,6 +1001,21 @@ subroutine selfTest() if (.not. to_flow(block_flow) == mixed_flow) error stop 'to_flow' end block basic_mixed + parse: block + + type(tDict), pointer :: dict + type(tList), pointer :: list + character(len=*), parameter :: & + lst = '[1, 2, 3, 4]', & + dct = '{a: 1, b: 2}' + + list => YAML_parse_str_asList(lst//IO_EOL) + if (list%asFormattedString() /= lst) error stop 'str_asList' + dict => YAML_parse_str_asDict(dct//IO_EOL) + if (dict%asFormattedString() /= dct) error stop 'str_asDict' + + end block parse + end subroutine selfTest #endif diff --git a/src/materialpoint.f90 b/src/materialpoint.f90 index a0f266afa..8ce0e15a1 100644 --- a/src/materialpoint.f90 +++ b/src/materialpoint.f90 @@ -18,6 +18,7 @@ module materialpoint use math use rotations use polynomials + use tables use lattice use material use phase @@ -60,6 +61,7 @@ subroutine materialpoint_initAll() call math_init() call rotations_init() call polynomials_init() + call tables_init() call lattice_init() #if defined(MESH) call discretization_mesh_init(restart=CLI_restartInc>0) diff --git a/src/phase.f90 b/src/phase.f90 index 12e5ea924..f7088b892 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -9,6 +9,7 @@ module phase use math use rotations use polynomials + use tables use IO use config use material diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index 9a58595ae..5970a5894 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -11,11 +11,7 @@ submodule(phase:thermal) externalheat source_thermal_externalheat_offset !< which source is my current thermal dissipation mechanism? type :: tParameters !< container type for internal constitutive parameters - real(pReal), dimension(:), allocatable :: & - t_n, & - f_T - integer :: & - nIntervals + type(tTable) :: f end type tParameters type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances) @@ -64,10 +60,7 @@ module function externalheat_init(source_length) result(mySources) associate(prm => param(ph)) src => sources%get_dict(so) - prm%t_n = src%get_as1dFloat('t_n') - prm%nIntervals = size(prm%t_n) - 1 - - prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n)) + prm%f = table(src,'t','f') Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) @@ -111,23 +104,13 @@ module function externalheat_f_T(ph,en) result(f_T) f_T integer :: & - so, interval - real(pReal) :: & - frac_time + so + so = source_thermal_externalheat_offset(ph) associate(prm => param(ph)) - do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (thermalState(ph)%p(so)%state(1,en) - prm%t_n(interval)) & - / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment - if ( (frac_time < 0.0_pReal .and. interval == 1) & - .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & - .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) & - f_T = prm%f_T(interval ) * (1.0_pReal - frac_time) + & - prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... - ! ...or extrapolate if outside of bounds - end do + f_T = prm%f%at(thermalState(ph)%p(so)%state(1,en)) end associate end function externalheat_f_T diff --git a/src/tables.f90 b/src/tables.f90 new file mode 100644 index 000000000..c62082705 --- /dev/null +++ b/src/tables.f90 @@ -0,0 +1,145 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @author Philip Eisenlohr, Michigan State University +!> @brief Tabular representation of variable data. +!-------------------------------------------------------------------------------------------------- +module tables + use prec + use IO + use YAML_parse + use YAML_types + + implicit none(type,external) + private + + type, public :: tTable + real(pReal), dimension(:), allocatable :: x,y + contains + procedure, public :: at => eval + end type tTable + + interface table + module procedure table_from_values + module procedure table_from_dict + end interface table + + public :: & + table, & + tables_init + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief Run self-test. +!-------------------------------------------------------------------------------------------------- +subroutine tables_init() + + print'(/,1x,a)', '<<<+- tables init -+>>>'; flush(IO_STDOUT) + + call selfTest() + +end subroutine tables_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Initialize a table from values. +!-------------------------------------------------------------------------------------------------- +function table_from_values(x,y) result(t) + + real(pReal), dimension(:), intent(in) :: x,y + type(tTable) :: t + + + if (size(x) < 1) call IO_error(603,ext_msg='missing tabulated x data') + if (size(y) < 1) call IO_error(603,ext_msg='missing tabulated y data') + if (size(x) /= size(y)) call IO_error(603,ext_msg='shape mismatch in tabulated data') + if (size(x) /= 1) then + if (any(x(2:size(x))-x(1:size(x)-1) <= 0.0_pReal)) & + call IO_error(603,ext_msg='ordinate data does not increase monotonically') + end if + + t%x = x + t%y = y + +end function table_from_values + + +!-------------------------------------------------------------------------------------------------- +!> @brief Initialize a table from a dictionary with values. +!-------------------------------------------------------------------------------------------------- +function table_from_dict(dict,x_label,y_label) result(t) + + type(tDict), intent(in) :: dict + character(len=*), intent(in) :: x_label, y_label + type(tTable) :: t + + + t = tTable(dict%get_as1dFloat(x_label),dict%get_as1dFloat(y_label)) + +end function table_from_dict + + +!-------------------------------------------------------------------------------------------------- +!> @brief Linearly interpolate/extrapolate tabular data. +!-------------------------------------------------------------------------------------------------- +pure function eval(self,x) result(y) + + class(tTable), intent(in) :: self + real(pReal), intent(in) :: x + real(pReal) :: y + + integer :: i + + + if (size(self%x) == 1) then + y = self%y(1) + else + i = max(1,min(findloc(self%x @brief Check correctness of table functionality. +!-------------------------------------------------------------------------------------------------- +subroutine selfTest() + + type(tTable) :: t + real(pReal), dimension(*), parameter :: & + x = real([ 1., 2., 3., 4.],pReal), & + y = real([ 1., 3., 2.,-2.],pReal), & + x_eval = real([ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0],pReal), & + y_true = real([-1.0, 0.0, 1.0, 2.0, 3.0, 2.5 ,2.0, 0.0,-2.0,-4.0,-6.0],pReal) + integer :: i + type(tDict), pointer :: dict + type(tList), pointer :: l_x, l_y + real(pReal) :: r + + + call random_number(r) + t = table(real([0.],pReal),real([r],pReal)) + if (dNeq(r,t%at(r),1.0e-9_pReal)) error stop 'table eval/mono' + + r = r-0.5_pReal + t = table(x+r,y) + do i = 1, size(x_eval) + if (dNeq(y_true(i),t%at(x_eval(i)+r),1.0e-9_pReal)) error stop 'table eval/values' + end do + + l_x => YAML_parse_str_asList('[1, 2, 3, 4]'//IO_EOL) + l_y => YAML_parse_str_asList('[1, 3, 2,-2]'//IO_EOL) + allocate(dict) + call dict%set('t',l_x) + call dict%set('T',l_y) + t = table(dict,'t','T') + do i = 1, size(x_eval) + if (dNeq(y_true(i),t%at(x_eval(i)))) error stop 'table eval/dict' + end do + +end subroutine selfTest + +end module tables