From 06f8697bd15f577c709725c06929d8a992482d11 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 7 Dec 2022 21:43:05 +0100 Subject: [PATCH 1/6] manage tabulated data with linear interpolation in-between points --- src/IO.f90 | 2 + src/Marc/DAMASK_Marc.f90 | 1 + src/Marc/materialpoint_Marc.f90 | 2 + src/materialpoint.f90 | 2 + src/tables.f90 | 141 ++++++++++++++++++++++++++++++++ 5 files changed, 148 insertions(+) create mode 100644 src/tables.f90 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/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/tables.f90 b/src/tables.f90 new file mode 100644 index 000000000..cdd24628b --- /dev/null +++ b/src/tables.f90 @@ -0,0 +1,141 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @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 (min(size(x),size(y))< 1) call IO_error(603,ext_msg='no data specified') + if (size(x) /= size(y)) call IO_error(603,ext_msg='non-matching shape of tabulated data') + if (size(x) /=1) then + if (any(x(1:size(x)-1) -x(2:size(x)) > 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 Evaluate a table. +!-------------------------------------------------------------------------------------------------- +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%x(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.,2.,2.,1.],pReal), & + x_eval = real([0.,.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.],pReal), & + y_true = real([0.,.5,1.,1.5,2.,2. ,2.,1.5,1.,.5, 0.],pReal) + integer :: i + type(tDict), pointer :: dict + type(tList), pointer :: l_x, l_y + real :: r + + + call random_number(r) + 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))) error stop 'table eval/values' + end do + + + l_x => YAML_parse_str_asList('[1, 2, 3, 4]') + l_y => YAML_parse_str_asList('[1, 2, 2, 1]') + 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 From 7b986f6f7744ee2b621c3bc2514201dfa938037a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 7 Dec 2022 21:55:20 +0100 Subject: [PATCH 2/6] using new table functionality --- src/phase.f90 | 1 + src/phase_thermal_externalheat.f90 | 27 +++++---------------------- 2 files changed, 6 insertions(+), 22 deletions(-) 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..b4b6976c8 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_n','f_T') 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 From cca5d6d495507c794fcdbe3b43462d74fda66b2c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 8 Dec 2022 16:57:26 +0100 Subject: [PATCH 3/6] impoved naming --- PRIVATE | 2 +- .../phase/thermal/source/externalheat_ramp-and-hold.yaml | 4 ++-- src/phase_thermal_externalheat.f90 | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) 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/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index b4b6976c8..5970a5894 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -60,7 +60,7 @@ module function externalheat_init(source_length) result(mySources) associate(prm => param(ph)) src => sources%get_dict(so) - prm%f = table(src,'t_n','f_T') + prm%f = table(src,'t','f') Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) From a3128d402186d218691e89e502c2b551e68a3209 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 8 Dec 2022 20:40:20 +0100 Subject: [PATCH 4/6] avoid failing tests --- src/tables.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/tables.f90 b/src/tables.f90 index cdd24628b..867433ab2 100644 --- a/src/tables.f90 +++ b/src/tables.f90 @@ -115,19 +115,19 @@ subroutine selfTest() integer :: i type(tDict), pointer :: dict type(tList), pointer :: l_x, l_y - real :: r + real(pReal) :: r call random_number(r) 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))) error stop 'table eval/values' + 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]') - l_y => YAML_parse_str_asList('[1, 2, 2, 1]') + l_x => YAML_parse_str_asList('[1, 2, 3, 4]'//IO_EOL) + l_y => YAML_parse_str_asList('[1, 2, 2, 1]'//IO_EOL) allocate(dict) call dict%set('t',l_x) call dict%set('T',l_y) From a9cb81b1cb8bdc26b1e592b4f2d70d5fb676d374 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 8 Dec 2022 20:40:47 +0100 Subject: [PATCH 5/6] functions are meant for proper files (ending with EOL) --- src/YAML_parse.f90 | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) 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 From 4830ea19c947ab34e2e80c6b99b36f4485a71470 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 9 Dec 2022 10:00:00 -0500 Subject: [PATCH 6/6] bug fix;more elaborate self test --- src/tables.f90 | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/tables.f90 b/src/tables.f90 index 867433ab2..c62082705 100644 --- a/src/tables.f90 +++ b/src/tables.f90 @@ -1,5 +1,6 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven +!> @author Philip Eisenlohr, Michigan State University !> @brief Tabular representation of variable data. !-------------------------------------------------------------------------------------------------- module tables @@ -50,11 +51,12 @@ function table_from_values(x,y) result(t) type(tTable) :: t - if (min(size(x),size(y))< 1) call IO_error(603,ext_msg='no data specified') - if (size(x) /= size(y)) call IO_error(603,ext_msg='non-matching shape of tabulated data') - if (size(x) /=1) then - if (any(x(1:size(x)-1) -x(2:size(x)) > 0.0_pReal)) & - call IO_error(603,ext_msg='ordinate data does not increase monotonically') + 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 @@ -79,7 +81,7 @@ end function table_from_dict !-------------------------------------------------------------------------------------------------- -!> @brief Evaluate a table. +!> @brief Linearly interpolate/extrapolate tabular data. !-------------------------------------------------------------------------------------------------- pure function eval(self,x) result(y) @@ -91,11 +93,11 @@ pure function eval(self,x) result(y) if (size(self%x) == 1) then - y = self%x(1) + y = self%y(1) else i = max(1,min(findloc(self%x YAML_parse_str_asList('[1, 2, 3, 4]'//IO_EOL) - l_y => YAML_parse_str_asList('[1, 2, 2, 1]'//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)