Merge remote-tracking branch 'origin/development' into restructure-numerics

.gitlab-ci.yml

@@ -218,8 +218,10 @@ grid_performance:
           git push
         fi
 
+
+###################################################################################################
 update_plots:
-  stage: statistics
+  stage: finalize
   script:
     - cd $(mktemp -d)
     - git clone -q git@git.damask.mpie.de:damask/statistics.git .
@@ -234,8 +236,6 @@ update_plots:
   only:
     - development
 
-
-###################################################################################################
 update_revision:
   stage: finalize
   before_script:

VERSION

@@ -1 +1 @@
-3.0.0-alpha7-668-gbfc2a76a0
+3.0.0-alpha7-671-g0cace1ef8

src/IO.f90

@@ -548,6 +548,8 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
 
 !--------------------------------------------------------------------------------------------------
 ! user errors
+    case (600)
+      msg = 'only one source entry allowed'
     case (603)
       msg = 'invalid data for table'
     case (610)

src/Marc/DAMASK_Marc.f90

@@ -153,11 +153,10 @@ end module DAMASK_interface
 #include "../phase_mechanical_plastic_dislotungsten.f90"
 #include "../phase_mechanical_plastic_nonlocal.f90"
 #include "../phase_mechanical_eigen.f90"
-#include "../phase_mechanical_eigen_cleavageopening.f90"
 #include "../phase_mechanical_eigen_thermalexpansion.f90"
 #include "../phase_thermal.f90"
-#include "../phase_thermal_dissipation.f90"
+#include "../phase_thermal_source_dissipation.f90"
-#include "../phase_thermal_externalheat.f90"
+#include "../phase_thermal_source_externalheat.f90"
 #include "../phase_damage.f90"
 #include "../phase_damage_isobrittle.f90"
 #include "../phase_damage_anisobrittle.f90"

src/grid/DAMASK_grid.f90

@@ -75,7 +75,7 @@ program DAMASK_grid
     cutBack = .false.,&
     sig
   integer :: &
-    i, j, m, field, &
+    i, j, field, &
     errorID = 0, &
     cutBackLevel = 0, &                                                                             !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
     stepFraction = 0, &                                                                             !< fraction of current time interval
@@ -110,15 +110,8 @@ program DAMASK_grid
     load, &
     num_solver, &
     num_grid, &
-    load_step, &
+    solver
-    solver, &
-    step_bc, &
-    step_mech, &
-    step_discretization
   type(tList), pointer :: &
-#ifdef __INTEL_LLVM_COMPILER
-    tensor, &
-#endif
     load_steps
   character(len=:), allocatable :: &
     fileContent, fname
@@ -210,13 +203,251 @@ program DAMASK_grid
     ID(field) = FIELD_DAMAGE_ID
   end if damageActive
 
+!--------------------------------------------------------------------------------------------------
+! doing initialization depending on active solvers
+  call spectral_utilities_init()
+
+  do field = 2, nActiveFields
+    select case (ID(field))
+
+      case (FIELD_THERMAL_ID)
+        call grid_thermal_spectral_init(num_grid)
+
+      case (FIELD_DAMAGE_ID)
+        call grid_damage_spectral_init(num_grid)
+
+    end select
+  end do
+
+  call mechanical_init(num_grid)
+  call config_numerics_deallocate()
 
 !--------------------------------------------------------------------------------------------------
-  load_steps => load%get_list('loadstep')
+! write header of output file
-  allocate(loadCases(load_steps%length))                                                            ! array of load cases
+  if (worldrank == 0) then
+    writeHeader: if (CLI_restartInc < 1) then
+      open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
+      write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded StagIterationsNeeded' ! statistics file
+    else writeHeader
+      open(newunit=statUnit,file=trim(getSolverJobName())//&
+                                  '.sta',form='FORMATTED', position='APPEND', status='OLD')
+    end if writeHeader
+  end if
 
+  writeUndeformed: if (CLI_restartInc < 1) then
+    print'(/,1x,a)', '... saving initial configuration ..........................................'
+    flush(IO_STDOUT)
+    call materialpoint_result(0,0.0_pREAL)
+  end if writeUndeformed
+
+  loadCases = parseLoadSteps(load%get_list('loadstep'))
+
+  loadCaseLooping: do l = 1, size(loadCases)
+    t_0 = t                                                                                         ! load case start time
+    guess = loadCases(l)%estimate_rate                                                              ! change of load case? homogeneous guess for the first inc
+
+    incLooping: do inc = 1, loadCases(l)%N
+      totalIncsCounter = totalIncsCounter + 1
+
+!--------------------------------------------------------------------------------------------------
+! forwarding time
+      Delta_t_prev = Delta_t                                                                        ! last time intervall that brought former inc to an end
+      if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then                                           ! linear scale
+        Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL)
+      else
+        Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) &
+                                 / (1.0_pREAL-loadCases(l)%r**loadCases(l)%N)
+      end if
+      Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL)                      ! depending on cut back level, decrease time step
+
+      skipping: if (totalIncsCounter <= CLI_restartInc) then                                        ! not yet at restart inc?
+        t = t + Delta_t                                                                             ! just advance time, skip already performed calculation
+        guess = .true.                                                                              ! QUESTION:why forced guessing instead of inheriting loadcase preference
+      else skipping
+        stepFraction = 0                                                                            ! fraction scaled by stepFactor**cutLevel
+
+        subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
+          t_remaining = loadCases(l)%t + t_0 - t
+          t = t + Delta_t                                                                           ! forward target time
+          stepFraction = stepFraction + 1                                                           ! count step
+
+!--------------------------------------------------------------------------------------------------
+! report beginning of new step
+          print'(/,1x,a)', '###########################################################################'
+          print'(1x,a,1x,es12.5,6(a,i0))', &
+                  'Time', t, &
+                  's: Increment ', inc,'/',loadCases(l)%N,&
+                  '-', stepFraction,'/',subStepFactor**cutBackLevel,&
+                  ' of load case ', l,'/',size(loadCases)
+          write(incInfo,'(4(a,i0))') &
+                  'Increment ',totalIncsCounter,'/',sum(loadCases%N),&
+                  '-', stepFraction,'/',subStepFactor**cutBackLevel
+          flush(IO_STDOUT)
+
+!--------------------------------------------------------------------------------------------------
+! forward fields
+          do field = 1, nActiveFields
+            select case(ID(field))
+              case(FIELD_MECH_ID)
+                call mechanical_forward (&
+                        cutBack,guess,Delta_t,Delta_t_prev,t_remaining, &
+                        deformation_BC = loadCases(l)%deformation, &
+                        stress_BC      = loadCases(l)%stress, &
+                        rotation_BC    = loadCases(l)%rot)
+
+              case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
+              case(FIELD_DAMAGE_ID);  call grid_damage_spectral_forward(cutBack)
+            end select
+          end do
+          if (.not. cutBack) call materialpoint_forward
+
+!--------------------------------------------------------------------------------------------------
+! solve fields
+          stagIter = 1
+          stagIterate = .true.
+          do while (stagIterate)
+
+            if (nActiveFields > 1) print'(/,1x,a,i0)', 'Staggered Iteration ',stagIter
+            do field = 1, nActiveFields
+              select case(ID(field))
+                case(FIELD_MECH_ID)
+                  solres(field) = mechanical_solution(incInfo)
+                case(FIELD_THERMAL_ID)
+                  solres(field) = grid_thermal_spectral_solution(Delta_t)
+                case(FIELD_DAMAGE_ID)
+                  solres(field) = grid_damage_spectral_solution(Delta_t)
+              end select
+
+              if (.not. solres(field)%converged) exit                                               ! no solution found
+
+            end do
+            stagIter = stagIter + 1
+            stagIterate =            stagIter <= stagItMax &
+                         .and.       all(solres(:)%converged) &
+                         .and. .not. all(solres(:)%stagConverged)                                   ! stationary with respect to staggered iteration
+          end do
+
+!--------------------------------------------------------------------------------------------------
+! check solution and either advance or retry with smaller timestep
+
+          if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) &                           ! converged
+               .and. .not. solres(1)%termIll) then                                                  ! and acceptable solution found
+            call mechanical_updateCoords()
+            Delta_t_prev = Delta_t
+            cutBack = .false.
+            guess = .true.                                                                          ! start guessing after first converged (sub)inc
+            if (worldrank == 0) then
+              write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
+                                solres(1)%converged, solres(1)%iterationsNeeded, StagIter
+              flush(statUnit)
+            end if
+          elseif (cutBackLevel < maxCutBack) then                                                   ! further cutbacking tolerated?
+            cutBack = .true.
+            stepFraction = (stepFraction - 1) * subStepFactor                                       ! adjust to new denominator
+            cutBackLevel = cutBackLevel + 1
+            t = t - Delta_t
+            Delta_t = Delta_t/real(subStepFactor,pREAL)                                             ! cut timestep
+            print'(/,1x,a)', 'cutting back '
+          else                                                                                      ! no more options to continue
+            if (worldrank == 0) close(statUnit)
+            call IO_error(950)
+          end if
+
+        end do subStepLooping
+
+        cutBackLevel = max(0, cutBackLevel - 1)                                                     ! try half number of subincs next inc
+
+        if (all(solres(:)%converged)) then
+          print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'converged'
+        else
+          print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged'
+        end if; flush(IO_STDOUT)
+
+        call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
+        if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+        if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then
+          print'(/,1x,a)', '... saving results ........................................................'
+          flush(IO_STDOUT)
+          call materialpoint_result(totalIncsCounter,t)
+        end if
+        if (sig) call signal_setSIGUSR1(.false.)
+        call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
+        if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+        if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then
+          do field = 1, nActiveFields
+            select case (ID(field))
+              case(FIELD_MECH_ID)
+                call mechanical_restartWrite()
+              case(FIELD_THERMAL_ID)
+                call grid_thermal_spectral_restartWrite()
+              case(FIELD_DAMAGE_ID)
+                call grid_damage_spectral_restartWrite()
+            end select
+          end do
+          call materialpoint_restartWrite()
+        end if
+        if (sig) call signal_setSIGUSR2(.false.)
+        call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
+        if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+        if (sig) exit loadCaseLooping
+      end if skipping
+
+    end do incLooping
+
+  end do loadCaseLooping
+
+
+!--------------------------------------------------------------------------------------------------
+! report summary of whole calculation
+  print'(/,1x,a)', '###########################################################################'
+  if (worldrank == 0) close(statUnit)
+
+  call quit(0)                                                                                      ! no complains ;)
+
+
+contains
+
+subroutine getMaskedTensor(values,mask,tensor)
+
+  real(pREAL), intent(out), dimension(3,3) :: values
+  logical,     intent(out), dimension(3,3) :: mask
+  type(tList), pointer :: tensor
+
+  type(tList), pointer :: row
+  integer :: i,j
+
+
+  values = 0.0_pREAL
+  do i = 1,3
+    row => tensor%get_list(i)
+    do j = 1,3
+      mask(i,j) = row%get_asStr(j) == 'x'
+      if (.not. mask(i,j)) values(i,j) = row%get_asReal(j)
+    end do
+  end do
+
+end subroutine getMaskedTensor
+
+
+function parseLoadsteps(load_steps) result(loadCases)
+
+  type(tList), intent(in), target :: load_steps
+  type(tLoadCase), allocatable, dimension(:) :: loadCases                                           !< array of all load cases
+
+  integer :: l,m
+  type(tDict), pointer :: &
+    load_step, &
+    step_bc, &
+    step_mech, &
+    step_discretization
+#ifdef __INTEL_LLVM_COMPILER
+  type(tList), pointer :: &
+    tensor
+#endif
+
+
+  allocate(loadCases(load_steps%length))
   do l = 1, load_steps%length
-
     load_step => load_steps%get_dict(l)
     step_bc   => load_step%get_dict('boundary_conditions')
     step_mech => step_bc%get_dict('mechanical')
@@ -316,226 +547,6 @@ program DAMASK_grid
 
     end if reportAndCheck
   end do
-
+end function parseLoadsteps
-!--------------------------------------------------------------------------------------------------
-! doing initialization depending on active solvers
-  call spectral_utilities_init()
-
-  do field = 2, nActiveFields
-    select case (ID(field))
-
-      case (FIELD_THERMAL_ID)
-        call grid_thermal_spectral_init(num_grid)
-
-      case (FIELD_DAMAGE_ID)
-        call grid_damage_spectral_init(num_grid)
-
-    end select
-  end do
-
-  call mechanical_init(num_grid)
-  call config_numerics_deallocate()
-
-!--------------------------------------------------------------------------------------------------
-! write header of output file
-  if (worldrank == 0) then
-    writeHeader: if (CLI_restartInc < 1) then
-      open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
-      write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                ! statistics file
-    else writeHeader
-      open(newunit=statUnit,file=trim(getSolverJobName())//&
-                                  '.sta',form='FORMATTED', position='APPEND', status='OLD')
-    end if writeHeader
-  end if
-
-  writeUndeformed: if (CLI_restartInc < 1) then
-    print'(/,1x,a)', '... saving initial configuration ..........................................'
-    flush(IO_STDOUT)
-    call materialpoint_result(0,0.0_pREAL)
-  end if writeUndeformed
-
-  loadCaseLooping: do l = 1, size(loadCases)
-    t_0 = t                                                                                         ! load case start time
-    guess = loadCases(l)%estimate_rate                                                              ! change of load case? homogeneous guess for the first inc
-
-    incLooping: do inc = 1, loadCases(l)%N
-      totalIncsCounter = totalIncsCounter + 1
-
-!--------------------------------------------------------------------------------------------------
-! forwarding time
-      Delta_t_prev = Delta_t                                                                        ! last time intervall that brought former inc to an end
-      if (dEq(loadCases(l)%r,1.0_pREAL,1.e-9_pREAL)) then                                           ! linear scale
-        Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL)
-      else
-        Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) &
-                                 / (1.0_pREAL-loadCases(l)%r**loadCases(l)%N)
-      end if
-      Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL)                      ! depending on cut back level, decrease time step
-
-      skipping: if (totalIncsCounter <= CLI_restartInc) then                                        ! not yet at restart inc?
-        t = t + Delta_t                                                                             ! just advance time, skip already performed calculation
-        guess = .true.                                                                              ! QUESTION:why forced guessing instead of inheriting loadcase preference
-      else skipping
-        stepFraction = 0                                                                            ! fraction scaled by stepFactor**cutLevel
-
-        subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
-          t_remaining = loadCases(l)%t + t_0 - t
-          t = t + Delta_t                                                                           ! forward target time
-          stepFraction = stepFraction + 1                                                           ! count step
-
-!--------------------------------------------------------------------------------------------------
-! report beginning of new step
-          print'(/,1x,a)', '###########################################################################'
-          print'(1x,a,1x,es12.5,6(a,i0))', &
-                  'Time', t, &
-                  's: Increment ', inc,'/',loadCases(l)%N,&
-                  '-', stepFraction,'/',subStepFactor**cutBackLevel,&
-                  ' of load case ', l,'/',size(loadCases)
-          write(incInfo,'(4(a,i0))') &
-                  'Increment ',totalIncsCounter,'/',sum(loadCases%N),&
-                  '-', stepFraction,'/',subStepFactor**cutBackLevel
-          flush(IO_STDOUT)
-
-!--------------------------------------------------------------------------------------------------
-! forward fields
-          do field = 1, nActiveFields
-            select case(ID(field))
-              case(FIELD_MECH_ID)
-                call mechanical_forward (&
-                        cutBack,guess,Delta_t,Delta_t_prev,t_remaining, &
-                        deformation_BC = loadCases(l)%deformation, &
-                        stress_BC      = loadCases(l)%stress, &
-                        rotation_BC    = loadCases(l)%rot)
-
-              case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
-              case(FIELD_DAMAGE_ID);  call grid_damage_spectral_forward(cutBack)
-            end select
-          end do
-          if (.not. cutBack) call materialpoint_forward
-
-!--------------------------------------------------------------------------------------------------
-! solve fields
-          stagIter = 0
-          stagIterate = .true.
-          do while (stagIterate)
-            do field = 1, nActiveFields
-              select case(ID(field))
-                case(FIELD_MECH_ID)
-                  solres(field) = mechanical_solution(incInfo)
-                case(FIELD_THERMAL_ID)
-                  solres(field) = grid_thermal_spectral_solution(Delta_t)
-                case(FIELD_DAMAGE_ID)
-                  solres(field) = grid_damage_spectral_solution(Delta_t)
-              end select
-
-              if (.not. solres(field)%converged) exit                                               ! no solution found
-
-            end do
-            stagIter = stagIter + 1
-            stagIterate =            stagIter < stagItMax &
-                         .and.       all(solres(:)%converged) &
-                         .and. .not. all(solres(:)%stagConverged)                                   ! stationary with respect to staggered iteration
-          end do
-
-!--------------------------------------------------------------------------------------------------
-! check solution and either advance or retry with smaller timestep
-
-          if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) &                           ! converged
-               .and. .not. solres(1)%termIll) then                                                  ! and acceptable solution found
-            call mechanical_updateCoords()
-            Delta_t_prev = Delta_t
-            cutBack = .false.
-            guess = .true.                                                                          ! start guessing after first converged (sub)inc
-            if (worldrank == 0) then
-              write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
-                                solres(1)%converged, solres(1)%iterationsNeeded
-              flush(statUnit)
-            end if
-          elseif (cutBackLevel < maxCutBack) then                                                   ! further cutbacking tolerated?
-            cutBack = .true.
-            stepFraction = (stepFraction - 1) * subStepFactor                                       ! adjust to new denominator
-            cutBackLevel = cutBackLevel + 1
-            t = t - Delta_t
-            Delta_t = Delta_t/real(subStepFactor,pREAL)                                             ! cut timestep
-            print'(/,1x,a)', 'cutting back '
-          else                                                                                      ! no more options to continue
-            if (worldrank == 0) close(statUnit)
-            call IO_error(950)
-          end if
-
-        end do subStepLooping
-
-        cutBackLevel = max(0, cutBackLevel - 1)                                                     ! try half number of subincs next inc
-
-        if (all(solres(:)%converged)) then
-          print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'converged'
-        else
-          print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged'
-        end if; flush(IO_STDOUT)
-
-        call MPI_Allreduce(signal_SIGUSR1,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
-        if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-        if (mod(inc,loadCases(l)%f_out) == 0 .or. sig) then
-          print'(/,1x,a)', '... saving results ........................................................'
-          flush(IO_STDOUT)
-          call materialpoint_result(totalIncsCounter,t)
-        end if
-        if (sig) call signal_setSIGUSR1(.false.)
-        call MPI_Allreduce(signal_SIGUSR2,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
-        if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-        if (mod(inc,loadCases(l)%f_restart) == 0 .or. sig) then
-          do field = 1, nActiveFields
-            select case (ID(field))
-              case(FIELD_MECH_ID)
-                call mechanical_restartWrite()
-              case(FIELD_THERMAL_ID)
-                call grid_thermal_spectral_restartWrite()
-              case(FIELD_DAMAGE_ID)
-                call grid_damage_spectral_restartWrite()
-            end select
-          end do
-          call materialpoint_restartWrite()
-        end if
-        if (sig) call signal_setSIGUSR2(.false.)
-        call MPI_Allreduce(signal_SIGINT,sig,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
-        if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-        if (sig) exit loadCaseLooping
-      end if skipping
-
-    end do incLooping
-
-  end do loadCaseLooping
-
-
-!--------------------------------------------------------------------------------------------------
-! report summary of whole calculation
-  print'(/,1x,a)', '###########################################################################'
-  if (worldrank == 0) close(statUnit)
-
-  call quit(0)                                                                                      ! no complains ;)
-
-
-contains
-
-subroutine getMaskedTensor(values,mask,tensor)
-
-  real(pREAL), intent(out), dimension(3,3) :: values
-  logical,     intent(out), dimension(3,3) :: mask
-  type(tList), pointer :: tensor
-
-  type(tList), pointer :: row
-  integer :: i,j
-
-
-  values = 0.0_pREAL
-  do i = 1,3
-    row => tensor%get_list(i)
-    do j = 1,3
-      mask(i,j) = row%get_asStr(j) == 'x'
-      if (.not. mask(i,j)) values(i,j) = row%get_asReal(j)
-    end do
-  end do
-
-end subroutine getMaskedTensor
 
 end program DAMASK_grid

src/grid/grid_damage_spectral.f90

@@ -75,7 +75,6 @@ subroutine grid_damage_spectral_init(num_grid)
   type(tDict), pointer, intent(in) :: num_grid
 
   integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
-  integer :: i, j, k, ce
   DM :: DM_damage
   real(pREAL), dimension(:,:,:), pointer :: phi                                                     ! 0-indexed
   Vec :: uBound, lBound
@@ -92,7 +91,6 @@ subroutine grid_damage_spectral_init(num_grid)
     petsc_options
 
 
-
   print'(/,1x,a)', '<<<+-  grid_spectral_damage init  -+>>>'
 
   print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
@@ -194,11 +192,7 @@ subroutine grid_damage_spectral_init(num_grid)
     phi_stagInc = phi_lastInc
   end if restartRead
 
-  ce = 0
+  call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3]))
-  do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1
-    ce = ce + 1
-    call homogenization_set_phi(phi(i,j,k),ce)
-  end do; end do; end do
 
   call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -216,7 +210,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   real(pREAL), intent(in) :: &
     Delta_t                                                                                         !< increment in time for current solution
 
-  integer :: i, j, k, ce
   type(tSolutionState) :: solution
   PetscInt  :: devNull
   PetscReal :: phi_min, phi_max, stagNorm
@@ -227,8 +220,6 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   SNESConvergedReason :: reason
 
 
-  solution%converged = .false.
-
 !--------------------------------------------------------------------------------------------------
 ! set module wide availabe data
   params%Delta_t = Delta_t
@@ -256,13 +247,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
   phi_stagInc = phi
 
-!--------------------------------------------------------------------------------------------------
+  call homogenization_set_phi(reshape(phi,[product(cells(1:2))*cells3]))
-! updating damage state
-  ce = 0
-  do k = 0, cells3-1;  do j = 0, cells(2)-1;  do i = 0,cells(1)-1
-    ce = ce + 1
-    call homogenization_set_phi(phi(i,j,k),ce)
-  end do; end do; end do
 
   call DMDAVecRestoreArrayF90(DM_damage,phi_PETSc,phi,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -283,7 +268,6 @@ subroutine grid_damage_spectral_forward(cutBack)
 
   logical, intent(in) :: cutBack
 
-  integer :: i, j, k, ce
   DM :: DM_damage
   real(pREAL),  dimension(:,:,:), pointer :: phi                                                    ! 0-indexed
   PetscErrorCode :: err_PETSc
@@ -295,11 +279,7 @@ subroutine grid_damage_spectral_forward(cutBack)
     CHKERRQ(err_PETSc)
 
   if (cutBack) then
-    ce = 0
+    call homogenization_set_phi(reshape(phi_lastInc,[product(cells(1:2))*cells3]))
-    do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells(1)
-      ce = ce + 1
-      call homogenization_set_phi(phi_lastInc(i,j,k),ce)
-    end do; end do; end do
     phi = phi_lastInc
     phi_stagInc = phi_lastInc
   else

src/grid/grid_mech_spectral_basic.f90

@@ -365,7 +365,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
   if (stress_BC%myType=='dot_P') P_aim = P_aim &
                                        + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t
 
-  F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, &                                      ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
+  F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, &                                ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
               rotation_BC%rotate(F_aim,active=.true.)),[9,cells(1),cells(2),cells3])
   call DMDAVecRestoreArrayF90(DM_mech,F_PETSc,F,err_PETSc)
   CHKERRQ(err_PETSc)

src/grid/grid_mech_spectral_polarization.f90

@@ -408,11 +408,11 @@ subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,D
   if (stress_BC%myType=='dot_P') P_aim = P_aim &
                                       + merge(.0_pREAL,stress_BC%values,stress_BC%mask)*Delta_t
 
-  F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, &                                      ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
+  F = reshape(utilities_forwardTensorField(Delta_t,F_lastInc,Fdot, &                                ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
                                            rotation_BC%rotate(F_aim,active=.true.)),&
               [9,cells(1),cells(2),cells3])
   if (guess) then
-     F_tau = reshape(Utilities_forwardField(Delta_t,F_tau_lastInc,F_taudot), &
+     F_tau = reshape(Utilities_forwardTensorField(Delta_t,F_tau_lastInc,F_taudot), &
                      [9,cells(1),cells(2),cells3])                                                  ! does not have any average value as boundary condition
    else
     do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1)

src/grid/grid_thermal_spectral.f90

@@ -74,7 +74,7 @@ subroutine grid_thermal_spectral_init(num_grid)
   type(tDict), pointer, intent(in) :: num_grid
 
   integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: cells3_global
-  integer :: i, j, k, ce
+  integer :: ce
   DM :: DM_thermal
   real(pREAL), dimension(:,:,:), pointer :: T                                                       ! 0-indexed
   integer(MPI_INTEGER_KIND) :: err_MPI
@@ -88,7 +88,6 @@ subroutine grid_thermal_spectral_init(num_grid)
     petsc_options
 
 
-
   print'(/,1x,a)', '<<<+-  grid_thermal_spectral init  -+>>>'
 
   print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
@@ -171,11 +170,8 @@ subroutine grid_thermal_spectral_init(num_grid)
     dotT_lastInc = 0.0_pREAL * T_lastInc
   end if restartRead
 
-  ce = 0
+  call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), &
-  do k = 0, cells3-1; do j = 0, cells(2)-1; do i = 0, cells(1)-1
+                                       [(0.0_pReal, ce = 1,product(cells(1:2))*cells3)])
-    ce = ce + 1
-    call homogenization_thermal_setField(T(i,j,k),0.0_pREAL,ce)
-  end do; end do; end do
 
   call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -204,8 +200,6 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   SNESConvergedReason :: reason
 
 
-  solution%converged = .false.
-
 !--------------------------------------------------------------------------------------------------
 ! set module wide availabe data
   params%Delta_t = Delta_t
@@ -233,13 +227,8 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
   T_stagInc = T
 
-!--------------------------------------------------------------------------------------------------
+  call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), &
-! updating thermal state
+                                       reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t)
-  ce = 0
-  do k = 0, cells3-1;  do j = 0, cells(2)-1;  do i = 0, cells(1)-1
-    ce = ce + 1
-    call homogenization_thermal_setField(T(i,j,k),(T(i,j,k)-T_lastInc(i+1,j+1,k+1))/params%Delta_t,ce)
-  end do; end do; end do
 
   call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -260,7 +249,6 @@ subroutine grid_thermal_spectral_forward(cutBack)
 
   logical, intent(in) :: cutBack
 
-  integer :: i, j, k, ce
   DM :: DM_thermal
   real(pREAL),  dimension(:,:,:), pointer :: T                                                      ! 0-indexed
   PetscErrorCode :: err_PETSc
@@ -272,11 +260,8 @@ subroutine grid_thermal_spectral_forward(cutBack)
   CHKERRQ(err_PETSc)
 
   if (cutBack) then
-    ce = 0
+    call homogenization_thermal_setField(reshape(T_lastInc,[product(cells(1:2))*cells3]), &
-    do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells(1)
+                                         reshape(dotT_lastInc,[product(cells(1:2))*cells3]))
-      ce = ce + 1
-      call homogenization_thermal_setField(T_lastInc(i,j,k),dotT_lastInc(i,j,k),ce)
-    end do; end do; end do
     T = T_lastInc
     T_stagInc = T_lastInc
   else

src/grid/spectral_utilities.f90

@@ -124,7 +124,7 @@ module spectral_utilities
     utilities_maskedCompliance, &
     utilities_constitutiveResponse, &
     utilities_calculateRate, &
-    utilities_forwardField, &
+    utilities_forwardTensorField, &
     utilities_updateCoords
 
 contains
@@ -864,7 +864,7 @@ end function utilities_calculateRate
 !> @brief forwards a field with a pointwise given rate, if aim is given,
 !> ensures that the average matches the aim
 !--------------------------------------------------------------------------------------------------
-function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
+function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim)
 
   real(pREAL), intent(in) :: &
     Delta_t                                                                                         !< Delta_t of current step
@@ -875,22 +875,22 @@ function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
     aim                                                                                             !< average field value aim
 
   real(pREAL),                       dimension(3,3,cells(1),cells(2),cells3) :: &
-    utilities_forwardField
+    utilities_forwardTensorField
   real(pREAL),                       dimension(3,3) :: fieldDiff                                    !< <a + adot*t> - aim
   integer(MPI_INTEGER_KIND) :: err_MPI
 
 
-  utilities_forwardField = field_lastInc + rate*Delta_t
+  utilities_forwardTensorField = field_lastInc + rate*Delta_t
   if (present(aim)) then                                                                            !< correct to match average
-    fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt
+    fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt
     call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
     if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
     fieldDiff = fieldDiff - aim
-    utilities_forwardField = utilities_forwardField &
+    utilities_forwardTensorField = utilities_forwardTensorField &
                                  - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
   end if
 
-end function utilities_forwardField
+end function utilities_forwardTensorField
 
 
 !--------------------------------------------------------------------------------------------------

src/homogenization.f90

@@ -145,9 +145,8 @@ module homogenization
       real(pREAL) :: f
     end function homogenization_f_T
 
-    module subroutine homogenization_thermal_setField(T,dot_T, ce)
+    module subroutine homogenization_thermal_setField(T,dot_T)
-      integer, intent(in) :: ce
+      real(pREAL), dimension(:),  intent(in) :: T, dot_T
-      real(pREAL),   intent(in) :: T, dot_T
     end subroutine homogenization_thermal_setField
 
     module function homogenization_damage_active() result(active)
@@ -170,10 +169,8 @@ module homogenization
       real(pREAL) :: f
     end function homogenization_f_phi
 
-    module subroutine homogenization_set_phi(phi,ce)
+    module subroutine homogenization_set_phi(phi)
-      integer, intent(in) :: ce
+      real(pREAL), dimension(:), intent(in) :: phi
-      real(pREAL),   intent(in) :: &
-        phi
     end subroutine homogenization_set_phi
 
   end interface

src/homogenization_damage.f90

@@ -151,20 +151,19 @@ end function homogenization_f_phi
 !--------------------------------------------------------------------------------------------------
 !> @brief Set damage field.
 !--------------------------------------------------------------------------------------------------
-module subroutine homogenization_set_phi(phi,ce)
+module subroutine homogenization_set_phi(phi)
 
-  integer, intent(in) :: ce
+  real(pREAL), dimension(:), intent(in) :: phi
-  real(pREAL), intent(in) :: phi
 
-  integer :: &
+  integer :: ho, en, ce
-    ho, &
-    en
 
 
+  do ce=lbound(phi,1), ubound(phi,1)
     ho = material_ID_homogenization(ce)
     en = material_entry_homogenization(ce)
-  damagestate_h(ho)%state(1,en) = phi
+    damagestate_h(ho)%state(1,en) = phi(ce)
-  current(ho)%phi(en) = phi
+    current(ho)%phi(en) = phi(ce)
+  end do
 
 end subroutine homogenization_set_phi
 

src/homogenization_thermal.f90

@@ -173,15 +173,20 @@ end function homogenization_f_T
 !--------------------------------------------------------------------------------------------------
 !> @brief Set thermal field and its rate (T and dot_T).
 !--------------------------------------------------------------------------------------------------
-module subroutine homogenization_thermal_setField(T,dot_T, ce)
+module subroutine homogenization_thermal_setField(T,dot_T)
 
-  integer, intent(in) :: ce
+  real(pREAL), dimension(:), intent(in) :: T, dot_T
-  real(pREAL), intent(in) :: T, dot_T
+
+  integer :: ho, en, ce
 
 
-  current(material_ID_homogenization(ce))%T(material_entry_homogenization(ce)) = T
+  do ce=max(lbound(T,1),lbound(dot_T,1)), min(ubound(T,1),ubound(dot_T,1))
-  current(material_ID_homogenization(ce))%dot_T(material_entry_homogenization(ce)) = dot_T
+    ho = material_ID_homogenization(ce)
+    en = material_entry_homogenization(ce)
+    current(ho)%T(en) = T(ce)
+    current(ho)%dot_T(en) = dot_T(ce)
     call thermal_partition(ce)
+  end do
 
 end subroutine homogenization_thermal_setField
 

src/phase.f90

@@ -49,6 +49,29 @@ module phase
     type(tState), dimension(:), allocatable :: p                                                    !< tState for each active source mechanism in a phase
   end type
 
+  enum, bind(c); enumerator :: &
+    UNDEFINED, &
+    MECHANICAL_PLASTICITY_NONE, &
+    MECHANICAL_PLASTICITY_ISOTROPIC, &
+    MECHANICAL_PLASTICITY_PHENOPOWERLAW, &
+    MECHANICAL_PLASTICITY_KINEHARDENING, &
+    MECHANICAL_PLASTICITY_DISLOTWIN, &
+    MECHANICAL_PLASTICITY_DISLOTUNGSTEN, &
+    MECHANICAL_PLASTICITY_NONLOCAL, &
+    MECHANICAL_EIGEN_THERMALEXPANSION, &
+    DAMAGE_ISOBRITTLE, &
+    DAMAGE_ANISOBRITTLE, &
+    THERMAL_SOURCE_DISSIPATION, &
+    THERMAL_SOURCE_EXTERNALHEAT
+  end enum
+
+
+  integer(kind(UNDEFINED)), dimension(:), allocatable :: &
+    mechanical_plasticity_type, &                                                                   !< plasticity of each phase
+    damage_type                                                                                     !< damage type of each phase
+  integer(kind(UNDEFINED)),  dimension(:,:), allocatable :: &
+    thermal_source_type, &
+    mechanical_eigen_kinematics_type
 
   character(len=2), allocatable, dimension(:) :: phase_lattice
   real(pREAL),      allocatable, dimension(:) :: phase_cOverA

src/phase_damage.f90

@@ -9,23 +9,14 @@ submodule(phase) damage
       l_c = 0.0_pREAL                                                                               !< characteristic length
   end type tDamageParameters
 
-  enum, bind(c); enumerator :: &
-    DAMAGE_UNDEFINED_ID, &
-    DAMAGE_ISOBRITTLE_ID, &
-    DAMAGE_ANISOBRITTLE_ID
-  end enum
-
   integer :: phase_damage_maxSizeDotState
 
-
+  type :: tFieldQuantities
-  type :: tDataContainer
     real(pREAL), dimension(:), allocatable :: phi
-  end type tDataContainer
+  end type tFieldQuantities
 
-  integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
-    phase_damage                                                                                    !< active sources mechanisms of each phase
 
-  type(tDataContainer), dimension(:), allocatable :: current
+  type(tFieldQuantities), dimension(:), allocatable :: current
 
   type(tDamageParameters), dimension(:), allocatable :: param
 
@@ -114,11 +105,11 @@ module subroutine damage_init()
 
   end do
 
-  allocate(phase_damage(phases%length), source = DAMAGE_UNDEFINED_ID)
+  allocate(damage_type(phases%length), source = UNDEFINED)
 
   if (damage_active) then
-    where(isobrittle_init()  ) phase_damage = DAMAGE_ISOBRITTLE_ID
+    where(isobrittle_init()  ) damage_type = DAMAGE_ISOBRITTLE
-    where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
+    where(anisobrittle_init()) damage_type = DAMAGE_ANISOBRITTLE
   end if
 
   phase_damage_maxSizeDotState = maxval(damageState%sizeDotState)
@@ -159,8 +150,8 @@ module function phase_damage_C66(C66,ph,en) result(C66_degraded)
   real(pREAL), dimension(6,6) :: C66_degraded
 
 
-  damageType: select case (phase_damage(ph))
+  damageType: select case (damage_type(ph))
-    case (DAMAGE_ISOBRITTLE_ID) damageType
+    case (DAMAGE_ISOBRITTLE) damageType
       C66_degraded = C66 * damage_phi(ph,en)**2
     case default damageType
       C66_degraded = C66
@@ -204,13 +195,14 @@ module function phase_f_phi(phi,co,ce) result(f)
     ph, &
     en
 
+
   ph = material_ID_phase(co,ce)
   en = material_entry_phase(co,ce)
 
-  select case(phase_damage(ph))
+  select case(damage_type(ph))
-    case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID)
+    case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE)
       f = 1.0_pREAL &
-        - 2.0_pREAL * phi*damageState(ph)%state(1,en)
+        - 2.0_pREAL * phi*damageState(ph)%state(1,en)                                               ! ToDo: MD: seems to be phi**2
     case default
       f = 0.0_pREAL
   end select
@@ -318,8 +310,8 @@ module subroutine damage_restartWrite(groupHandle,ph)
   integer, intent(in) :: ph
 
 
-  select case(phase_damage(ph))
+  select case(damage_type(ph))
-    case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID)
+    case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE)
       call HDF5_write(damageState(ph)%state,groupHandle,'omega_damage')
   end select
 
@@ -332,8 +324,8 @@ module subroutine damage_restartRead(groupHandle,ph)
   integer, intent(in) :: ph
 
 
-  select case(phase_damage(ph))
+  select case(damage_type(ph))
-    case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID)
+    case(DAMAGE_ISOBRITTLE,DAMAGE_ANISOBRITTLE)
   call HDF5_read(damageState(ph)%state0,groupHandle,'omega_damage')
   end select
 
@@ -350,15 +342,15 @@ module subroutine damage_result(group,ph)
   integer,          intent(in) :: ph
 
 
-  if (phase_damage(ph) /= DAMAGE_UNDEFINED_ID) &
+  if (damage_type(ph) /= UNDEFINED) &
     call result_closeGroup(result_addGroup(group//'damage'))
 
-  sourceType: select case (phase_damage(ph))
+  sourceType: select case (damage_type(ph))
 
-    case (DAMAGE_ISOBRITTLE_ID) sourceType
+    case (DAMAGE_ISOBRITTLE) sourceType
       call isobrittle_result(ph,group//'damage/')
 
-    case (DAMAGE_ANISOBRITTLE_ID) sourceType
+    case (DAMAGE_ANISOBRITTLE) sourceType
       call anisobrittle_result(ph,group//'damage/')
 
   end select sourceType
@@ -381,9 +373,9 @@ function phase_damage_collectDotState(ph,en) result(broken)
 
   if (damageState(ph)%sizeState > 0) then
 
-    sourceType: select case (phase_damage(ph))
+    sourceType: select case (damage_type(ph))
 
-      case (DAMAGE_ANISOBRITTLE_ID) sourceType
+      case (DAMAGE_ANISOBRITTLE) sourceType
         call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! ToDo: use M_d
 
     end select sourceType
@@ -446,9 +438,9 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
 
   if (damageState(ph)%sizeState == 0) return
 
-   sourceType: select case (phase_damage(ph))
+   sourceType: select case (damage_type(ph))
 
-    case (DAMAGE_ISOBRITTLE_ID) sourceType
+    case (DAMAGE_ISOBRITTLE) sourceType
       call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en)
       broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
       if (.not. broken) then

src/phase_mechanical.f90

@@ -3,21 +3,6 @@
 !----------------------------------------------------------------------------------------------------
 submodule(phase) mechanical
 
-
-  enum, bind(c); enumerator :: &
-    PLASTIC_UNDEFINED_ID, &
-    PLASTIC_NONE_ID, &
-    PLASTIC_ISOTROPIC_ID, &
-    PLASTIC_PHENOPOWERLAW_ID, &
-    PLASTIC_KINEHARDENING_ID, &
-    PLASTIC_DISLOTWIN_ID, &
-    PLASTIC_DISLOTUNGSTEN_ID, &
-    PLASTIC_NONLOCAL_ID, &
-    EIGEN_UNDEFINED_ID, &
-    EIGEN_CLEAVAGE_OPENING_ID, &
-    EIGEN_THERMAL_EXPANSION_ID
-  end enum
-
   type(tTensorContainer), dimension(:), allocatable :: &
     ! current value
     phase_mechanical_Fe, &
@@ -37,9 +22,6 @@ submodule(phase) mechanical
     phase_mechanical_S0
 
 
-  integer(kind(PLASTIC_undefined_ID)), dimension(:),   allocatable :: &
-    phase_plasticity                                                                                !< plasticity of each phase
-
   interface
 
     module subroutine eigen_init(phases)
@@ -286,7 +268,7 @@ module subroutine mechanical_init(phases, num_mech)
   call elastic_init(phases)
 
   allocate(plasticState(phases%length))
-  allocate(phase_plasticity(phases%length),source = PLASTIC_UNDEFINED_ID)
+  allocate(mechanical_plasticity_type(phases%length),source = UNDEFINED)
   call plastic_init()
   do ph = 1,phases%length
     plasticState(ph)%state0 = plasticState(ph)%state
@@ -363,24 +345,24 @@ module subroutine mechanical_result(group,ph)
 
   call results(group,ph)
 
-  select case(phase_plasticity(ph))
+  select case(mechanical_plasticity_type(ph))
 
-    case(PLASTIC_ISOTROPIC_ID)
+    case(MECHANICAL_PLASTICITY_ISOTROPIC)
       call plastic_isotropic_result(ph,group//'mechanical/')
 
-    case(PLASTIC_PHENOPOWERLAW_ID)
+    case(MECHANICAL_PLASTICITY_PHENOPOWERLAW)
       call plastic_phenopowerlaw_result(ph,group//'mechanical/')
 
-    case(PLASTIC_KINEHARDENING_ID)
+    case(MECHANICAL_PLASTICITY_KINEHARDENING)
       call plastic_kinehardening_result(ph,group//'mechanical/')
 
-    case(PLASTIC_DISLOTWIN_ID)
+    case(MECHANICAL_PLASTICITY_DISLOTWIN)
       call plastic_dislotwin_result(ph,group//'mechanical/')
 
-    case(PLASTIC_DISLOTUNGSTEN_ID)
+    case(MECHANICAL_PLASTICITY_DISLOTUNGSTEN)
       call plastic_dislotungsten_result(ph,group//'mechanical/')
 
-    case(PLASTIC_NONLOCAL_ID)
+    case(MECHANICAL_PLASTICITY_NONLOCAL)
       call plastic_nonlocal_result(ph,group//'mechanical/')
 
   end select

src/phase_mechanical_eigen.f90

@@ -3,15 +3,7 @@ submodule(phase:mechanical) eigen
   integer, dimension(:), allocatable :: &
     Nmodels
 
-  integer(kind(EIGEN_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
-    model
-  integer(kind(EIGEN_UNDEFINED_ID)),  dimension(:), allocatable :: &
-    model_damage
-
   interface
-    module function damage_anisobrittle_init() result(myKinematics)
-      logical, dimension(:), allocatable :: myKinematics
-    end function damage_anisobrittle_init
 
     module function thermalexpansion_init(kinematics_length) result(myKinematics)
       integer, intent(in) :: kinematics_length
@@ -60,17 +52,12 @@ module subroutine eigen_init(phases)
     Nmodels(ph) = kinematics%length
   end do
 
-  allocate(model(maxval(Nmodels),phases%length), source = EIGEN_undefined_ID)
+  allocate(mechanical_eigen_kinematics_type(maxval(Nmodels),phases%length), source = UNDEFINED)
 
   if (maxval(Nmodels) /= 0) then
-    where(thermalexpansion_init(maxval(Nmodels))) model = EIGEN_thermal_expansion_ID
+    where(thermalexpansion_init(maxval(Nmodels))) mechanical_eigen_kinematics_type = MECHANICAL_EIGEN_THERMALEXPANSION
   end if
 
-  allocate(model_damage(phases%length),  source = EIGEN_UNDEFINED_ID)
-
-  where(damage_anisobrittle_init())  model_damage = EIGEN_cleavage_opening_ID
-
-
 end subroutine eigen_init
 
 
@@ -173,17 +160,9 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
   dLi_dFi = 0.0_pREAL
 
 
-  plasticType: select case (phase_plasticity(ph))
-    case (PLASTIC_isotropic_ID) plasticType
-      call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en)
-      Li = Li + my_Li
-      dLi_dS = dLi_dS + my_dLi_dS
-      active = .true.
-  end select plasticType
-
   KinematicsLoop: do k = 1, Nmodels(ph)
-    kinematicsType: select case (model(k,ph))
+    kinematicsType: select case (mechanical_eigen_kinematics_type(k,ph))
-      case (EIGEN_thermal_expansion_ID) kinematicsType
+      case (MECHANICAL_EIGEN_THERMALEXPANSION) kinematicsType
         call thermalexpansion_LiAndItsTangent(my_Li, my_dLi_dS, ph,en)
         Li = Li + my_Li
         dLi_dS = dLi_dS + my_dLi_dS
@@ -191,13 +170,21 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
     end select kinematicsType
   end do KinematicsLoop
 
-  select case (model_damage(ph))
+  plasticType: select case (mechanical_plasticity_type(ph))
-    case (EIGEN_cleavage_opening_ID)
+    case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType
+      call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,ph,en)
+      Li = Li + my_Li
+      dLi_dS = dLi_dS + my_dLi_dS
+      active = .true.
+  end select plasticType
+
+  damageType: select case (damage_type(ph))
+    case (DAMAGE_ANISOBRITTLE)
       call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en)
       Li = Li + my_Li
       dLi_dS = dLi_dS + my_dLi_dS
       active = .true.
-  end select
+  end select damageType
 
   if (.not. active) return
 

src/phase_mechanical_eigen_cleavageopening.f90

@@ -1,30 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-submodule(phase:eigen) cleavageopening
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-module function damage_anisobrittle_init() result(myKinematics)
-
-  logical, dimension(:), allocatable :: myKinematics
-
-
-  myKinematics = kinematics_active2('anisobrittle')
-  if (count(myKinematics) == 0) return
-
-  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen:cleavageopening init  -+>>>'
-  print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
-
-end function damage_anisobrittle_init
-
-
-end submodule cleavageopening

src/phase_mechanical_elastic.f90

@@ -199,8 +199,8 @@ module function phase_homogenizedC66(ph,en) result(C)
   integer,      intent(in)    :: ph, en
 
 
-  plasticType: select case (phase_plasticity(ph))
+  plasticType: select case (mechanical_plasticity_type(ph))
-    case (PLASTIC_DISLOTWIN_ID) plasticType
+    case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
      C = plastic_dislotwin_homogenizedC(ph,en)
     case default plasticType
      C = elastic_C66(ph,en)

src/phase_mechanical_plastic.f90

@@ -211,17 +211,17 @@ contains
 module subroutine plastic_init
 
 
-  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plasticity init  -+>>>'
 
-  where(plastic_none_init())              phase_plasticity = PLASTIC_NONE_ID
+  where(plastic_none_init())              mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONE
-  where(plastic_isotropic_init())         phase_plasticity = PLASTIC_ISOTROPIC_ID
+  where(plastic_isotropic_init())         mechanical_plasticity_type = MECHANICAL_PLASTICITY_ISOTROPIC
-  where(plastic_phenopowerlaw_init())     phase_plasticity = PLASTIC_PHENOPOWERLAW_ID
+  where(plastic_phenopowerlaw_init())     mechanical_plasticity_type = MECHANICAL_PLASTICITY_PHENOPOWERLAW
-  where(plastic_kinehardening_init())     phase_plasticity = PLASTIC_KINEHARDENING_ID
+  where(plastic_kinehardening_init())     mechanical_plasticity_type = MECHANICAL_PLASTICITY_KINEHARDENING
-  where(plastic_dislotwin_init())         phase_plasticity = PLASTIC_DISLOTWIN_ID
+  where(plastic_dislotwin_init())         mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTWIN
-  where(plastic_dislotungsten_init())     phase_plasticity = PLASTIC_DISLOTUNGSTEN_ID
+  where(plastic_dislotungsten_init())     mechanical_plasticity_type = MECHANICAL_PLASTICITY_DISLOTUNGSTEN
-  where(plastic_nonlocal_init())          phase_plasticity = PLASTIC_NONLOCAL_ID
+  where(plastic_nonlocal_init())          mechanical_plasticity_type = MECHANICAL_PLASTICITY_NONLOCAL
 
-  if (any(phase_plasticity == PLASTIC_undefined_ID)) call IO_error(201)
+  if (any(mechanical_plasticity_type == UNDEFINED)) call IO_error(201)
 
 end subroutine plastic_init
 
@@ -251,7 +251,7 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
     i, j
 
 
-  if (phase_plasticity(ph) == PLASTIC_NONE_ID) then
+  if (mechanical_plasticity_type(ph) == MECHANICAL_PLASTICITY_NONE) then
     Lp      = 0.0_pREAL
     dLp_dFi = 0.0_pREAL
     dLp_dS  = 0.0_pREAL
@@ -259,24 +259,24 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
 
     Mp = matmul(matmul(transpose(Fi),Fi),S)
 
-    plasticType: select case (phase_plasticity(ph))
+    plasticType: select case (mechanical_plasticity_type(ph))
 
-      case (PLASTIC_ISOTROPIC_ID) plasticType
+      case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType
         call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTIC_PHENOPOWERLAW_ID) plasticType
+      case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType
         call phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTIC_KINEHARDENING_ID) plasticType
+      case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType
         call kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTIC_NONLOCAL_ID) plasticType
+      case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
         call nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTIC_DISLOTWIN_ID) plasticType
+      case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
         call dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
-      case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
+      case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType
         call dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
     end select plasticType
@@ -308,28 +308,28 @@ module function plastic_dotState(subdt,ph,en) result(dotState)
     dotState
 
 
-  if (phase_plasticity(ph) /= PLASTIC_NONE_ID) then
+  if (mechanical_plasticity_type(ph) /= MECHANICAL_PLASTICITY_NONE) then
     Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                        phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en))
 
-    plasticType: select case (phase_plasticity(ph))
+    plasticType: select case (mechanical_plasticity_type(ph))
 
-      case (PLASTIC_ISOTROPIC_ID) plasticType
+      case (MECHANICAL_PLASTICITY_ISOTROPIC) plasticType
         dotState = isotropic_dotState(Mp,ph,en)
 
-      case (PLASTIC_PHENOPOWERLAW_ID) plasticType
+      case (MECHANICAL_PLASTICITY_PHENOPOWERLAW) plasticType
         dotState = phenopowerlaw_dotState(Mp,ph,en)
 
-      case (PLASTIC_KINEHARDENING_ID) plasticType
+      case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType
         dotState = plastic_kinehardening_dotState(Mp,ph,en)
 
-      case (PLASTIC_DISLOTWIN_ID) plasticType
+      case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
         dotState = dislotwin_dotState(Mp,ph,en)
 
-      case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
+      case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType
         dotState = dislotungsten_dotState(Mp,ph,en)
 
-      case (PLASTIC_NONLOCAL_ID) plasticType
+      case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
         call nonlocal_dotState(Mp,subdt,ph,en)
         dotState = plasticState(ph)%dotState(:,en)
 
@@ -349,15 +349,15 @@ module subroutine plastic_dependentState(ph,en)
     en
 
 
-  plasticType: select case (phase_plasticity(ph))
+  plasticType: select case (mechanical_plasticity_type(ph))
 
-    case (PLASTIC_DISLOTWIN_ID) plasticType
+    case (MECHANICAL_PLASTICITY_DISLOTWIN) plasticType
       call dislotwin_dependentState(ph,en)
 
-    case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
+    case (MECHANICAL_PLASTICITY_DISLOTUNGSTEN) plasticType
       call dislotungsten_dependentState(ph,en)
 
-    case (PLASTIC_NONLOCAL_ID) plasticType
+    case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
       call nonlocal_dependentState(ph,en)
 
   end select plasticType
@@ -384,19 +384,19 @@ module function plastic_deltaState(ph, en) result(broken)
 
   broken = .false.
 
-  select case (phase_plasticity(ph))
+  select case (mechanical_plasticity_type(ph))
-    case (PLASTIC_NONLOCAL_ID,PLASTIC_KINEHARDENING_ID)
+    case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING)
 
       Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                          phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                   phase_mechanical_S(ph)%data(1:3,1:3,en))
 
-      plasticType: select case (phase_plasticity(ph))
+      plasticType: select case (mechanical_plasticity_type(ph))
 
-        case (PLASTIC_KINEHARDENING_ID) plasticType
+        case (MECHANICAL_PLASTICITY_KINEHARDENING) plasticType
           call plastic_kinehardening_deltaState(Mp,ph,en)
 
-        case (PLASTIC_NONLOCAL_ID) plasticType
+        case (MECHANICAL_PLASTICITY_NONLOCAL) plasticType
           call plastic_nonlocal_deltaState(Mp,ph,en)
 
       end select plasticType

src/phase_mechanical_plastic_nonlocal.f90

@@ -1252,7 +1252,7 @@ function rhoDotFlux(timestep,ph,en)
       !* The entering flux from my neighbor will be distributed on my slip systems according to the
       !* compatibility
       if (neighbor_n > 0) then
-      if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID .and. &
+      if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL .and. &
           any(dependentState(ph)%compatibility(:,:,:,n,en) > 0.0_pREAL)) then
 
         forall (s = 1:ns, t = 1:4)
@@ -1298,7 +1298,7 @@ function rhoDotFlux(timestep,ph,en)
       !* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
       !* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
       if (opposite_n > 0) then
-      if (phase_plasticity(np) == PLASTIC_NONLOCAL_ID) then
+      if (mechanical_plasticity_type(np) == MECHANICAL_PLASTICITY_NONLOCAL) then
 
         normal_me2neighbor_defConf = math_det33(Favg) &
                                    * matmul(math_inv33(transpose(Favg)),geom(ph)%IPareaNormal(1:3,n,en))  ! normal of the interface in (average) deformed configuration (pointing en => neighbor)

src/phase_thermal.f90

@@ -15,19 +15,11 @@ submodule(phase) thermal
   type(tSourceState),  allocatable, dimension(:) :: &
     thermalState
 
-  enum, bind(c); enumerator :: &
+  type :: tFieldQuantities
-    THERMAL_UNDEFINED_ID ,&
-    THERMAL_DISSIPATION_ID, &
-    THERMAL_EXTERNALHEAT_ID
-  end enum
-
-  type :: tDataContainer             ! ?? not very telling name. Better: "fieldQuantities" ??
     real(pREAL), dimension(:), allocatable :: T, dot_T
-  end type tDataContainer
+  end type tFieldQuantities
-  integer(kind(THERMAL_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
-    thermal_source
 
-  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(en) reads quite good
+  type(tFieldQuantities), dimension(:), allocatable :: current
 
   type(tThermalParameters), dimension(:), allocatable :: param
 
@@ -36,36 +28,36 @@ submodule(phase) thermal
 
   interface
 
-    module function dissipation_init(source_length) result(mySources)
+    module function source_dissipation_init(maxNsources) result(isMySource)
-      integer, intent(in) :: source_length
+      integer, intent(in) :: maxNsources
-      logical, dimension(:,:), allocatable :: mySources
+      logical, dimension(:,:), allocatable :: isMySource
-    end function dissipation_init
+    end function source_dissipation_init
 
-    module function externalheat_init(source_length) result(mySources)
+    module function source_externalheat_init(maxNsources) result(isMySource)
-      integer, intent(in) :: source_length
+      integer, intent(in) :: maxNsources
-      logical, dimension(:,:), allocatable :: mySources
+      logical, dimension(:,:), allocatable :: isMySource
-    end function externalheat_init
+    end function source_externalheat_init
 
 
-    module subroutine externalheat_dotState(ph, en)
+    module subroutine source_externalheat_dotState(ph, en)
       integer, intent(in) :: &
         ph, &
         en
-    end subroutine externalheat_dotState
+    end subroutine source_externalheat_dotState
 
-    module function dissipation_f_T(ph,en) result(f_T)
+    module function source_dissipation_f_T(ph,en) result(f_T)
       integer, intent(in) :: &
         ph, &
         en
       real(pREAL) :: f_T
-    end function dissipation_f_T
+    end function source_dissipation_f_T
 
-    module function externalheat_f_T(ph,en)  result(f_T)
+    module function source_externalheat_f_T(ph,en)  result(f_T)
       integer, intent(in) :: &
         ph, &
         en
       real(pREAL) :: f_T
-    end function externalheat_f_T
+    end function source_externalheat_f_T
 
  end interface
 
@@ -129,11 +121,11 @@ module subroutine thermal_init(phases)
 
   end do
 
-  allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
+  allocate(thermal_source_type(maxval(thermal_Nsources),phases%length), source = UNDEFINED)
 
   if (maxval(thermal_Nsources) /= 0) then
-    where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
+    where(source_dissipation_init (maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_DISSIPATION
-    where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
+    where(source_externalheat_init(maxval(thermal_Nsources))) thermal_source_type = THERMAL_SOURCE_EXTERNALHEAT
   end if
 
   thermal_source_maxSizeDotState = 0
@@ -151,7 +143,7 @@ end subroutine thermal_init
 
 
 !----------------------------------------------------------------------------------------------
-!< @brief Calculate thermal source.
+!< @brief Calculate thermal source (forcing term).
 !----------------------------------------------------------------------------------------------
 module function phase_f_T(ph,en) result(f)
 
@@ -165,13 +157,13 @@ module function phase_f_T(ph,en) result(f)
   f = 0.0_pREAL
 
   do so = 1, thermal_Nsources(ph)
-   select case(thermal_source(so,ph))
+   select case(thermal_source_type(so,ph))
 
-     case (THERMAL_DISSIPATION_ID)
+     case (THERMAL_SOURCE_DISSIPATION)
-       f = f + dissipation_f_T(ph,en)
+       f = f + source_dissipation_f_T(ph,en)
 
-     case (THERMAL_EXTERNALHEAT_ID)
+     case (THERMAL_SOURCE_EXTERNALHEAT)
-       f = f + externalheat_f_T(ph,en)
+       f = f + source_externalheat_f_T(ph,en)
 
    end select
 
@@ -183,22 +175,22 @@ end function phase_f_T
 !--------------------------------------------------------------------------------------------------
 !> @brief tbd.
 !--------------------------------------------------------------------------------------------------
-function phase_thermal_collectDotState(ph,en) result(broken)
+function phase_thermal_collectDotState(ph,en) result(ok)
 
   integer, intent(in) :: ph, en
-  logical :: broken
+  logical :: ok
 
   integer :: i
 
 
-  broken = .false.
+  ok = .true.
 
   SourceLoop: do i = 1, thermal_Nsources(ph)
 
-    if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) &
+    if (thermal_source_type(i,ph) == THERMAL_SOURCE_EXTERNALHEAT) &
-      call externalheat_dotState(ph,en)
+      call source_externalheat_dotState(ph,en)
 
-    broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
+    ok = ok .and. .not. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,en)))
 
   end do SourceLoop
 
@@ -241,28 +233,27 @@ module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_)
   logical :: converged_
 
 
-  converged_ = .not. integrateThermalState(Delta_t,ph,en)
+  converged_ = integrateThermalState(Delta_t,ph,en)
 
 end function phase_thermal_constitutive
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief integrate state with 1st order explicit Euler method
+!> @brief Integrate state with 1st order explicit Euler method.
 !--------------------------------------------------------------------------------------------------
-function integrateThermalState(Delta_t, ph,en) result(broken)
+function integrateThermalState(Delta_t, ph,en) result(converged)
 
   real(pREAL), intent(in) :: Delta_t
   integer, intent(in) :: ph, en
-  logical :: &
+  logical :: converged
-    broken
 
   integer :: &
     so, &
     sizeDotState
 
 
-  broken = phase_thermal_collectDotState(ph,en)
+  converged = phase_thermal_collectDotState(ph,en)
-  if (broken) return
+  if (converged) then
 
     do so = 1, thermal_Nsources(ph)
       sizeDotState = thermalState(ph)%p(so)%sizeDotState
@@ -270,6 +261,8 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
                                                       + thermalState(ph)%p(so)%dotState(1:sizeDotState,en) * Delta_t
     end do
 
+  end if
+
 end function integrateThermalState
 
 
@@ -318,7 +311,7 @@ end subroutine thermal_forward
 
 
 !----------------------------------------------------------------------------------------------
-!< @brief Get temperature (for use by non-thermal physics)
+!< @brief Get temperature (for use by non-thermal physics).
 !----------------------------------------------------------------------------------------------
 pure module function thermal_T(ph,en) result(T)
 
@@ -332,7 +325,7 @@ end function thermal_T
 
 
 !----------------------------------------------------------------------------------------------
-!< @brief Get rate of temperature (for use by non-thermal physics)
+!< @brief Get rate of temperature (for use by non-thermal physics).
 !----------------------------------------------------------------------------------------------
 module function thermal_dot_T(ph,en) result(dot_T)
 

src/phase_thermal_source_dissipation.f90

@@ -5,7 +5,7 @@
 !> @brief material subroutine for thermal source due to plastic dissipation
 !> @details to be done
 !--------------------------------------------------------------------------------------------------
-submodule(phase:thermal) dissipation
+submodule(phase:thermal) source_dissipation
 
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pREAL) :: &
@@ -22,10 +22,10 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function dissipation_init(source_length) result(mySources)
+module function source_dissipation_init(maxNsources) result(isMySource)
 
-  integer, intent(in)                  :: source_length
+  integer, intent(in)                  :: maxNsources
-  logical, dimension(:,:), allocatable :: mySources
+  logical, dimension(:,:), allocatable :: isMySource
 
   type(tDict), pointer :: &
     phases, &
@@ -35,26 +35,29 @@ module function dissipation_init(source_length) result(mySources)
   class(tList), pointer :: &
     sources
   character(len=:), allocatable :: refs
-  integer :: so,Nmembers,ph
+  integer :: ph,Nmembers,so,Nsources
 
 
-  mySources = thermal_active('dissipation',source_length)
+  isMySource = thermal_active('dissipation',maxNsources)
-  if (count(mySources) == 0) return
+  if (count(isMySource) == 0) return
 
-  print'(/,1x,a)', '<<<+-  phase:thermal:dissipation init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal:source_dissipation init  -+>>>'
-  print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
   allocate(param(phases%length))
 
   do ph = 1, phases%length
+    Nsources = count(isMySource(:,ph))
+    if (Nsources == 0) cycle
+    if (Nsources > 1) call IO_error(600,ext_msg='dissipation')
+    Nmembers = count(material_ID_phase == ph)
     phase => phases%get_dict(ph)
-    if (count(mySources(:,ph)) == 0) cycle !ToDo: error if > 1
     thermal => phase%get_dict('thermal')
     sources => thermal%get_list('source')
     do so = 1, sources%length
-      if (mySources(so,ph)) then
+      if (isMySource(so,ph)) then
         associate(prm  => param(ph))
           src => sources%get_dict(so)
           print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
@@ -62,22 +65,21 @@ module function dissipation_init(source_length) result(mySources)
           if (len(refs) > 0) print'(/,1x,a)', refs
 
           prm%kappa = src%get_asReal('kappa')
-          Nmembers = count(material_ID_phase == ph)
           call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
-
         end associate
+        exit
       end if
     end do
   end do
 
 
-end function dissipation_init
+end function source_dissipation_init
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Ninstancess dissipation rate
 !--------------------------------------------------------------------------------------------------
-module function dissipation_f_T(ph,en) result(f_T)
+module function source_dissipation_f_T(ph,en) result(f_T)
 
   integer, intent(in) :: ph, en
   real(pREAL) :: &
@@ -91,6 +93,6 @@ module function dissipation_f_T(ph,en) result(f_T)
     f_T = prm%kappa*sum(abs(Mp*mechanical_L_p(ph,en)))
   end associate
 
-end function dissipation_f_T
+end function source_dissipation_f_T
 
-end submodule dissipation
+end submodule source_dissipation

src/phase_thermal_source_externalheat.f90

@@ -4,14 +4,14 @@
 !> @author Philip Eisenlohr, Michigan State University
 !> @brief material subroutine for variable heat source
 !--------------------------------------------------------------------------------------------------
-submodule(phase:thermal) externalheat
+submodule(phase:thermal) source_externalheat
 
 
   integer,           dimension(:),   allocatable :: &
-    source_thermal_externalheat_offset                                                              !< which source is my current thermal dissipation mechanism?
+    source_ID                                                                                       !< index in phase source list corresponding to this source
 
   type :: tParameters                                                                               !< container type for internal constitutive parameters
-    type(tTable) :: f
+    type(tTable) :: f                                                                               !< external heat power as (tabulated) function of time
   end type tParameters
 
   type(tParameters), dimension(:), allocatable  :: param                                            !< containers of constitutive parameters (len Ninstances)
@@ -24,10 +24,10 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-module function externalheat_init(source_length) result(mySources)
+module function source_externalheat_init(maxNsources) result(isMySource)
 
-  integer, intent(in)                  :: source_length
+  integer, intent(in)                  :: maxNsources
-  logical, dimension(:,:), allocatable :: mySources
+  logical, dimension(:,:), allocatable :: isMySource
 
   type(tDict), pointer :: &
     phases, &
@@ -37,28 +37,31 @@ module function externalheat_init(source_length) result(mySources)
   type(tList), pointer :: &
     sources
   character(len=:), allocatable :: refs
-  integer :: so,Nmembers,ph
+  integer :: ph,Nmembers,so,Nsources
 
 
-  mySources = thermal_active('externalheat',source_length)
+  isMySource = thermal_active('externalheat',maxNsources)
-  if (count(mySources) == 0) return
+  if (count(isMySource) == 0) return
 
-  print'(/,1x,a)', '<<<+-  phase:thermal:externalheat init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal:source_externalheat init  -+>>>'
-  print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(isMySource); flush(IO_STDOUT)
 
 
   phases => config_material%get_dict('phase')
   allocate(param(phases%length))
-  allocate(source_thermal_externalheat_offset  (phases%length), source=0)
+  allocate(source_ID(phases%length), source=0)
 
   do ph = 1, phases%length
+    Nsources = count(isMySource(:,ph))
+    if (Nsources == 0) cycle
+    if (Nsources > 1) call IO_error(600,ext_msg='externalheat')
+    Nmembers = count(material_ID_phase == ph)
     phase => phases%get_dict(ph)
-    if (count(mySources(:,ph)) == 0) cycle
     thermal => phase%get_dict('thermal')
     sources => thermal%get_list('source')
     do so = 1, sources%length
-      if (mySources(so,ph)) then
+      if (isMySource(so,ph)) then
-        source_thermal_externalheat_offset(ph) = so
+        source_ID(ph) = so
         associate(prm  => param(ph))
           src => sources%get_dict(so)
           print'(1x,a,i0,a,i0)', 'phase ',ph,' source ',so
@@ -66,41 +69,36 @@ module function externalheat_init(source_length) result(mySources)
           if (len(refs) > 0) print'(/,1x,a)', refs
 
           prm%f = table(src,'t','f')
-
-          Nmembers = count(material_ID_phase == ph)
           call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
         end associate
+        exit
       end if
     end do
   end do
 
-end function externalheat_init
+end function source_externalheat_init
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief rate of change of state
 !> @details state only contains current time to linearly interpolate given heat powers
 !--------------------------------------------------------------------------------------------------
-module subroutine externalheat_dotState(ph, en)
+module subroutine source_externalheat_dotState(ph, en)
 
   integer, intent(in) :: &
     ph, &
     en
 
-  integer :: &
-    so
 
-  so = source_thermal_externalheat_offset(ph)
+  thermalState(ph)%p(source_ID(ph))%dotState(1,en) = 1.0_pREAL                                         ! state is current time
 
-  thermalState(ph)%p(so)%dotState(1,en) = 1.0_pREAL                                                 ! state is current time
+end subroutine source_externalheat_dotState
-
-end subroutine externalheat_dotState
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local heat generation rate
 !--------------------------------------------------------------------------------------------------
-module function externalheat_f_T(ph,en) result(f_T)
+module function source_externalheat_f_T(ph,en) result(f_T)
 
   integer, intent(in) :: &
     ph, &
@@ -108,16 +106,11 @@ module function externalheat_f_T(ph,en) result(f_T)
   real(pREAL) :: &
     f_T
 
-  integer :: &
-    so
-
-
-  so = source_thermal_externalheat_offset(ph)
 
   associate(prm => param(ph))
-    f_T = prm%f%at(thermalState(ph)%p(so)%state(1,en))
+    f_T = prm%f%at(thermalState(ph)%p(source_ID(ph))%state(1,en))
   end associate
 
-end function externalheat_f_T
+end function source_externalheat_f_T
 
-end submodule externalheat
+end submodule source_externalheat