Dispersion Model (turtle based)

breed [inerts inert]
breed [reactants reactant]
breed [products product]
turtles-own [conv-xcor conv-ycor diff-xcor diff-ycor res-time]
globals [effluent effluent-reactant effluent-product effluent-rtd-list effluent-rtd reactant-profile reactant-profile-patches]


;;GLOBAL VARIABLES:
  ;;effluent = agentset of turtles leaving the right system boundary at each time point
  ;;effluent-reactant = effluent with reactant breed
  ;;effluent-product = effluent with product breed
  ;;effluent-rtd-list = list of lists containing the residence times of all agents in the effluent. At each tick, a new list for the current effluent is appended to the end, and the first list is removed if the sampling window is exceeded.
  ;;effluent-rtd = single list of residence times of all effluent within the sampling window
  ;;reactant-profile = list of total reactant count in the current model for each value of pxcor, ordered by pxcor
  ;;reactant-profile-patches = list of total reactant count imported from the patch-based model for each value of pxcor, ordered by pxcor


;;TURTLE VARIABLES:
  ;;conv-xcor/conv-ycor = the projected x-y coordinates of each turtle after convection
  ;;diff-xcor/diff-ycor = the projected x-y coordinates of each turtle after convection and dispersion
  ;;res-time = the residence time of an individual turtle, incremented at each tick

to setup
  clear-all
    
  ;;Set default shapes
  set-default-shape turtles "circle"
  
  ;;Initialize sets and lists
  set reactant-profile ([pycor] of patches with [pycor = 0])  ;;Create an array of zeroes with length equal to reactor length
  set effluent (turtles with [xcor = max-pxcor])
  set effluent-reactant 0
  set effluent-product 0
  set effluent-rtd-list []
  set effluent-rtd []
  
  ;;Fill reactor 
  ifelse fill-reactor? 
    [fill-reactor]
    [ask patches with [pxcor = min-pxcor] [feed-reactor]]
  
  reset-ticks  
end 

to fill-reactor
;;Fills entire reactor with reactant at setup.  This avoids time delay due to convection, but still incurs delay due to residence time.

  ;;Ask patches with xcor equal to multiple of flow discretization to sprout reactant
  ask patches with [pxcor = feed-velocity * floor (pxcor / feed-velocity)]
    [sprout-reactants feed-concentration
      [set color blue
        set size 1
        set heading 90
        set res-time 0]]
end 

to go
;;Note that for units, 1 tick = 1 second
  
  run-kinetics
  convective-transport
  dispersive-transport
  analyze-effluent
  update-plots
  move
  ask patches with [pxcor = min-pxcor] [feed-reactor]
  calculate-profile
  tick
end 

to feed-reactor  ;;Patch procedure
;;ADDS MOLECULAR PACKETS AT REACTOR INLET  
  
  ;;Add reactant to feed
  sprout-reactants feed-concentration
  [set color blue
    set size 1
    set heading 90
    set res-time 0]
  
  
  ;;Add tracer during step experiments
  if tracer-step? 
    [sprout-inerts tracer-density
      [set color yellow
        set size 2
        set heading 90
        set res-time 0]]
end 

to tracer-pulse
  ;;Add single pulse of tracer to feed
  
  ask patches with [pxcor = min-pxcor]
    [sprout-inerts tracer-density
      [set color yellow
        set size 2
        set heading 90
        set res-time 0]]
end 

to run-kinetics
;;RUNS REACTION KINETICS ACCORDING TO SELECTED MECHANISM


  ;;Conc. Independent: Reaction proceeds with equal probability regardless of local concentration
  if reaction-type = "Concentration Independent"
    [ask reactants 
      [if random-float 100 < reaction-probability
        [set breed products
         set color red]]]


  ;;Conc. Dependent:  Reaction probability scales proportionally with the local reactant concentration.
  ;;Reactant concentration is based upon reactants within the collision radius.
   if reaction-type = "Concentration Dependent"
    [ask reactants 
      [if random-float 100 < (reaction-probability * ( count (other reactants in-radius collision-radius) ))
        [set breed products
         set color red]]]

 
  ;;Bimolecular.  Reaction does not proceed unless another reactant is within the collision radius.
   if reaction-type = "Bimolecular"
      ;;Ask reactants to do: I check if I am within the collision-radius of another reactant, if so, 
      ;;and we probabilistically have the correct orientation, we react.  
      ;;When we react, I ask the other reactant to die, then I become the product.
 
    [ask reactants 
      [if (count (other reactants in-radius collision-radius) > 0) and random-float 100 < reaction-probability  
        [ask one-of other reactants in-radius collision-radius 
          [die]
         set breed products
         set color red]]]
end 

to convective-transport
;;BULK FLOW TRANSPORT MECHANISM
    
  ;;Ask turtles to predict their downstream position
  ask turtles [
    set heading 90
    set conv-xcor (xcor + feed-velocity) ]
end 

to dispersive-transport
;;DISPERSIVE TRANSPORT MECHANISM
;;Ask turtles to predict their post-diffusion position relative to their predicted post-convection position.
  
  ;;Tell all turtles: I take a random walk of path-length equal to "diffusivity." 
  ask turtles [ 
    
    ;;If random diffusion is on, the random walk step size is a random float less than the diffusivity.  This helps smooth the data, but isn't a true fixed-step random walk.
    ifelse not random-diffusion?
      [ifelse random 2 = 1
        [set diff-xcor (conv-xcor + (diffusivity))]
        [set diff-xcor (conv-xcor - (diffusivity))]]
    
      [ifelse random 2 = 1
        [set diff-xcor (conv-xcor + (random-float diffusivity))]
        [set diff-xcor (conv-xcor - (random-float diffusivity))]]
    
    
    if diff-xcor < min-pxcor
      [set diff-xcor min-pxcor]]
end 

to analyze-effluent
;;BOOKKEEPING OF ELEMENTS LEAVING THE REACTOR
  
  
  ;;IDENTIFY PARTICLES LEAVING THE REACTOR
  
  ;;Store all molecules leaving the reactor as effluent
  set effluent (turtles with [diff-xcor > max-pxcor])
 
  ;;Store effluent reactant (moving average does not store past values to save memory)
  set effluent-reactant (((sampling-window - 1) / sampling-window) * effluent-reactant) + ((1 / sampling-window) * (count effluent with [breed = reactants]) )

  ;;Store effluent product (moving average does not store past values to save memory)
  set effluent-product (((sampling-window - 1) / sampling-window) * effluent-product) + ((1 / sampling-window) * (count effluent with [breed = products]) )
   
   
  ;;STORE EFFLUENT RESIDENCE TIMES
  
  ;;Add effluent residence times to sliding window. effluent-rtd-list is a list of sets of residence times for the effluent molecules at each tick.
  set effluent-rtd-list lput ([res-time] of effluent) effluent-rtd-list
  
  ;;Removes the oldest set of residence times.
  if (length effluent-rtd-list) > sampling-window
    [set effluent-rtd-list remove-item 0 effluent-rtd-list]
    
  ;;Combines residence time lists into a single list for plotting.
  set effluent-rtd effluent-rtd-merged
end 

to move
;;EXECUTES ACTUAL MOTION

  ;;Remove effluent from reactor
  ask effluent [die]
  
  ;;Turtles reach their projected position after convective + diffusive transport, then increment residence time
  ask turtles [
    set xcor diff-xcor
    set res-time (res-time + 1)]
end 

to calculate-profile
  
  ;;Ask each patch along the longitudinal axis of the reactor to store the current total number of reactants with the same x-coordinate
  ;;This value is then appended to the moving average reactant-profile across the sampling window

  let centerline (patches with [pycor = 0])
  let inst-reactant-profile (map ([sum ([count reactants-here] of patches with [pxcor = [pxcor] of ?])]) (sort-on [pxcor] centerline))
  set reactant-profile (  map [((1 / sampling-window) * ?1) + (((sampling-window - 1) / sampling-window) * ?2)] inst-reactant-profile reactant-profile)
end 




;;End of main model code.  Remaining code is for analysis and import/export.
;;===================================================================================


;;IMPORT EXPORT CODE FOR SYNCING WITH PATCH-BASED MODEL

to import-patch-model
;;FUNCTION IMPORTS PARAMETERS AND RESULTS FROM PATCH-BASED MODEL
;;THEN SETS CURRENT MODEL TO THOSE OPERATING PARAMETERS 
  
  ;;Open file containing data from turtle-based model
  file-open "patch-results.csv"
  
  ;;If file is at the end (i.e. already imported once), close it and reopen it.
  if file-at-end? 
    [file-close-all
      file-open "patch-results.csv"]
  
  ;;Store each of the variables from the other model.  List order is known a priori, not ideal but... no tengo tiempo!
  let dim                 read-from-string (file-read-line)
  resize-world (item 0 dim) (item 1 dim) (item 2 dim) (item 3 dim)   ;;Set dimensions to those of patch model
  set feed-velocity       read-from-string (file-read-line)          ;;Sets feed-velocity to patch-based model's value
  set diffusivity         round (read-from-string (file-read-line))  ;;Sets diffusivity to patch-based model's value.
  set reaction-probability (100 * read-from-string (file-read-line)) ;;Sets reaction-probability to patch-based model's rate-constant
  set reactant-profile-patches read-from-string (file-read-line)          ;;Import results (conc. profile) from patch-based model
  
                                                                          ;;Close file, resetting line count.
  file-close-all
end 

to export-results
;;WRITE OPERATING PARAMETERS AND CONCENTRATION PROFILE TO A FILE FROM WHICH IT CAN BE IMPORTED TO PATCH MODEL

  ;;Normalize reactant concentration profile to stable concentration near end of reactor. I've shifted this right because the bounds on random walk cause an anomaly near the inlet.
  let reactant-profile-norm (map [? / (item (world-width - 10) reactant-profile)] reactant-profile )
  
  ;;Create file (or delete old one and reopen)   
  ifelse file-exists? "turtle-results.csv"
    [file-delete "turtle-results.csv"
      file-open "turtle-results.csv" ]
    [file-open "turtle-results.csv"]
  
  ;;Write variables to file.  I've left out labels because I can't figure out how to make the import step recognize them.
  file-print (list min-pxcor max-pxcor min-pycor max-pycor)
  file-print feed-velocity
  file-print diffusivity
  file-print (reaction-probability / 100)  ;;Converted to percentage probability
  file-print reaction-type
  file-print reactant-profile-norm     
  
  ;;Save and close file  
  file-flush
  file-close
end 




;; REPORTERS - ONLY USED FOR ANALYSIS

to-report space-time
;;CALCULATE SPACE TIME (convection/reactor volume)
 
  ifelse feed-velocity > 0
    [report (((2 * max-pxcor) + 1) / feed-velocity)]
    [report "inf"]
end 

to-report conversion
;;CALCULATE CONVERSION  

  let reactant-in (count reactants with [res-time = 0])
  report (( reactant-in - effluent-reactant ) / reactant-in)
end 

to-report peclet []
;;CALCULATE PECLET NUMBER (dimensionless ratio = convection/dispersion)
 
  if length effluent-rtd > 1 
    [ifelse non-dim-res-time?
      [if length effluent-rtd > 1
        [ifelse variance effluent-rtd > 0 
          [report (2 / (variance effluent-rtd))]
          [report "inf."] ]]
      [ifelse variance effluent-rtd != 0
        [report ((2 *  (space-time ^ 2)) / (variance effluent-rtd)) ]
        [report "inf."]]]
end 

to-report effluent-rtd-merged
  
;;CALCULATE RESIDENCE TIME DISTRIBUTION   
 
  ;;Takes list of lists and combines them into a single list of residence times for the given sliding window
  let rtd-list reduce [sentence ?1 ?2] effluent-rtd-list
  
  ;;Nondimensionalize residence times
  ifelse non-dim-res-time? 
    [report map [ ? / (space-time)] rtd-list]
    [report rtd-list]
end 

to-report mean-sq-error [profile1 profile2]
  
  ;;Rescales patch profile based on concentration ten patches from the reactor exit
  let scaled-profile2 (map [? *  (item (world-width - 10) profile1)] profile2)
  
  ;;Filters out first ten and last five entries from profile.  These are the regions where noise has a strong influence.
  let filtered-profile1 sublist profile1 (10) ((length profile1) - 5)
  let filtered-profile2 sublist scaled-profile2 (10) ((length scaled-profile2) - 5)
  
  ;;Calculates sum of mean squared error between profiles
  let squared-errors (map [(?1 - ?2) ^ 2] filtered-profile1 filtered-profile2)

  ;;Report sum of squared error
  report (sum squared-errors)   
end