Dispersion Model (patch based)

patches-own [reactant-in product-in reactant product rxn-rate]
globals [centerline reactant-profile reactant-profile-turtle]


;;GLOBAL VARIABLES:
  ;;centerline = agentset of patches with pycor = 0
  ;;reactant-profile = ordered list of reactant concentrations along the longitudinal axis (centerline) of the current reactor
  ;;reactant-profile-turtle = ordered list of reactant concentrations along the longitudinal axis (centerline) of the imported turtle-based reactor


;;PATCH VARIABLES:
  ;;reactant,product = local concentration of reactant/product within the present patch
  ;;reactant-in,product-in = local concentration of reactant/product in the patch directly upstream (will become my concentration after convection)
  ;;rxn-rate = local change in reactant/product concentration, expressed in mol/L/s

to setup
  clear-all
  
  if fill-reactor?
    [ask patches [set reactant feed-concentration] 
     recolor]
  
  reset-ticks  
end 

to go
;;1 tick = 1 second

  ask patches with [pxcor = min-pxcor] [feed-reactor]
  convective-transport
  diffusive-transport
  ask patches with [pxcor > min-pxcor] [run-kinetics]
  calculate-performance
  recolor
  update-plots
  tick
end 

to feed-reactor  ;;Patch procedure
;;Ask patches to set their concentration to that of the feed.
  
    set reactant feed-concentration
    set product 0
end 

to run-kinetics  ;;Patch procedure
;;Calculate instantaneous deterministic rxn rate at each patch according to local concentration
;;Ask each patch to add the quantity reacted per tick to their previous concentration 

   set rxn-rate (rate-constant * ( reactant ^ rxn-order )) ;mol reacted / L / s
   set reactant ( reactant - rxn-rate )
   set product ( product + rxn-rate )
end 

to convective-transport
;;RUNS CONVECTION PROCESS FOR PATCHES  
  
;;Set ticks per convection step
  let steps-per-tick  (( feed-velocity * (world-width - 1)  ) / reactor-length)
  
  ;;Ask patches to undergo convection steps-per-tick times per tick
  let counter 0
  while [counter < steps-per-tick]
    [flow
      set counter (counter + 1)]
end 

to flow
;;ACTUAL CONVECTIVE TRANSPORT OCCURS HERE
  ;;Fluid moves one patch per tick.  
  ;;Patches store the incoming reactant and product before setting their own in order to make sure that every patch inherits the current value of its upstream neighbor.
  
  ;;Store concentration of incoming fluid
  ask patches with [pxcor > min-pxcor] [
    set reactant-in [reactant] of patch-at -1 0
    set product-in [product] of patch-at -1 0]
  
  ;;Set new concentrations    
  ask patches with [pxcor > min-pxcor] [
    set reactant reactant-in
    set product product-in]
end 

to diffusive-transport
;;DIFFUSIVE MECHANISM FOR THE PATCH-BASED MODEL
  
  ;;Ask patches to undergo diffusion "diffusivity" times per tick.  This keeps diffusion on the same time scale as convection, as they are both in units of "executions per tick".
  let counter 0
  while [counter < diffusivity]
    [diffuse reactant 1
      diffuse product 1
      set counter (counter + 1)]
end 

to calculate-performance
;;CALCULATES PERFORMANCE METRICS FOR THE REACTOR
  
  ;;Defines "centerline" as patches lying along the longitudinal axis.  Convenient for calculating profiles
  set centerline patches with [pycor = 0]
  
  ;;Store reactant concentration profile
  set reactant-profile (map [[reactant] of ?] (sort-on [pxcor] centerline))   
end 

to recolor
;;RECOLORS PATCHES AT EACH TICK. COLOR DEPENDS UPON VISUALIZATION MODE
  
  if visualization-mode = "reactant-profile"
    [ask patches [set pcolor scale-color blue reactant 0 (feed-concentration * 2)]] ;;Scales color between blue (high reactant concentration) and black (no/low reactant)
  
  if visualization-mode = "product-profile"
    [ask patches [set pcolor scale-color red product 0 (feed-concentration * 2)]] ;;Scales color between red (high product concentration) and black (no/low products)
  
  if visualization-mode = "blended-profile"
    [ask patches 
      [ifelse (reactant + product) = 0
        [set pcolor black ]                                                     ;;If patch is empty, set color black
        [let reactant-frac (reactant / (reactant + product))                                     ;;Reactant mole fraction 
          set pcolor rgb ((1 - reactant-frac) * 255) 0 (reactant-frac * 255)]]] ;;Scales color between blue (reactants) and red (products)
end 

to import-turtle-model
;;FUNCTION IMPORTS PARAMETERS AND RESULTS FROM TURTLE-BASED MODEL
;;THEN SETS CURRENT MODEL TO THOSE OPERATING PARAMETERS 
  
  
  ;;Open file containing data from turtle-based model
  file-open "turtle-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 "turtle-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 turtle model
  set feed-velocity       read-from-string (file-read-line)                 ;;Sets feed-velocity to turtle-based model's value
  set diffusivity         (2 * round (read-from-string (file-read-line)))   ;;Sets diffusivity to turtle-based model's value, rounded because patch model takes integer values.  The 8x is an approximate scaling factor between the two mechanisms..
  set rate-constant       read-from-string (file-read-line)                 ;;Sets rate constant to turtle-based model's reaction-probability 
  let reaction-type       file-read-line                                    ;;Stores the reaction type (i.e. conc. ind./dep. or bimolecular) as string
  set reactant-profile-turtle read-from-string (file-read-line)                 ;;Import results (conc. profile) from turtle-based model
  
  
                                                                                ;;Set appropriate rxn order based on turtle model's mechanism
  import-rxn-order reaction-type
  
  ;;Close file, resetting line count.
  file-close-all
end 

to import-rxn-order [reaction-type]
;;CONVERTS TURTLE-BASED MODEL MECHANISM TO RXN ORDER.  BASED ON INSIGHT FROM THIS MODEL!
  
  show reaction-type
  
  ifelse reaction-type = "Bimolecular" or reaction-type = "Concentration Independent"
    [set rxn-order 1]
    [set rxn-order 2]
end 

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

;;Normalize reactant concentration profile 
  let raw-profile (map [([reactant] of ?)] (sort-on [pxcor] centerline))
  let normalized-profile (map [? / (item (world-width - 10) raw-profile)] raw-profile)
  
  
  ;;Create file (or delete old one and reopen)   
  ifelse file-exists? "patch-results.csv"
    [file-delete "patch-results.csv"
      file-open "patch-results.csv" ]
    [file-open "patch-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 (rate-constant)  
  file-print normalized-profile     
  
  ;;Save and close file  
  file-flush
  file-close
end 

to-report conversion
;;CALCULATES FRACTIONAL CONVERSION
  
  ;;Conversion is defined as: ((reactant-in - reactant-out) / reactant-in )
  let outlet (mean [reactant] of patches with [pxcor = max-pxcor])
  let inlet  (mean [reactant] of patches with [pxcor = min-pxcor])
  ifelse inlet > 0
    [report (( inlet - outlet ) / inlet)]
    [report 0]
end 

to-report mean-sq-error [profile1 profile2]
;;REPORTS MEAN SQUARED ERROR BETWEEN PATCH AND TURTLE BASED MODELS

  ;;Turtle concentration profile is scaled to a feed of 1 before export.  Here it's rescaled based on the concentration at the point 10 patches from the end of the reactor.
  let scaled-patch-conc ([reactant] of (item (world-width - 10) (sort-on [pxcor] centerline)))
  let scaled-profile2 map [? * scaled-patch-conc ] 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