oregonator.xmds

<?xml version="1.0"?>
<simulation>
  
<!-- $Id: oregonator_body.part 996 2004-08-03 05:32:16Z cochrane $ -->

<!--  Copyright (C) 2000-2004                                           -->
<!--                                                                    -->
<!--  Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane -->
<!--                                                                    -->
<!--  This file is part of xmds.                                        -->
<!--                                                                    -->
<!--  This program is free software; you can redistribute it and/or     -->
<!--  modify it under the terms of the GNU General Public License       -->
<!--  as published by the Free Software Foundation; either version 2    -->
<!--  of the License, or (at your option) any later version.            -->
<!--                                                                    -->
<!--  This program is distributed in the hope that it will be useful,   -->
<!--  but WITHOUT ANY WARRANTY; without even the implied warranty of    -->
<!--  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     -->
<!--  GNU General Public License for more details.                      -->
<!--                                                                    -->
<!--  You should have received a copy of the GNU General Public License -->
<!--  along with this program; if not, write to the Free Software       -->
<!--  Foundation, Inc., 59 Temple Place - Suite 330, Boston,            -->
<!--  MA  02111-1307, USA.                                              -->

  <name> oregonator </name>      <!-- the name of the simulation -->
  
  <author> Paul Cochrane </author>  <!-- the author of the simulation -->
  <description>
    <!-- a description of what the simulation is supposed to do -->
    Example simulation of the Oregonator model of oscillating chemical
    kinetics equations.
    Calculates concentrations of components participating 
    in the Oregonator model.  The reaction scheme is

		A + Y ----> X + R        (1)
		X + Y ----> 2 R          (2)
		A + X ----> 2 X + 2 Z    (3)
		2 X ----> A + R          (4)
		B + Z ----> 1/2 Y        (5)

    Rate laws are

	d[X]/dt = k1 [A] [Y] - k2 [X] [Y] + k3 [A] [X] - 2 k4 [X]^2
	d[Y]/dt = -k1 [A] [Y] - k2 [X] [Y] + (1/2) k5 [B] [Z]
	d[Z]/dt = 2 k3 [A] [X] - k5 [B] [Z]

    which are transformed to

	d[X]_d(tau) = (k1 [A] [Y] - k2 [X] [Y] + k3 [A] [X] - 2 k4 [X]^2)/(k5 [B])
	d[Y]_d(tau) = (-k1 [A] [Y] - k2 [X] [Y] + (1/2) k5 [B] [Z])/(k5 [B])
	d[Z]_d(tau) = (2 k3 [A] [X] - k5 [B] [Z])/(k5 [B])

    where tau = k5*[B]*t is a unitless time variable.  The
    program integrates these equations, assuming that concentrations
    of A and B are constant.

    Adapted for xmds from "Mathematica computer programs for physical
    chemistry", William H. Cropper, Springer Verlag (1998)
  </description>
  
  <!-- Global system parameters and functionality -->
  <prop_dim> tau </prop_dim>    <!-- name of main propagation dim -->
  
  <error_check> yes </error_check>   <!-- defaults to yes -->
  <use_wisdom> yes </use_wisdom>     <!-- defaults to no -->
  <benchmark> yes </benchmark>       <!-- defaults to no -->
  <use_prefs> yes </use_prefs>       <!-- defaults to yes -->
  
  <!-- Global variables for the simulation -->
  <globals>
  <![CDATA[
    // rate constants
    //   Units are  L mol^-1 s^-1
    const double k1 = 0.005;
    const double k2 = 1.0;
    const double k3 = 1.0;
    const double k4 = 1.0;
    const double k5 = 0.0004;
    
    // the constant concentrations of the components A and B
    // in mol L^-1
    const double a = 0.01;
    const double b = 1.0;

    // the initial concentrations of the components X, Y and Z
    const double CX0 = 0.0;
    const double CY0 = 0.0;
    const double CZ0 = 0.00025;
  ]]>
  </globals>
  
  <!-- Field to be integrated over -->
  <field>
    <name> main </name>
    <samples> 1 </samples>       <!-- sample 1st point of dim? -->
    
    <vector>
      <name> main </name>
      <type> double </type>           <!-- data type of vector -->
      <components> CX CY CZ </components>       <!-- names of components -->
      <![CDATA[
        CX = CX0;
	CY = CY0;
	CZ = CZ0;
      ]]>
    </vector>
  </field>
  
  <!-- The sequence of integrations to perform -->
  <sequence>
    <integrate>
      <algorithm> RK4IP </algorithm> <!-- RK4EX, RK4IP, SIEX, SIIP -->
      <interval> 50 </interval>   <!-- how far in main dim? -->
      <lattice> 5000 </lattice>     <!-- no. points in main dim -->
      <samples> 500 </samples> <!-- no. pts in output moment group -->
      
      <![CDATA[
        dCX_dtau = (k1*a*CY - k2*CX*CY + k3*a*CX - 2.0*k4*CX*CX)/(k5*b);
	dCY_dtau = (-k1*a*CY - k2*CX*CY + 0.5*k5*b*CZ)/(k5*b);
	dCZ_dtau = (2.0*k3*a*CX - k5*b*CZ)/(k5*b);
      ]]>
    </integrate>
  </sequence>
  
  <!-- The output to generate -->
  <output format="ascii">
    <group>
      <sampling>
        <lattice> 500 </lattice>           <!-- no. points to sample -->
        <moments> X X10 Y Z </moments>           <!-- names of moments -->
        <![CDATA[
          X = CX;
	  X10 = CX*10;  // scaled X concentration for plotting
	  Y = CY;
	  Z = CZ;
        ]]>
      </sampling>
    </group>
  </output>
  
</simulation>