Merge branch '201-external-heat-syntax' into 'development'

Resolve "external heat: Syntax"

Closes #201 and #107

See merge request damask/DAMASK!683

PRIVATE

@@ -1 +1 @@
-Subproject commit 599d5accf4f5249257972cef997e8078e857c5a1
+Subproject commit e9254133c1e9ea3855a4fd17078d4c6f7d8728b1

examples/config/phase/thermal/source/externalheat_ramp-and-hold.yaml

@@ -1,4 +1,4 @@
 type: externalheat
 
-f_T: [1, 1, 0, 0]
+f: [1, 1, 0, 0]
-t_n: [0, 500, 500.001, 1000]
+t: [0, 500, 500.001, 1000]

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

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"

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.)

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
 

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)

src/phase.f90

@@ -9,6 +9,7 @@ module phase
   use math
   use rotations
   use polynomials
+  use tables
   use IO
   use config
   use material

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 :: &
+    type(tTable) :: f
-      t_n, &
-      f_T
-    integer :: &
-     nIntervals
   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%f = table(src,'t','f')
-          prm%nIntervals = size(prm%t_n) - 1
-
-          prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n))
 
           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
+    so
-  real(pReal) :: &
+
-    frac_time
 
   so = source_thermal_externalheat_offset(ph)
 
   associate(prm => param(ph))
-    do interval = 1, prm%nIntervals                                                                 ! scan through all rate segments
+    f_T = prm%f%at(thermalState(ph)%p(so)%state(1,en))
-      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
   end associate
 
 end function externalheat_f_T

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<x,.true.,dim=1,back=.true.),size(self%x)-1))
+    y = self%y(i) &
+      + (x-self%x(i)) * (self%y(i+1)-self%y(i)) / (self%x(i+1)-self%x(i))
+  end if
+
+end function eval
+
+
+!--------------------------------------------------------------------------------------------------
+!> @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