Merge lp:~jamesmadd/dolfin-adjoint/timestepping into lp:dolfin-adjoint

[James Maddison]

Timestepping library, with manual, tests, and dolfin-adjoint integration

[lp:~jamesmadd/dolfin-adjoint/timestepping]

644. By James Maddison <email address hidden> on 2013-05-08

Corrections to the manual

645. By James Maddison <email address hidden> on 2013-05-09

Do not delete the checkpoints directory

646. By James Maddison <email address hidden> on 2013-05-09

Create the checkpoints directory

647. By James Maddison <email address hidden> on 2013-05-09

Do not delete checkpoint files

648. By James Maddison <email address hidden> on 2013-05-09

Fix a parallel deadlock when checkpointing to disk

[Patrick Farrell (pefarrell)]

Hi James!

Sorry for the delay.

Very minor comments. All of these are optional suggestions, feel free to disregard.

a) Can you make your test script edit the PYTHONPATH/sys.path to always use the local copy of timestepping?

b) Is there any way to run the tests in parallel?

c) Actually, more generally, is there a compelling reason not to use the dolfin-adjoint test tool?

d) Is it safe to always put checkpoints in checkpoints~/? That seems a bit strange. Is it safe if more than one process is executing the same code on the same machine at once?

e) Can you make the makefile in the manual delete the .dvi and the .ps? (Is there a reason you don't use pdflatex?)

f) In the manual, can you hypersetup those hideous red boxes away? Maybe:

\definecolor{DarkBlue}{rgb}{0.00,0.00,0.55}
\hypersetup{
    linkcolor = DarkBlue,
    anchorcolor = DarkBlue,
    citecolor = DarkBlue,
    filecolor = DarkBlue,
    pagecolor = DarkBlue,
    urlcolor = DarkBlue,
    colorlinks = true,
}

g) In the manual, you have a "is non-linear then the a non-linear equation .. must be solved." Get rid of the the.

h) In the manual, maybe use the same lstlisting settings as we have in the paper?

i) When I run the tests with the test harness (precise, fenics from PPA), I get

Successes: 51
Failures: 2
Failures:
/data/pfarrell/src/dolfin-adjoint/timestepping/tests/short/burgers_leapfrog_da
/data/pfarrell/src/dolfin-adjoint/timestepping/tests/short/linearised_shallow_water_cn/shallow_water

burgers_leapfrog_da fails with

-4.89248513609469882e+01 2.22399876292911358e-12
Traceback (most recent call last):
  File "./burgers_leapfrog_da", line 114, in <module>
    assert(err < 7.0e-13)
AssertionError

shallow_water fails with

Assembling form with rank 2
4.10821370684035065e-05
4.10821370682097189e-05
1.93787585804627849e-16
8.51321880127267620e-07
8.51321895311346568e-07
1.51840789485756049e-14
Traceback (most recent call last):
  File "./shallow_water", line 73, in <module>
    assert(err < 7.0e-15)
AssertionError

j) How do you feel about refactoring timestepping.py into multiple files .. ?

[James Maddison (jamesmadd)]

a, e, f, g, h, i, j are all straightforward, so I'll fix these.

Regarding b, c, d:

> d) Is it safe to always put checkpoints in checkpoints~/? That seems a bit
> strange. Is it safe if more than one process is executing the same code on the
> same machine at once?

I'm not actually sure how to handle this one. It's convenient to have the checkpoint files in a single directory, so that it's easy to clean them up afterwards, but it's inconvenient in that you can only run one test in a directory at once. Do you have any suggestions?

[lp:~jamesmadd/dolfin-adjoint/timestepping]

649. By James Maddison <email address hidden> on 2013-05-09

Always pick up local library versions when running the timestepping tests

650. By James Maddison <email address hidden> on 2013-05-09

Fix link style, lstlisting style, and a typo

651. By James Maddison <email address hidden> on 2013-05-09

Relax some test tolerances

652. By James Maddison <email address hidden> on 2013-05-09

Fix bibliography formatting

653. By James Maddison <email address hidden> on 2013-05-09

Start refactoring the timestepping library

654. By James Maddison <email address hidden> on 2013-05-09

Fix to test script

655. By James Maddison <email address hidden> on 2013-05-09

Separate out exception classes

656. By James Maddison <email address hidden> on 2013-05-09

Separate out C++ embedding

[Patrick Farrell (pefarrell)]

Maybe you could use the tempfile.mkdtemp function to create a guaranteed-unique temporary directory, and put them in that?

While looking at what dolfin-adjoint does, I see it doesn't do anything better than your timestepping library at the moment ..

[lp:~jamesmadd/dolfin-adjoint/timestepping]

657. By James Maddison <email address hidden> on 2013-05-10

Separate out version code and FEniCS patches

658. By James Maddison <email address hidden> on 2013-05-10

Separate out VTK I/O

659. By James Maddison <email address hidden> on 2013-05-10

Separate out FEniCS utilities, quadrature, and time levels

[James Maddison (jamesmadd)]

> Maybe you could use the tempfile.mkdtemp function to create a guaranteed-
> unique temporary directory, and put them in that?

That might be a problem if a model leaves large amounts of data in /tmp (e.g. due to a crash or termination).

[lp:~jamesmadd/dolfin-adjoint/timestepping]

660. By James Maddison <email address hidden> on 2013-05-10

Separate out time dependent functions

[James Maddison (jamesmadd)]

c) Actually, more generally, is there a compelling reason not to use the dolfin-adjoint test tool?

Does the tool support parallel tests with MPI?

[lp:~jamesmadd/dolfin-adjoint/timestepping]

661. By James Maddison <email address hidden> on 2013-05-10

Separate out static declarations

662. By James Maddison <email address hidden> on 2013-05-10

Separate out system solvers and caches

663. By James Maddison <email address hidden> on 2013-05-10

Separate out pre-assembly

664. By James Maddison <email address hidden> on 2013-05-10

Final refactoring

[James Maddison (jamesmadd)]

a, e, f, g, h, i, j are complete. python/timestepping/time_systems.py is still a large source file, but contains four (large) quite closely connected classes.

[lp:~jamesmadd/dolfin-adjoint/timestepping]

665. By James Maddison <email address hidden> on 2013-05-10

Fix circular dependency

666. By James Maddison <email address hidden> on 2013-05-12

- Give classes in time_systems.py more descriptive names
- More graceful resolution of a namespace clash with dolfin-adjoint
- Make checkpointing less likely to lead to deadlocks in parallel
- Make dependency processing less likely to lead to deadlocks in parallel
- Work around an obscure bug with FEniCS 1.0 and SciPy 0.12.0
- Add is_static method to pre-assembled forms
- Make static bi-linear form pre-assembly more repeatable

667. By James Maddison <email address hidden> on 2013-05-12

Delete manual test script

668. By James Maddison <email address hidden> on 2013-05-12

Delete redundant exclusion of 'test'

669. By James Maddison <email address hidden> on 2013-05-13

Add multi-process support to the test script

[James Maddison (jamesmadd)]

> b) Is there any way to run the tests in parallel?

This has also been addressed. Run the test script with -n [PROCS] (PROCS = 0 for all system cores).

[lp:~jamesmadd/dolfin-adjoint/timestepping]

670. By James Maddison <email address hidden> on 2013-05-13

Minor bugfix to parallel testing

671. By James Maddison <email address hidden> on 2013-05-13

Parallel testing optimisation

672. By James R. Maddison <email address hidden> on 2013-05-13

Add dates of first addition of FEniCS derived code

673. By James R. Maddison <email address hidden> on 2013-05-14

Documentation and constructor argument fix

674. By James R. Maddison <email address hidden> on 2013-05-15

Add some code addition dates

=== added directory 'timestepping'
=== added file 'timestepping/COPYING'
--- timestepping/COPYING	1970-01-01 00:00:00 +0000
+++ timestepping/COPYING	2013-05-15 14:57:23 +0000
@@ -0,0 +1,674 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    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 3 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, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    <program>  Copyright (C) <year>  <name of author>
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.

=== added file 'timestepping/COPYING.LESSER'
--- timestepping/COPYING.LESSER	1970-01-01 00:00:00 +0000
+++ timestepping/COPYING.LESSER	2013-05-15 14:57:23 +0000
@@ -0,0 +1,165 @@
+                   GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+  This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+  0. Additional Definitions.
+
+  As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+  "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+  An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+  A "Combined Work" is a work produced by combining or linking an
+Application with the Library.  The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+  The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+  The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+  1. Exception to Section 3 of the GNU GPL.
+
+  You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+  2. Conveying Modified Versions.
+
+  If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+   a) under this License, provided that you make a good faith effort to
+   ensure that, in the event an Application does not supply the
+   function or data, the facility still operates, and performs
+   whatever part of its purpose remains meaningful, or
+
+   b) under the GNU GPL, with none of the additional permissions of
+   this License applicable to that copy.
+
+  3. Object Code Incorporating Material from Library Header Files.
+
+  The object code form of an Application may incorporate material from
+a header file that is part of the Library.  You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+   a) Give prominent notice with each copy of the object code that the
+   Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the object code with a copy of the GNU GPL and this license
+   document.
+
+  4. Combined Works.
+
+  You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+   a) Give prominent notice with each copy of the Combined Work that
+   the Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the Combined Work with a copy of the GNU GPL and this license
+   document.
+
+   c) For a Combined Work that displays copyright notices during
+   execution, include the copyright notice for the Library among
+   these notices, as well as a reference directing the user to the
+   copies of the GNU GPL and this license document.
+
+   d) Do one of the following:
+
+       0) Convey the Minimal Corresponding Source under the terms of this
+       License, and the Corresponding Application Code in a form
+       suitable for, and under terms that permit, the user to
+       recombine or relink the Application with a modified version of
+       the Linked Version to produce a modified Combined Work, in the
+       manner specified by section 6 of the GNU GPL for conveying
+       Corresponding Source.
+
+       1) Use a suitable shared library mechanism for linking with the
+       Library.  A suitable mechanism is one that (a) uses at run time
+       a copy of the Library already present on the user's computer
+       system, and (b) will operate properly with a modified version
+       of the Library that is interface-compatible with the Linked
+       Version.
+
+   e) Provide Installation Information, but only if you would otherwise
+   be required to provide such information under section 6 of the
+   GNU GPL, and only to the extent that such information is
+   necessary to install and execute a modified version of the
+   Combined Work produced by recombining or relinking the
+   Application with a modified version of the Linked Version. (If
+   you use option 4d0, the Installation Information must accompany
+   the Minimal Corresponding Source and Corresponding Application
+   Code. If you use option 4d1, you must provide the Installation
+   Information in the manner specified by section 6 of the GNU GPL
+   for conveying Corresponding Source.)
+
+  5. Combined Libraries.
+
+  You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+   a) Accompany the combined library with a copy of the same work based
+   on the Library, uncombined with any other library facilities,
+   conveyed under the terms of this License.
+
+   b) Give prominent notice with the combined library that part of it
+   is a work based on the Library, and explaining where to find the
+   accompanying uncombined form of the same work.
+
+  6. Revised Versions of the GNU Lesser General Public License.
+
+  The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+  Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+  If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.

=== added directory 'timestepping/manual'
=== added file 'timestepping/manual/bibliography.bib'
--- timestepping/manual/bibliography.bib	1970-01-01 00:00:00 +0000
+++ timestepping/manual/bibliography.bib	2013-05-15 14:57:23 +0000
@@ -0,0 +1,51 @@
+% This file was created with JabRef 2.9b2.
+% Encoding: UTF-8
+
+@BOOK{donea2003,
+  title = {Finite element methods for flow problems},
+  publisher = {John Wiley \& Sons},
+  year = {2003},
+  author = {Donea, J. and Huerta, A.},
+  owner = {james},
+  quality = {1},
+  timestamp = {2013.05.08}
+}
+
+@BOOK{hairer1987,
+  title = {Solving ordinary differential equations I: Nonstiff problems},
+  publisher = {Springer},
+  year = {1987},
+  author = {Hairer, E. and N{\o}rsett, S. P. and Wanner, G.},
+  owner = {james},
+  timestamp = {2013.05.07}
+}
+
+@BOOK{iserles2009,
+  title = {A first course in the numerical analysis of differential equations},
+  publisher = {Cambridge University Press},
+  year = {2009},
+  author = {Iserles, A.},
+  edition = {2nd},
+  owner = {james},
+  timestamp = {2013.05.07}
+}
+
+@BOOK{mitchell1980,
+  title = {The finite difference method in partial differential equations},
+  publisher = {John Wiley \& Sons},
+  year = {1980},
+  author = {Mitchell, A. R. and Griffiths, D. F.},
+  owner = {james},
+  timestamp = {2013.05.07}
+}
+
+@BOOK{logg2012,
+  title = {Automated solution of differential equations by the finite element
+	method: The FEniCS book},
+  publisher = {Springer},
+  year = {2012},
+  editor = {Logg, A. and Mardal, K.-A. and Wells, G. N.},
+  owner = {james},
+  timestamp = {2013.01.09}
+}
+

=== added directory 'timestepping/manual/examples'
=== added file 'timestepping/manual/examples/TimeFunction'
--- timestepping/manual/examples/TimeFunction	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/TimeFunction	2013-05-15 14:57:23 +0000
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit square
+mesh = UnitSquareMesh(10, 10)
+# P1 function space
+space_p1 = FunctionSpace(mesh, "CG", 1)
+# P2_{DG} function space
+space_p2dg = FunctionSpace(mesh, "DG", 2)
+
+# Define two different sets of time levels
+levels_1 = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+levels_2 = TimeLevels(levels = [n], cycle_map = {},
+  last_past_level = n - 1)
+
+# Define time dependent functions on the levels
+F1 = TimeFunction(levels_1, space_p1, name = "F1")
+F2 = TimeFunction(levels_2, space_p2dg, name = "F2")
+
+# An equation involving the two time dependent functions
+test_p1 = TestFunction(space_p1)
+eq = inner(test_p1, F1[n + 1]) * dx == inner(test_p1, F2[n]) * dx
\ No newline at end of file

=== added file 'timestepping/manual/examples/TimeLevel'
--- timestepping/manual/examples/TimeLevel	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/TimeLevel	2013-05-15 14:57:23 +0000
@@ -0,0 +1,49 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from timestepping import *
+from fractions import Fraction
+
+# Create a TimeLevel with an offset of 2
+n1 = TimeLevel(2)
+# Create a TimeLevel with an offset of -1/2
+n2 = TimeLevel(-Fraction(1, 2))  
+# Create a FinalTimeLevel with an offset of 0
+N1 = FinalTimeLevel()            
+
+# Create a TimeLevel with an offset of 2
+n3 = n + 2
+# Create a TimeLevel with an offset of -3/2
+n4 = n - Fraction(3, 2)
+
+# Create a FinalTimeLevel with an offset of 1
+N2 = N + 1
+
+print n1 == n2  # False
+print n1 > n2   # True
+print n1 == n3  # True
+print n1 >= n3  # True
+
+print N1 == N2  # False
+print N1 <= N2  # True
+
+assert(not n1 == n2)
+assert(n1 > n2)
+assert(n1 == n3)
+assert(n1 >= n3)
+
+assert(not N1 == N2)
+assert(N1 <= N2)
\ No newline at end of file

=== added file 'timestepping/manual/examples/TimeLevels'
--- timestepping/manual/examples/TimeLevels	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/TimeLevels	2013-05-15 14:57:23 +0000
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from timestepping import *
+from fractions import *
+
+# Create a TimeLevels with one past level and one future level, with
+# the past level data replaced by the future level data during the
+# timestep variable cycle. Suitable for use with forward Euler,
+# backward Euler, or Crank-Nicolson schemes.
+levels_1 = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+
+# Create a TimeLevels with one future level
+levels_2 = TimeLevels(levels = [n], cycle_map = {},
+  last_past_level = n - 1)
+
+# Create a TimeLevels with one past level and two future levels,
+# with the past level data replaced by the latest future level data
+# during the timestep variable cycle. Suitable for use with a
+# second order Runge-Kutta scheme.
+levels_3 = TimeLevels(levels = [n, n + Fraction(1, 2), n + 1],
+  cycle_map = {n:n + 1})
+
+# Create a TimeLevels with three past levels and one future level,
+# with the past levels replaced by the nearest later level during
+# the timestep variable cycle. Suitable for use with a third order
+# Adams-Bashforth scheme.
+levels_4 = TimeLevels(levels = [n - 2, n - 1, n, n + 1],
+  cycle_map = {n - 2:n - 1, n - 1:n, n:n + 1})
\ No newline at end of file

=== added file 'timestepping/manual/examples/add_solve'
--- timestepping/manual/examples/add_solve	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/add_solve	2013-05-15 14:57:23 +0000
@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = Constant(0.05)
+bc1 = DirichletBC(space, 1.0, "on_boundary && near(x[0], 0.0)")
+bc2 = DirichletBC(space, 0.0, "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = Constant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = Function(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
\ No newline at end of file

=== added file 'timestepping/manual/examples/add_solve_statics'
--- timestepping/manual/examples/add_solve_statics	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/add_solve_statics	2013-05-15 14:57:23 +0000
@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = StaticConstant(0.05)
+bc1 = StaticDirichletBC(space, 1.0,
+  "on_boundary && near(x[0], 0.0)")
+bc2 = StaticDirichletBC(space, 0.0,
+  "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = StaticConstant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = StaticFunction(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
\ No newline at end of file

=== added file 'timestepping/manual/examples/diffusion_fe'
--- timestepping/manual/examples/diffusion_fe	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/diffusion_fe	2013-05-15 14:57:23 +0000
@@ -0,0 +1,65 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = StaticConstant(0.05)
+bc1 = StaticDirichletBC(space, 1.0,
+  "on_boundary && near(x[0], 0.0)")
+bc2 = StaticDirichletBC(space, 0.0,
+  "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = StaticConstant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = StaticFunction(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
+
+# Assemble the TimeSystem
+system = system.assemble()
+
+# Timestep the model
+t = 0.0
+while t * (1.0 + 1.0e-9) < 1.0:
+  system.timestep()
+  t += float(dt)
+# Finalise
+system.finalise()
\ No newline at end of file

=== added file 'timestepping/manual/examples/diffusion_fe_adjoined'
--- timestepping/manual/examples/diffusion_fe_adjoined	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/diffusion_fe_adjoined	2013-05-15 14:57:23 +0000
@@ -0,0 +1,77 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = StaticConstant(0.05)
+bc1 = StaticDirichletBC(space, 1.0,
+  "on_boundary && near(x[0], 0.0)")
+bc2 = StaticDirichletBC(space, 0.0,
+  "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = StaticConstant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = StaticFunction(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
+
+# Assemble the TimeSystem, enabling the adjoint. Set the
+# functional to be equal to spatial integral of the final u.
+system = system.assemble(adjoint = True, functional = u[N] * dx)
+
+# Timestep the model
+t = 0.0
+while t * (1.0 + 1.0e-9) < 1.0:
+  system.timestep()
+  t += float(dt)
+# Finalise
+system.finalise()
+
+# Perform a model constrained derivative calculation
+dJdm = system.compute_gradient(nu)
+
+# Verify the stored forward model data
+system.verify_checkpoints()
+# Verify the computed derivative using a Taylor remainder
+# convergence test
+orders = system.taylor_test(nu, grad = dJdm)
+# Check the convergence order
+assert((orders > 1.99).all())
\ No newline at end of file

=== added file 'timestepping/manual/examples/statics'
--- timestepping/manual/examples/statics	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/statics	2013-05-15 14:57:23 +0000
@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# Create a constant which is declared as static
+c = StaticConstant(1.0, name = "c")
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+# Create a function which is declared as static
+F = StaticFunction(space, name = "F")
+# Initialise the static function
+F.interpolate(Expression("x[0]"))
+
+# Create a Dirichlet boundary condition which is declared as static
+bc = StaticDirichletBC(space, 0.0, "on_boundary")
\ No newline at end of file

=== added file 'timestepping/manual/examples/time_discretisations'
--- timestepping/manual/examples/time_discretisations	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/time_discretisations	2013-05-15 14:57:23 +0000
@@ -0,0 +1,267 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+mesh = UnitIntervalMesh(10)
+space = FunctionSpace(mesh, "CG", 1)
+test = TestFunction(space)
+
+T_ic = StaticFunction(space, name = "T_ic")
+T_ic.assign(Constant(1.0))
+x = StaticFunction(space, name = "x")
+x.interpolate(Expression("x[0]"))
+def F(test, T):
+  return inner(test, x) * dx
+dt = StaticConstant(0.5)
+solver_parameters = {"linear_solver":"lu"}
+
+def check(name, system, T_f, ns_i = 0, tol = 0.0):
+  system = system.assemble()
+  system.timestep(ns = 2)
+  system.finalise()
+  if not isinstance(T_f, dolfin.Function):
+    lT_f = Function(space, name = "T_f")
+    pa_solve(inner(test, lT_f) * dx == inner(test, T_f) * dx,
+      lT_f, solver_parameters = solver_parameters)
+    T_f = lT_f
+  err = (T_f.vector() - (T_ic.vector() + 0.5 * (ns_i + 2) * x.vector())).norm("linf")
+  print "%s: %.17e" % (name, err)
+  assert(err < tol)
+  return
+
+# First order Adams-Bashforth
+
+levels    = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+levels_dT = TimeLevels(levels = [n], cycle_map = {},
+  last_past_level = n - 1)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[n]) * dx == dt * F(test, T[n]),
+  dT[n], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n]), (1.0, dT[n])),
+  T[n + 1])
+
+check("AB1", system, T[N], tol = 3.0e-15)
+
+# Second order Adams-Bashforth
+
+levels    = TimeLevels(levels = [n, n + 1, n + 2],
+  cycle_map = {n:n + 1, n + 1:n + 2}, last_past_level = n + 1)
+levels_dT = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1},
+  last_past_level = n)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[0]) * dx == dt * F(test, T[0]),
+  dT[0], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (1.0, dT[0])),
+  T[1])
+
+system.add_solve(inner(test, dT[n + 1]) * dx ==
+  dt * F(test, T[n + 1]),
+  dT[n + 1], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n + 1]),
+                                   (3.0 / 2.0, dT[n + 1]),
+                                   (-1.0 / 2.0, dT[n])),
+  T[n + 2])
+
+check("AB2", system, T[N + 1], ns_i = 1, tol = 4.0e-15)
+
+# Third order Adams-Bashforth
+
+from fractions import Fraction
+hf = Fraction(1, 2)
+
+levels    = TimeLevels(levels = [n, n + hf, n + 1, n + 2, n + 3],
+  cycle_map = {n:n + 1, n + 1:n + 2, n + 2:n + 3},
+  last_past_level = n + 2)
+levels_dT = TimeLevels(levels = [n, n + hf, n + 1, n + 2],
+  cycle_map = {n:n + 1, n + 1:n + 2}, last_past_level = n + 1)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[0]) * dx == dt * F(test, T[0]),
+  dT[0], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (0.5, dT[0])),
+  T[hf])
+
+system.add_solve(inner(test, dT[hf]) * dx == dt * F(test, T[hf]),
+  dT[hf], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (1.0, dT[hf])),
+  T[1])
+
+system.add_solve(inner(test, dT[1]) * dx == dt * F(test, T[1]),
+  dT[1], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[1]),
+                                   (3.0 / 2.0, dT[1]),
+                                   (-1.0 / 2.0, dT[0])),
+  T[2])
+
+system.add_solve(inner(test, dT[n + 2]) * dx ==
+  dt * F(test, T[n + 2]),
+  dT[n + 2], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n + 2]),
+                                   (23.0 / 12.0, dT[n + 2]),
+                                   (-4.0 / 3.0, dT[n + 1]),
+                                   (5.0 / 12.0, dT[n])),
+  T[n + 3])
+
+check("AB3", system, T[N + 2], ns_i = 2, tol = 6.0e-15)
+
+# Leapfrog
+
+levels    = TimeLevels(levels = [n, n + 1, n + 2],
+  cycle_map = {n:n + 1, n + 1:n + 2}, last_past_level = n + 1)
+levels_dT = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1},
+  last_past_level = n)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[0]) * dx == dt * F(test, T[0]),
+  dT[0], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (1.0, dT[0])),
+  T[1])
+
+system.add_solve(inner(test, dT[n + 1]) * dx ==
+  dt * F(test, T[n + 1]),
+  dT[n + 1], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n]),
+                                   (2.0, dT[n + 1])),
+  T[n + 2])
+
+check("LF", system, T[N + 1], ns_i = 1, tol = 4.0e-15)
+
+# Backward Euler
+
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, space, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, T[n + 1]) * dx ==
+  inner(test, T[n]) * dx + dt * F(test, T[n + 1]),
+  T[n + 1], solver_parameters = solver_parameters)
+
+check("AM1", system, T[N], tol = 3.0e-15)
+
+# Crank-Nicolson (Implicit trapezium rule)
+
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, space, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, T[n + 1]) * dx ==
+  inner(test, T[n]) * dx +
+  dt * (0.5 * F(test, T[n + 1]) + 0.5 * F(test, T[n])),
+  T[n + 1], solver_parameters = solver_parameters)
+
+check("AM2", system, T[N], tol = 3.0e-15)
+
+# Crank-Nicolson (Implicit midpoint rule)
+
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, space, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, T[n + 1]) * dx ==
+  inner(test, T[n]) * dx +
+  dt * F(test, 0.5 * T[n + 1] + 0.5 * T[n]),
+  T[n + 1], solver_parameters = solver_parameters)
+
+check("IMR", system, T[N], tol = 3.0e-15)
+
+# Second order explicit Runge-Kutta
+
+from fractions import Fraction
+hf = Fraction(1, 2)
+
+levels   = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+levels_F = TimeLevels(levels = [n, n + hf], cycle_map = {},
+  last_past_level = n - hf)
+T    = TimeFunction(levels,   space, name = "T")
+F_s  = TimeFunction(levels_F, space, name = "F_s")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, F_s[n]) * dx == F(test, T[n]),
+  F_s[n], solver_parameters = solver_parameters)
+system.add_solve(inner(test, F_s[n + hf]) * dx ==
+  F(test, T[n] + 0.5 * dt * F_s[n]),
+  F_s[n + hf], solver_parameters = solver_parameters)
+  
+system.add_solve(LinearCombination((1.0, T[n]),
+                                   (dt, F_s[n + hf])),
+  T[n + 1])
+
+check("RK2", system, T[N], tol = 5.0e-16)
+
+# P1_DG
+
+spaces = space * space
+tests = TestFunction(spaces)
+test1, test2 = split(tests)
+
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, spaces, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(inner(tests, T[0]) * dx == inner(test2, T_ic) * dx,
+  T[0], solver_parameters = solver_parameters)
+
+m_11 = m_22 = 1.0 / 3.0
+m_12 = m_21 = 1.0 / 6.0
+system.add_solve(
+  inner(test1, 0.5 * T[n + 1][1] + 0.5 * T[n + 1][0] - T[n][1])*dx +
+  inner(test2, 0.5 * T[n + 1][1] - 0.5 * T[n + 1][0]) * dx ==
+  dt*(m_11 * F(test1, T[n + 1][0]) + m_12 * F(test1, T[n + 1][1]) +
+      m_21 * F(test2, T[n + 1][0]) + m_22 * F(test2, T[n + 1][1])),
+  T[n + 1], solver_parameters = {"linear_solver":"lu"})
+  
+check("P1_DG", system, split(T[n])[1], tol = 1.0e-14)
\ No newline at end of file

=== added file 'timestepping/manual/examples/wrappers'
--- timestepping/manual/examples/wrappers	1970-01-01 00:00:00 +0000
+++ timestepping/manual/examples/wrappers	2013-05-15 14:57:23 +0000
@@ -0,0 +1,28 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from dolfin import *
+from timestepping import *
+
+# A named Constant
+c = Constant(1.0, name = "c")
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+# A named Function
+F = Function(space, name = "F")
\ No newline at end of file

=== added file 'timestepping/manual/makefile'
--- timestepping/manual/makefile	1970-01-01 00:00:00 +0000
+++ timestepping/manual/makefile	2013-05-15 14:57:23 +0000
@@ -0,0 +1,24 @@
+LATEXFLAGS=-halt-on-error
+BIBTEXFLAGS=
+
+TEX = $(wildcard *.tex)
+PDF = $(addsuffix .pdf, $(basename $(TEX)))
+
+default: $(PDF)
+
+%.pdf: %.tex bibliography.bib
+	pdflatex $(LATEXFLAGS) $<
+	bibtex $(BIBTEXFLAGS) $(basename $<) || true
+	pdflatex $(LATEXFLAGS) $<
+	pdflatex $(LATEXFLAGS) $<
+	pdflatex $(LATEXFLAGS) $<
+	pdflatex $(LATEXFLAGS) $<
+
+clean:
+	rm -f $(PDF)
+	rm -f $(addsuffix .aux, $(basename $(TEX)))
+	rm -f $(addsuffix .bbl, $(basename $(TEX)))
+	rm -f $(addsuffix .blg, $(basename $(TEX)))
+	rm -f $(addsuffix .log, $(basename $(TEX)))
+	rm -f $(addsuffix .out, $(basename $(TEX)))
+	rm -f $(addsuffix .toc, $(basename $(TEX)))

=== added file 'timestepping/manual/manual.tex'
--- timestepping/manual/manual.tex	1970-01-01 00:00:00 +0000
+++ timestepping/manual/manual.tex	2013-05-15 14:57:23 +0000
@@ -0,0 +1,1636 @@
+\documentclass[a4paper]{book}
+
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{array}
+\usepackage{graphicx}
+\usepackage{hyperref}
+\usepackage{listings}
+\usepackage{longtable}
+\usepackage{natbib}
+\usepackage{nicefrac}
+\usepackage{ulem}
+\usepackage[svgnames]{xcolor}
+
+\hypersetup{
+  colorlinks = true,
+  citecolor = DarkBlue,
+  linkcolor = DarkBlue,
+  urlcolor = DarkBlue
+}
+
+\lstset{
+  language = Python,
+  basicstyle = \ttfamily\footnotesize,
+  frame = single,
+  commentstyle = \color{Grey},
+  keywordstyle = \color{blue},
+  stringstyle = \color{DarkViolet},
+  keepspaces = true,
+  showspaces = false,
+  showstringspaces = false
+}
+
+\newcommand{\version}{1.2.0}
+
+\begin{document}
+
+\begin{titlepage}
+\begin{center}
+
+\Huge{\verb+timestepping+ Manual} \\[0.02\textheight]
+\large{The high-level representation of transient finite element models with DOLFIN} \\[0.02\textheight]
+\Huge{Version \version} \\[0.15\textheight]
+\vfill
+\large{\today}
+
+\end{center}
+\end{titlepage}
+
+\tableofcontents
+
+\chapter{Introduction}
+
+\section{Overview}
+
+The \verb+timestepping+ Python library is an extension for DOLFIN, adding a
+high-level representation for time dependent finite element models.
+Specifically, the \verb+timestepping+ library enables the rapid development of
+efficient timestepping finite element models using a symbolic representation of
+the entire model discretisation. The FEniCS system (for which DOLFIN acts as a
+front-end) uses automated code generation to convert a symbolic representation
+of a finite element spatial discretisation into a working model. The
+\verb+timestepping+ library builds on this principle, adding a syntax to enable
+a symbolic representation of a temporal discretisation.
+
+A machine interpretable representation of a discretisation is amenable to
+automated analysis and manipulation. In particular, if a symbolic representation
+of a model time discretisation is available then a number of optimisation
+strategies can be applied automatically. This enables time independent data,
+which may be expensive to recompute, to be computed once and cached ahead of the
+execution of the main model time loop. Moreover, the FEniCS system supplies
+symbolic manipulation tools which enable one to differentiate and adjoin
+spatially discretised equations. This functionality can be combined with a
+high-level representation of the temporal discretisation to enable the automatic
+derivation of an associated discrete adjoint model. Since the resulting adjoint
+model itself has a symbolic representation, it is also amenable to the
+application of optimisation strategies. Combined, the approach allows the entire
+model to be specified using only a symbolic representation of model equations,
+yields an efficient implementation of the forward model code, and enables the
+automated derivation and efficient implementation of an associated discrete
+adjoint model code.
+
+Since the \verb+timestepping+ library builds upon the FEniCS system, it inherits
+benefits associated with the application of automated code generation
+techniques. In particular, automated code generation can be used to separate the
+specification of a numerical model from the details of its implementation. The
+\verb+timestepping+ library inherits much of the portability of the FEniCS
+system: architectures and parallelisation strategies supported by the FEniCS
+system can be used in combination with the \verb+timestepping+ library.
+
+This manual describes the functionality and interfaces provided by the
+\linebreak \verb+timestepping+ Python library. This library makes extensive use
+of the DOLFIN Python library, but the functionality and interfaces provided by
+DOLFIN will not be described in detail here. For further details of DOLFIN and
+the FEniCS system the reader is instead referred to \citet{logg2012}. In
+addition, some familiarity with the Python language, the finite element method,
+simple time discretisations, and adjoint modelling is assumed.
+
+\section{General principles}\label{sect:principles}
+
+The \verb+timestepping+ library works on the principle that a timestepping
+model can be broken into three key stages:
+\begin{enumerate}
+  \item The initialisation stage: A series of equations which are solved
+        exactly once.
+  \item The timestepping stage: A series of equations which are solved
+        repeatedly, and potentially a very large number of times.
+  \item The finalisation stage: A series of equations which are solved
+        exactly once.
+\end{enumerate}
+The library further assumes that the timestepping stage can be sub-divided into
+two key stages:
+\begin{enumerate}
+  \item The timestep solve stage: A series of timestep equations are solved,
+        using earlier time data as input and yielding future time data as
+        output.
+  \item The timestep variable cycle stage: Past time data is replaced using
+        later time data.
+\end{enumerate}
+The interfaces provided by the library mimic this structure: the library
+provides a syntax for the description of each of these key stages. The
+\verb+timestepping+ library may not be suitable for problems which are
+incompatible with this representation.
+
+The implementation of the \verb+timestepping+ library makes a series of
+additional assumptions regarding the performance demands of the timestepping
+model. While these are not in principle essential to the development of time
+discretisation specific optimisations, they have nevertheless generally been
+applied in the development of the library. The assumptions are:
+\begin{enumerate}
+  \item Model initialisation is cheap.
+  \item Work performed on every model timestep is expensive.
+  \item Caching of data is preferable to recalculation.
+\end{enumerate}
+Since it is assumed that model initialisation is cheap, it is generally assumed
+that the cost associated with the manipulation and analysis of model equations
+is negligible. The timestepping model itself is optimised, but the process of
+generating the optimised model is not. It is further generally assumed that
+all work performed by a model timestep should be minimised, and therefore
+that aggressive caching strategies should be applied. It is assumed that
+sufficient memory is available for such caching\footnote{Storage of the forward
+model solution, for the purposes of performing an adjoint calculation, is a
+notable exception to this principle.}
+
+\section{Dependencies}
+
+The \verb+timestepping+ library has been tested with the following dependency
+versions:
+
+\begin{center}
+\begin{tabular}{| c | c | }
+\hline
+Dependency  & Tested dependency versions \\
+\hline
+Python  & 2.7.3, 2.7.4 \\
+\hline
+NumPy  & 1.6.1, 1.7.1 \\
+SciPy  & 0.9.0, 0.11.0 \\
+VTK    & 5.8.0 \\
+PETSc  & 3.1.0p8, 3.2.0p5 \\
+\hline
+DOLFIN   & 1.0.0, 1.1.1+, 1.2.0 \\
+FIAT     & 1.0.0, 1.1.0 \\
+FErari   & 1.0.0 \\
+FFC      & 1.0.0, 1.1.0, 1.2.0 \\
+Instant  & 1.0.0, 1.1.0, 1.2.0 \\
+SyFi     & 1.0 \\
+UFC      & 2.0.5, 2.1.1, 2.2.0 \\
+UFL      & 1.0.0, 1.1.0, 1.2.0 \\
+Viper    & 1.0.0 \\
+\hline
+\end{tabular}
+\end{center}
+
+\section{Library interfaces}
+
+The \verb+timestepping+ library can be divided into a primary high-level
+interface and a secondary low-level interface. The high-level interface requires
+interaction with high-level data types representing features such as model time
+levels, functions, and model equations, and provides a syntax enabling the
+description and implementation of a model with the structure discussed in
+section \ref{sect:principles}. The low-level interface provides additional
+advanced features, for example relating to the definition of explicit time
+dependency, but requires interaction with the internal structure of the
+\verb+timestepping+ library.
+
+This manual documents only the primary high-level interface. Further
+documentation relating to all functions and classes in the \verb+timestepping+
+library can be accessed via the online documentation. For example, in a Python
+shell:
+\begin{lstlisting}
+>>> import timestepping
+>>> help(timestepping)
+\end{lstlisting}
+will display the online documentation.
+
+\chapter{The high-level interface}
+
+This chapter details the primary high-level interface provided by the \linebreak
+\verb+timestepping+ library. This interface provides a syntax enabling the
+description and implementation of a model with the structure discussed in
+section \ref{sect:principles}.
+
+\section{Accessing the library}
+
+The \verb+timestepping+ library is intended to be used in combination with
+DOLFIN. The \verb+timestepping+ library overrides some aspects of the DOLFIN
+interface, and hence the standard way of accessing the library is via:
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+\end{lstlisting}
+
+\section{Changes to DOLFIN interfaces}\label{sect:wrappers}
+
+The \verb+timestepping+ library provides \verb+Constant+ and \verb+Function+
+functions, which wrap the DOLFIN \verb+dolfin.Constant+ and
+\verb+dolfin.Function+ constructors. The arguments accepted by these functions
+are identical to the arguments accepted by the DOLFIN constructors, except that
+each accepts an optional string \verb+name+ keyword argument defining the
+\verb+Constant+ or \verb+Function+ name.
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# A named Constant
+c = Constant(1.0, name = "c")
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+# A named Function
+F = Function(space, name = "F")
+\end{lstlisting}
+
+\section{Time level abstraction}
+
+Model time levels are represented in the following way:
+\begin{enumerate}
+  \item Initialisation stage time levels: Via an integer or
+        \verb+fractions.Fraction+.
+  \item Timestep stage time levels: Via a \verb+TimeLevel+ object.
+  \item Finalisation stage time levels: Via a \verb+FinalTimeLevel+ object.
+\end{enumerate}
+
+In the initialisation stage time levels are represented using an integer
+or a \verb+fractions.Fraction+. The value zero is used to indicate the
+transition point between the initialisation stage and the timestepping stage.
+
+In the timestepping stage time levels are represented using a \verb+TimeLevel+
+object. A \verb+TimeLevel+ has an associated integer or
+\verb+fractions.Fraction+ ``offset'', enabling the logical comparison of
+different \verb+TimeLevel+ objects. The \verb+timestepping+ library provides
+an abstract handle \verb+n+, which is a \linebreak \verb+TimeLevel+ assigned an
+offset value of zero.
+
+In the finalisation stage time levels are represented using a
+\verb+FinalTimeLevel+ object. A \verb+FinalTimeLevel+ has an associated integer
+or \verb+fractions.Fraction+ ``offset'' enabling the logical comparison of
+different \verb+FinalTimeLevel+ objects. The \verb+timestepping+ library
+provides an abstract handle \verb+N+, which is a \linebreak
+\verb+FinalTimeLevel+ assigned an offset value of zero. An offset value of zero
+is used to indicate the transition point between the timestepping stage and the
+finalisation stage.
+
+In the initialisation stage levels represented using an integer or \linebreak
+\verb+fractions.Fraction+ are well defined, but those represented using a
+\verb+TimeLevel+ or \verb+FinalTimeLevel+ are not. Immediately after
+initialisation levels represented using an integer, \verb+fractions.Fraction+,
+or \verb+TimeLevel+ are well defined, but those represented using a
+\verb+FinalTimeLevel+ are not. After the first timestep levels represented using
+a \verb+TimeLevel+ are well defined, but those represented using an integer,
+\verb+fractions.Fraction+, or \verb+FinalTimeLevel+ are not. After finalisation
+levels represented using a \verb+TimeLevel+ or \verb+FinalTimeLevel+ are well
+defined, but those represented using an integer or \verb+fractions.Fraction+ are
+not.
+
+Immediately after initialisation levels represented using an integer or
+\linebreak \verb+fractions.Fraction+ with value \verb+i+ are equivalent to
+levels represented using a \verb+TimeLevel+ with value \verb=n + i=. After
+finalisation levels represented using a \verb+TimeLevel+ with value \verb=n + i=
+are equivalent to those represented using a \linebreak \verb+FinalTimeLevel+
+with value \verb=N + i=.
+
+New \verb+TimeLevel+ and \verb+FinalTimeLevel+ objects can be instantiated
+directly. The high-level constructor arguments are:
+\begin{lstlisting}
+  n_custom = TimeLevel(arg)
+  N_custom = FinalTimeLevel(arg)
+\end{lstlisting}
+where the optional argument is an integer or \verb+fractions.Fraction+ defining
+the offset. If this argument is not supplied then a default value of zero is
+used. New \verb+TimeLevel+ and \verb+FinalTimeLevel+ objects can also be
+instantiated via the addition or subtraction of an integer or
+\verb+fractions.Fraction+ from an existing \verb+TimeLevel+ or
+\verb+FinalTimeLevel+:
+\begin{lstlisting}
+n_custom = n + offset_1
+N_custom = N - offset_2
+\end{lstlisting}
+
+\subsection*{Examples}
+\begin{lstlisting}
+from timestepping import *
+from fractions import Fraction
+
+# Create a TimeLevel with an offset of 2
+n1 = TimeLevel(2)
+# Create a TimeLevel with an offset of -1/2
+n2 = TimeLevel(-Fraction(1, 2))
+# Create a FinalTimeLevel with an offset of 0
+N1 = FinalTimeLevel()
+
+# Create a TimeLevel with an offset of 2
+n3 = n + 2
+# Create a TimeLevel with an offset of -3/2
+n4 = n - Fraction(3, 2)
+
+# Create a FinalTimeLevel with an offset of 1
+N2 = N + 1
+
+print n1 == n2  # False
+print n1 > n2   # True
+print n1 == n3  # True
+print n1 >= n3  # True
+
+print N1 == N2  # False
+print N1 <= N2  # True
+\end{lstlisting}
+
+\section{Function time levels}
+
+A \verb+TimeLevels+ object represents a unique set of model time levels,
+combined with a ``cycle map'' defining how data is to be re-assigned during the
+timestep variable cycle stage. A \verb+TimeLevels+ also defines the ``present''
+point, defining which levels are to be treated as being in the past and which
+levels are to be treated as being in the future.
+
+A \verb+TimeLevels+ object is instantiated via:
+\begin{lstlisting}
+time_levels = TimeLevels(levels, cycle_map, last_past_level)
+\end{lstlisting}
+where \verb+levels+ is a list containing \verb+TimeLevel+ objects,
+\verb+cycle_map+ is a dictionary with \verb+TimeLevel+ keys and values, and the
+optional \verb+last_past_level+ is a \verb+TimeLevel+.
+
+\verb+levels+ defines the model time levels which are to be represented by the
+\verb+TimeLevels+ object. The elements of \verb+levels+ must be unique.
+
+\verb+cycle_map+ defines the timestep variable cycle. At the end of a model
+timestep data corresponding to each key in \verb+cycle_map+ is replaced with
+data corresponding to the associated value in the \verb+cycle_map+. The
+replacement of data is performed in order, with earlier (lower offset) keys in
+\verb+cycle_map+ replaced before later keys. All keys in \verb+cycle_map+ must
+be earlier (have a lower offset) than their associated values. The elements in
+\verb+cycle_map.values()+ must be unique.
+
+The optional \verb+last_past_level+ defines the ``present'' point associated
+with the \verb+TimeLevels+. Time levels before or equal to
+\verb+last_past_level+ (with a lower or equal offset) are treated as being in
+the past, while other levels are treated as being in the future. If
+\verb+last_past_level+ is not supplied then a default value of \verb+n+ is used.
+All keys in \verb+cycle_map+ must be in the past.
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from timestepping import *
+from fractions import *
+
+# Create a TimeLevels with one past level and one future level, with
+# the past level data replaced by the future level data during the
+# timestep variable cycle. Suitable for use with forward Euler,
+# backward Euler, or Crank-Nicolson schemes.
+levels_1 = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+
+# Create a TimeLevels with one future level
+levels_2 = TimeLevels(levels = [n], cycle_map = {},
+  last_past_level = n - 1)
+
+# Create a TimeLevels with one past level and two future levels,
+# with the past level data replaced by the latest future level data
+# during the timestep variable cycle. Suitable for use with a
+# second order Runge-Kutta scheme.
+levels_3 = TimeLevels(levels = [n, n + Fraction(1, 2), n + 1],
+  cycle_map = {n:n + 1})
+
+# Create a TimeLevels with three past levels and one future level,
+# with the past levels replaced by the nearest later level during
+# the timestep variable cycle. Suitable for use with a third order
+# Adams-Bashforth scheme.
+levels_4 = TimeLevels(levels = [n - 2, n - 1, n, n + 1],
+  cycle_map = {n - 2:n - 1, n - 1:n, n:n + 1})
+\end{lstlisting}
+
+\section{Time dependent functions}
+
+A \verb+TimeFunction+ object defines a series of \verb+dolfin.Function+ objects,
+each specified on an individual time level.
+
+A \verb+TimeFunction+ object is instantiated via:
+\begin{lstlisting}
+tfn = TimeFunction(levels, *args, name = name, **kwargs)
+\end{lstlisting}
+where \verb+levels+ is a \verb+TimeLevels+, \verb+args+ and \verb+kwargs+ are
+arguments as accepted by a \verb+dolfin.Function+ constructor, and \verb+name+
+is an optional string representing the function name.
+
+A symbolic representation of a single function, defined on a single time level,
+can be accessed by indexing into the \verb+TimeFunction+ using an integer,
+\verb+fractions.Fraction+, \verb+TimeLevel+, or \verb+FinalTimeLevel+. For
+example:
+\begin{lstlisting}
+tfn[0]      # Representation of the function at time level 0
+tfn[n + 1]  # Representation of the function at time level n + 1
+tfn[N - 2]  # Representation of the function at time level N - 2
+\end{lstlisting}
+The index provided must match a level that, given the \verb+TimeLevels+ used to
+instantiate the \verb+TimeFunction+, is well defined at some point in the model
+execution.
+
+Prior to initialisation a \verb+dolfin.Function+ returned by indexing into the
+\verb+TimeFunction+ can be used for expressing equations symbolically, but
+cannot be used to access model data. Immediately after initialisation the
+\linebreak \verb+dolfin.Function+ returned by indexing into the
+\verb+TimeFunction+ with an integer, \verb+fractions.Fraction+ or
+\verb+TimeLevel+ can be used to access model data, provided the time level is in
+the past. After the first timestep the \verb+dolfin.Function+ returned by
+indexing into the \verb+TimeFunction+ with a \verb+TimeLevel+ can be used to
+access model data, provided the time level is in the past. After finalisation
+the \verb+dolfin.Function+ returned by indexing into the \verb+TimeFunction+
+with a \verb+TimeLevel+ or a \verb+FinalTimeLevel+ can be used to access model
+data. In all other cases, the data associated with \verb+dolfin.Function+
+objects returned by indexing into the \verb+TimeFunction+ are in an undefined
+state.
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit square
+mesh = UnitSquareMesh(10, 10)
+# P1 function space
+space_p1 = FunctionSpace(mesh, "CG", 1)
+# P2_{DG} function space
+space_p2dg = FunctionSpace(mesh, "DG", 2)
+
+# Define two different sets of time levels
+levels_1 = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+levels_2 = TimeLevels(levels = [n], cycle_map = {},
+  last_past_level = n - 1)
+
+# Define time dependent functions on the levels
+F1 = TimeFunction(levels_1, space_p1, name = "F1")
+F2 = TimeFunction(levels_2, space_p2dg, name = "F2")
+
+# An equation involving the two time dependent functions
+test_p1 = TestFunction(space_p1)
+eq = inner(test_p1, F1[n + 1]) * dx == inner(test_p1, F2[n]) * dx
+\end{lstlisting}
+
+\section{Model equations}
+
+A \verb+TimeSystem+ object is used to keep a record of model equations.
+Model equations can be registered with a \verb+TimeSystem+ by calling the
+\verb+add_assignment+ or \verb+add_solve+ methods.
+
+A \verb+TimeSystem+ is instantiated without arguments:
+\begin{lstlisting}
+system = TimeSystem()
+\end{lstlisting}
+
+\subsection{Registering assignments}
+
+An assignment can be registered via the \verb+add_assignment+ method:
+\begin{lstlisting}
+system.add_assignment(y, x)
+\end{lstlisting}
+Here \verb+x+ is a \verb+dolfin.Function+ corresponding to a time level of a
+\verb+TimeFunction+, and \verb+y+ is a float, a \verb+dolfin.Constant+, a
+\verb+dolfin.Function+, a \verb+LinearCombination+ (explained below), or a
+general \verb+ufl.expr.Expr+ (an arbitrary expression). This call corresponds to
+an equation in which \verb+x+ is set equal to the value of \verb+y+. The value
+of \verb+y+ must correspond to a either a constant value or a function defined
+in the function space of \verb+x+. The assignment must also be linear, with
+\verb+y+ independent of \verb+x+.
+
+A \verb+LinearCombination+ is instantiated via:
+\begin{lstlisting}
+lc = LinearCombination(*args)
+\end{lstlisting}
+where the arguments \verb+args+ are an arbitrary list of tuples
+\verb+(alpha, y)+. Each \verb+alpha+ is a float, a \verb+dolfin.Constant+, a
+\verb+dolfin.Function+, or a general \linebreak \verb+ufl.expr.Expr+. The
+\verb+LinearCombination+ then represents the value of $\sum_i \alpha_i y_i$.
+Note that the \verb+LinearCombination+ constructor does not check that all the
+\verb+alpha+ are independent of all the \verb+y+, and consequently does not
+check that the object represents a true linear combination. All the \verb+y+
+must be in the same function space, and all the \verb+alpha+ must correspond
+either to a constant value or a function defined in the function space of the
+\verb+y+.
+
+Constant real \verb+dolfin.Function+ objects (i.e. functions defined in a
+DOLFIN function space \verb+FunctionSpace(mesh, "R", 0)+) are, when used in
+general expressions (\verb+y+ in the \verb+add_solve+ method or an \verb+alpha+
+in the \verb+LinearCombination+ constructor), treated as corresponding to a
+single constant value.
+
+\subsection{Registering finite element variational problems}
+
+An equation representing a finite element spatial discretisation can be
+registered via the \verb+add_solve+ method. The arguments of this method are
+largely based upon the arguments accepted by the \verb+dolfin.solve+ function.
+The following are the high-level interfaces offered by this method:
+
+\begin{lstlisting}
+system.add_solve(y, x)
+\end{lstlisting}
+where \verb+y+ and \verb+x+ are as expected by the \verb+add_assignment+ method.
+This calls the \verb+add_assignment+ method.
+
+\begin{lstlisting}
+system.add_solve(a == L, x, bcs,
+  solver_parameters = solver_parameters,
+  adjoint_solver_parameters = adjoint_solver_parameters)
+\end{lstlisting}
+where \verb+a+ is a rank 2 \verb+dolfin.Form+, \verb+L+ is a rank 1 \linebreak
+\verb+dolfin.Form+, and \verb+x+ is a \verb+dolfin.Function+ corresponding to a
+time level of a \linebreak \verb+TimeFunction+. This registers a (possibly
+non-linear) variational problem. Alternatively, one may use:
+\begin{lstlisting}
+system.add_solve(eq, x, bcs,
+  solver_parameters = solver_parameters,
+  adjoint_solver_parameters = adjoint_solver_parameters)
+\end{lstlisting}
+where \verb+eq+ is a (possibly non-linear) \verb+dolfin.Equation+ depending on
+\verb+x+, and \verb+x+ is a \verb+dolfin.Function+ corresponding to a time level
+of a \verb+TimeFunction+. This registers a (possibly non-linear) variational
+problem represented by the equation \verb+eq+.
+
+If a registered variational problem is non-linear then the problem is solved
+using a Newton solver. The optional \verb+bcs+ is a \verb+dolfin.DirichletBC+ or
+a list of \verb+dolfin.DirichletBC+ objects, defining Dirichlet boundary
+conditions to be applied when solving the variational problem. The optional
+\verb+solver_parameters+ is a dictionary of DOLFIN solver parameters for use in
+solving the variational problem, while the optional
+\verb+adjoint_solver_parameters+ is a dictionary of DOLFIN solver parameters for
+use in solving an associated adjoint variational problem. If
+\verb+solver_parameters+ is not supplied then the default parameters defined in
+\verb+dolfin.parameters+ (at the time of the \verb+TimeSystem.assemble+ call --
+see section \ref{sect:assembly}) are used. If \verb+adjoint_solver_parameters+
+is not supplied then the forward solver parameters are used in solving an
+associated adjoint variational problem.
+
+\subsection{Permitted equation dependencies}
+
+The solution \verb+dolfin.Function+ of a registered equation, \verb+x+,
+corresponds either to an initialisation stage time level, a timestep stage time
+level, or a finalisation stage time level. Any time level dependencies of the
+equation must be in this stage. For example, any time level dependencies of an
+equation which solves for an initialisation stage time level must themselves
+correspond to an initialisation stage time level. The following is a valid
+registration:
+\begin{lstlisting}
+system.add_solve(F[-1], F[0])
+\end{lstlisting}
+(provide that \verb+F+ is defined on the relevant levels), while the following
+is an \emph{invalid} registration:
+\begin{lstlisting}
+system.add_solve(F[n], F[0])
+\end{lstlisting}
+
+If the solution \verb+dolfin.Function+ of a registered equation, \verb+x+,
+corresponds to an initialisation stage time level, then the time level must be
+in the past. If \verb+x+ corresponds to a timestep stage or finalisation stage
+time level, then the time level must be in the future. No dependencies of a
+registered equation can correspond to time levels later than \verb+x+.
+
+Note that the order of equation registration is irrelevant, and that
+circular dependencies are not permitted -- dependency resolution is applied in
+the analysis and optimisation step (see section \ref{sect:assembly}).
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = Constant(0.05)
+bc1 = DirichletBC(space, 1.0, "on_boundary && near(x[0], 0.0)")
+bc2 = DirichletBC(space, 0.0, "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = Constant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = Function(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
+\end{lstlisting}
+See also section \ref{sect:statics} for a modification of this example,
+declaring static data.
+
+\section{Static data}\label{sect:statics}
+
+Many optimisations of the timestepping model require the identification of
+static data, which does not vary throughout the model execution. Additional
+functions are provided to facilitate such optimisations.
+
+The \verb+StaticConstant+ function returns a \verb+dolfin.Constant+ that is
+declared as ``static'. The arguments are identical to the \verb+Constant+
+wrapper (see section \ref{sect:wrappers}). The \verb+StaticFunction+ function
+returns a \verb+dolfin.Function+ that is declared as ``static'. The arguments
+are identical to the \verb+Function+ wrapper (see section
+\ref{sect:wrappers}). The \verb+StaticDirichletBC+ class defines a
+\verb+dolfin.DirichletBC+ that is declared as ``static''. The constructor
+arguments are identical to the \verb+dolfin.DirichletBC+ constructor.
+
+Objects that are declared as ``static'' in this way should not be modified
+during the model execution -- doing so leads to undefined behaviour. Objects
+declared as static may be modified before the initialisation stage.
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# Create a constant which is declared as static
+c = StaticConstant(1.0, name = "c")
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+# Create a function which is declared as static
+F = StaticFunction(space, name = "F")
+# Initialise the static function
+F.interpolate(Expression("x[0]"))
+
+# Create a Dirichlet boundary condition which is declared as static
+bc = StaticDirichletBC(space, 0.0, "on_boundary")
+\end{lstlisting}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = StaticConstant(0.05)
+bc1 = StaticDirichletBC(space, 1.0,
+  "on_boundary && near(x[0], 0.0)")
+bc2 = StaticDirichletBC(space, 0.0,
+  "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = StaticConstant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = StaticFunction(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
+\end{lstlisting}
+
+\section{Analysis and optimisation}\label{sect:assembly}
+
+The \verb+assemble+ method of a \verb+TimeSystem+ analyses the model, performs
+time discretisation specific optimisations, and (by default) performs the
+initialisation stage, returning an object suitable for execution of the model
+time loop. An adjoint model is enabled by supplying additional arguments to the
+\verb+assemble+ method.
+
+The \verb+assemble+ method returns an \verb+ForwardModel+ if no adjoint model is
+enabled, and an \verb+ManagedModel+ otherwise. Note that there is a fairly
+strong assumption that at most one \verb+ForwardModel+ or
+\verb+ManagedModel+ is instantiated at once. Instantiating two or more at the
+same time leads to undefined behaviour.
+
+\subsection{Forward model only}\label{sect:forward_assembly}
+
+If an adjoint model is not enabled, the \verb+assemble+ method has high-level
+arguments:
+\begin{lstlisting}
+asystem = system.assemble(initialise = initialise)
+\end{lstlisting}
+This returns an \verb+ForwardModel+. If the optional \verb+initialise+ is not
+supplied or is \verb+True+ then, immediately after the \verb+assemble+ call, the
+initialisation stage is complete. Otherwise, the initialisation stage can be
+performed via a subsequent call to the \verb+initialise+ method:
+\begin{lstlisting}
+asystem.initialise()
+\end{lstlisting}
+
+\subsection{Enabling an adjoint model}\label{sect:adjoint_assembly}
+
+An adjoint model is enabled by calling the \verb+assemble+ method with the
+optional \verb+adjoint+ boolean argument set equal to \verb+True+. The remaining
+high-level arguments are:
+\begin{lstlisting}
+asystem = system.assemble(adjoint = adjoint,
+  disk_period = disk_period, functional = functional,
+  initialise = initialise)
+\end{lstlisting}
+This returns an \verb+ManagedModel+. If the optional \verb+initialise+ is not
+supplied or is \verb+True+ then, immediately after the \verb+assemble+ call, the
+initialisation stage is complete. Otherwise, the initialisation stage can be
+performed via a subsequent call to the \verb+initialise+ method:
+\begin{lstlisting}
+asystem.initialise()
+\end{lstlisting}
+The optional integer \verb+disk_period+ defines how often the forward model
+solution is to be checkpointed to disk. Data is checkpointed every
+\verb+disk_period+ timesteps, and is stored in a directory \verb+checkpoints~+
+relative to the current working directory. If \verb+disk_period+ is not supplied
+then the entire forward model solution is stored in memory. In the high-level
+interface the optional \verb+functional+ is a rank 0 form.
+\verb+dolfin.Function+ dependencies of the functional which correspond to levels
+of a \verb+TimeFunction+ must correspond to finalisation stage time levels. If
+\verb+functional+ is not supplied to the \verb+assemble+ method then a
+functional can instead be specified via the \verb+set_functional+ method:
+\begin{lstlisting}
+asystem.set_functional(functional)
+\end{lstlisting}
+
+\section{Timestepping}
+
+The model is timestepped by calling the \verb+timestep+ method of an \linebreak
+\verb+ForwardModel+ or \verb+ManagedModel+:
+\begin{lstlisting}
+asystem.timestep(s = s)
+\end{lstlisting}
+where the optional integer \verb+s+ is the number of timesteps to perform. If
+\verb+s+ is not supplied then a single timestep is performed. The
+\verb+timestep+ method performs both the timestep solve and timestep cycle
+stages. At the model end the finalisation stage is performed by calling the
+\verb+finalise+ method:
+\begin{lstlisting}
+asystem.finalise()
+\end{lstlisting}
+
+The initialisation stage should be complete before calling the \verb+timestep+
+method (see section \ref{sect:assembly}). If an adjoint model is not enabled
+then the model can, after a \verb+finalise+ call, be restarted via a call to the
+\verb+initialise+ method (see section \ref{sect:forward_assembly}). In this case
+the \verb+timestep+ and \verb+finalise+ methods should not be called between the
+\verb+finalise+ and \verb+initialise+ calls. If the model is not restarted, or
+if an adjoint model is enabled, then the \verb+timestep+ and \verb+finalise+
+methods should not be called after the \verb+finalise+ call.
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = StaticConstant(0.05)
+bc1 = StaticDirichletBC(space, 1.0,
+  "on_boundary && near(x[0], 0.0)")
+bc2 = StaticDirichletBC(space, 0.0,
+  "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = StaticConstant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = StaticFunction(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
+
+# Assemble the TimeSystem
+system = system.assemble()
+
+# Timestep the model
+t = 0.0
+while t * (1.0 + 1.0e-9) < 1.0:
+  system.timestep()
+  t += float(dt)
+# Finalise
+system.finalise()
+\end{lstlisting}
+
+\section{Adjoint modelling}
+
+Given a forward model described using a \verb+TimeSystem+, the associated
+discrete adjoint model can be derived in the analysis and optimisation step (see
+section \ref{sect:adjoint_assembly}). The resulting \verb+ManagedModel+
+returned can be used to perform a model constrained derivative calculation.
+
+\subsection{Checkpointing and storage}
+
+An adjoint model calculation requires knowledge of the forward model solution
+and, as a result, forward model data must be stored. It is often infeasible to
+store all required data in memory, and hence disk checkpointing and recovery
+must be used. In this procedure forward model data, sufficient to restart the
+forward model, are checkpointed periodically, and forward model data required
+between checkpoints are re-computed and stored. Forward model checkpointing and
+storage is configured in the analysis and optimisation step (see section
+\ref{sect:adjoint_assembly}).
+
+After the forward model finalisation stage (after the \linebreak
+\verb+ManagedModel.finalise+ call) stored forward model data can be verified via
+the \verb+verify_checkpoints+ method:
+\begin{lstlisting}
+asystem.verify_checkpoints(tolerance = tolerance)
+\end{lstlisting}
+where the optional \verb+tolerance+ argument is a positive float defining the
+tolerance used in comparisons. The \verb+verify_checkpoints+ method raises an
+exception on failure. For serial calculations the forward model data should be
+exactly reproducible, and verification should succeed with a tolerance of zero.
+For parallel calculations the verification may fail with a zero or small value
+for the tolerance. Note that the \verb+verify_checkpoints+ method re-runs the
+entire forward model, and hence this should only be used for debugging purposes.
+
+\subsection{Derivative calculation}
+
+A model constrained derivative can be computed via the \verb+compute_gradient+
+method. The high-level arguments for this method are:
+\begin{lstlisting}
+grad = asystem.compute_gradient(parameters = parameters,
+  project = project,
+  project_solver_parameters = project_solver_parameters)
+\end{lstlisting}
+where \verb+parameters+ is a \verb+dolfin.Constant+, a \verb+dolfin.Function+,
+or a list of \verb+dolfin.Constant+ or \verb+dolfin.Function+ objects, the
+optional \verb+project+ is a \linebreak boolean, and the optional
+\verb+project_solver_parameters+ is a dictionary of solver parameters. If
+\verb+parameters+ is a \verb+dolfin.Constant+ or \verb+dolfin.Function+ then it
+must have been declared as static (see section \ref{sect:statics}). If
+\verb+parameters+ is a list of \verb+dolfin.Constant+ or \verb+dolfin.Function+
+objects then all elements must have been declared as static. This method returns
+the model constrained derivative of the registered functional (see section
+\ref{sect:adjoint_assembly}). If \verb+parameters+ is a list, then the
+\verb+i+th element of the return value contains the derivative with respect to
+the \verb+i+th element of \verb+parameters+. The derivative with respect to a
+\verb+Constant+ parameter is returned as a \verb+Constant+. If \verb+project+ is
+\verb+False+ (the default) then the derivative with respect to a \verb+Function+
+parameter is returned as a \verb+GenericVector+. If \verb+project+ is
+\verb+True+ then the derivative with respect to a \verb+Function+ parameter is
+returned as a \verb+(GenericVector, Function)+ pair, where the \verb+Function+
+contains the derivative projected onto the parameter trial space, and where
+\verb+project_solver_parameters+ defines solver options for the associated mass
+matrix inversion.
+
+The \verb+compute_gradient+ method must be called after the forward model
+finalisation stage (after the \verb+ManagedModel.finalise+ call), and a
+functional must be registered.
+
+\subsection{Adjoint verification}
+
+Let $\vec{m}$ correspond to the degrees of freedom associated with a parameter,
+and define a functional
+$J \left( \vec{x} \left( \vec{m} \right), \vec{m} \right)$, where $\vec{x}$
+corresponds to the degrees of freedom associated with the forward model solution.
+Then:
+\begin{equation}
+  \left| J \left( \vec{x} \left( \vec{m} + \epsilon \vec{\delta m} \right), \vec{m} + \epsilon \vec{\delta m} \right)
+    - J \left( \vec{x} \left( \vec{m} \right), \vec{m} \right) \right|
+    = {\cal O} \epsilon,
+\end{equation}
+converges asymptotically at first order in perturbations of the parameter
+$\vec{m}$. However:
+\begin{equation}
+  \left| J \left( \vec{x} \left( \vec{m} + \epsilon \vec{\delta m} \right), \vec{m} + \epsilon \vec{\delta m} \right)
+    - J \left( \vec{x} \left( \vec{m} \right), \vec{m} \right)
+    - \epsilon \left< \frac{d J}{d \vec{m}}, \vec{\delta m} \right> \right|
+    = {\cal O} \epsilon^2,
+\end{equation}
+converges asymptotically at second order. Hence the result of a model
+constrained derivative calculation, yielding $d J / d \vec{m}$, can be verified
+using repeated forward model calculations, each with a small perturbation of the
+parameter $\vec{m}$.
+
+The \verb+taylor_test+ method performs such a verification. The high-level
+arguments for this method are:
+\begin{lstlisting}
+orders = asystem.taylor_test(parameter, grad = grad,
+  ntest = ntest, fact  = fact)
+\end{lstlisting}
+\verb+parameter+ is a \verb+dolfin.Constant+ or \verb+dolfin.Function+ that has
+been declared as static (see section \ref{sect:statics}) and the optional
+\verb+grad+ is the corresponding derivative returned by the
+\verb+compute_gradient+ method. If \verb+grad+ is not supplied then the
+model constrained derivative is computed using the adjoint model. The optional
+\verb+ntest+ argument is the number of forward model calculations to perform in
+the verification. This must be greater than or equal to 2, and a default value
+of 6 is used if this is not supplied. The optional \verb+fact+ is a positive
+float and sets the magnitude of the perturbations. Let $N$ be the value of
+\verb+ntest+ and \verb+f+ be the value of \verb+fact+. If \verb+parameter+ is a
+\verb+dolfin.Constant+ with value $m$ then the perturbations applied are
+$\epsilon_i \delta m$, where:
+\begin{subequations}
+  \begin{equation}
+    \epsilon_i = 2^{N - i}, \quad i \in \left\{1, \ldots, N \right\},
+  \end{equation}
+  \begin{equation}
+    \delta m = \left\{ \begin{array}{l} f \left| m \right| \quad \textrm{if } \left| m \right| > 0 \\ f \qquad \quad \textrm{otherwise} \end{array} \right. .
+  \end{equation}
+\end{subequations}
+If \verb+parameter+ is a \verb+dolfin.Function+ with degrees of freedom
+with associated values $\vec{m}$ then the perturbations applies are
+$\epsilon_i \vec{\delta m}$, where:
+\begin{subequations}
+  \begin{equation}
+    \epsilon_i = 2^{N - i}, \quad i \in \left\{1, \ldots, N \right\},
+  \end{equation}
+  \begin{equation}
+    \vec{\delta m} = \left\{ \begin{array}{l} \vec{F} \left|\left| \vec{m} \right|\right|_\infty \quad \textrm{if } \left|\left| \vec{m} \right|\right|_\infty > 0 \\ \vec{F} \qquad \qquad \textrm{otherwise} \end{array} \right.
+  \end{equation}
+\end{subequations}
+where the values of $\vec{F}$ are randomly generated values in the range
+$\left[ -f/2, +f/2 \right)$, with a uniform probability distribution. If
+\verb+fact+ is not supplied then a default value of $10^{-4}$ is used.
+
+The \verb+taylor_test+ method returns a list of relative orders of
+convergence. The elements of the returned list have values:
+\begin{equation}
+  o_i = \textrm{log}_2 \frac{r_{i + 1}}{r_i},
+\end{equation}
+where:
+\begin{subequations}
+  \begin{equation}
+    \epsilon_i = 2^{N - i}, \quad i \in \left\{1, \ldots, N \right\},
+  \end{equation}
+  \begin{equation}
+    r_i = \left| J \left( \vec{x} \left( \vec{m} + \epsilon_i \vec{\delta m} \right), \vec{m} + \epsilon_i \vec{\delta m} \right)
+      - J \left( \vec{x} \left( \vec{m} \right), \vec{m} \right)
+      - \epsilon_i \left< \frac{d J}{d \vec{m}}, \vec{\delta m} \right> \right|.
+  \end{equation}
+\end{subequations}
+The elements of the returned list should be close to 2 in a successful
+verification.
+
+It is recommended that new adjoint models should be verified using a Taylor
+remainder convergence test. Since the test involves repeated running of the
+forward model it may be necessary to test the adjoint on a small problem,
+for example using a reduced number of model timesteps. The Taylor remainder
+convergence test may fail if:
+\begin{enumerate}
+  \item The perturbations are too large, meaning that the asymptotic convergence
+        of the remainder is not observed.
+  \item The perturbations are too small, meaning that the convergence test is
+        polluted by machine precision errors.
+  \item The forward model is not differentiable, meaning that a model
+        constrained derivative does not exists.
+  \item Iterative solver tolerances are too large.
+  \item The adjoint model is inconsistent with the forward model.
+\end{enumerate}
+The first two cases can be avoided by decreasing or increasing \verb+fact+
+(assuming that these two cases do not both apply simultaneously, in which
+case a Taylor remainder convergence test cannot be used). The final case may
+occur if forward model data is modified manually during the model execution, or
+may indicate a bug in the \verb+timestepping+ library.
+
+\subsection*{Examples}
+
+\begin{lstlisting}
+from dolfin import *
+from timestepping import *
+
+# Define a simple structured mesh on the unit interval
+mesh = UnitIntervalMesh(10)
+# P1 function space
+space = FunctionSpace(mesh, "CG", 1)
+
+# Model parameters and boundary conditions
+dt = StaticConstant(0.05)
+bc1 = StaticDirichletBC(space, 1.0,
+  "on_boundary && near(x[0], 0.0)")
+bc2 = StaticDirichletBC(space, 0.0,
+  "on_boundary && near(x[0], 1.0)")
+bcs = [bc1, bc2]
+nu = StaticConstant(0.01)
+
+# Define time levels
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+# A time dependent function
+u = TimeFunction(levels, space, name = "u")
+
+# Initialise a TimeSystem
+system = TimeSystem()
+
+# Add an initial assignment
+u_ic = StaticFunction(space, name = "u_ic")
+u_ic.assign(Constant(0.0))
+bc1.apply(u_ic.vector())
+system.add_solve(u_ic, u[0])
+# Register a simple diffusion equation, discretised in time
+# using forward Euler
+test = TestFunction(space)
+system.add_solve(
+  inner(test, (1.0 / dt) * (u[n + 1] - u[n])) * dx ==
+    -nu * inner(grad(test), grad(u[n])) * dx,
+  u[n + 1], bcs,
+  solver_parameters = {"linear_solver":"lu"})
+
+# Assemble the TimeSystem, enabling the adjoint. Set the
+# functional to be equal to spatial integral of the final u.
+system = system.assemble(adjoint = True, functional = u[N] * dx)
+
+# Timestep the model
+t = 0.0
+while t * (1.0 + 1.0e-9) < 1.0:
+  system.timestep()
+  t += float(dt)
+# Finalise
+system.finalise()
+
+# Perform a model constrained derivative calculation
+dJdm = system.compute_gradient(nu)
+
+# Verify the stored forward model data
+system.verify_checkpoints()
+# Verify the computed derivative using a Taylor remainder
+# convergence test
+orders = system.taylor_test(nu, grad = dJdm)
+# Check the convergence order
+assert((orders > 1.99).all())
+\end{lstlisting}
+
+\chapter{Time discretisation examples}
+
+This chapter describes the implementation of a number of simple time
+discretisations using the \verb+timestepping+ library. The provided examples are
+largely illustrative, and may not be stable.
+
+All examples in this chapter consider
+the discretisation of the following system:
+\begin{subequations}
+  \begin{equation}\label{eqn:t_timestep}
+    \partial_t T = F(T),
+  \end{equation}
+  \begin{equation}\label{eqn:t_init}
+    \left. T \right|_{t = 0} = T_0.
+  \end{equation}
+\end{subequations}
+It is assumed that $T_0$ is defined using a DOLFIN variable \verb+T_0+, which
+may for example be a \verb+dolfin.Constant+ or \verb+dolfin.Function+. It is
+further assumed that some Python function for approximating
+$\int_\Omega \phi_i F(T)$ over the space $\Omega$ and with test functions
+$\phi_i$ is provided, with the following form:
+\begin{lstlisting}
+def F(test, T):
+  rhs = ... # Insert definition here
+  return rhs
+\end{lstlisting}
+where \verb+test+ is a spatial test function.
+
+\section{Adams-Bashforth}
+
+Adams-Bashforth schemes take the form:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t \sum_{i = 1}^N \gamma_i F( T^{n + 1 - i}),
+\end{equation}
+where the $\gamma_i$ are constants. Adams-Bashforth schemes therefore require
+the values of $F(T)$ from $N$ previous time levels. Given a choice of $N$, the
+$\gamma_i$ are chosen so as to yield $N$th order accuracy. If $N > 1$ then the
+scheme is not self starting, and care is required when initialising in order to
+maintain the accuracy of the method. Further details of Adams-Bashforth schemes
+can be found in \citet[section 2.1]{iserles2009} and
+\citet[section III.1]{hairer1987}.
+
+\subsection{Forward Euler}\label{sect:fe}
+
+For $N = 1$ the single coefficient $\gamma_1$ is chosen to be:
+\begin{equation}
+  \gamma_1 = 1,
+\end{equation}
+yielding:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t F(T^n).
+\end{equation}
+This is the forward Euler scheme, and has first order accuracy.
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+levels    = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+levels_dT = TimeLevels(levels = [n], cycle_map = {},
+  last_past_level = n - 1)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[n]) * dx == dt * F(test, T[n]),
+  dT[n], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n]), (1.0, dT[n])),
+  T[n + 1])
+\end{lstlisting}
+where \verb+space+ defines the function space for $T$, \verb+dt+ defines the
+timestep size, and \verb+solver_parameters+ is a dictionary of solver
+parameters.
+
+\subsection{Second order Adams-Bashforth}\label{sect:ab2}
+
+For $N = 2$ the $\gamma_i$ are chosen to be:
+\begin{align}
+  \gamma_1 = \frac{3}{2}, & & \gamma_2 = -\frac{1}{2},
+\end{align}
+yielding:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t \left[ \frac{3}{2} F(T^n) - \frac{1}{2} F(T^{n - 1}) \right].
+\end{equation}
+This scheme cannot be used to perform the first model timestep as $T^{-1}$ is
+not available. In order to ensure global second order accuracy in time the
+method used to perform the first timestep, and compute $T^1$, must be at least
+first order accurate. Hence, for example, a single forward Euler step
+\ref{sect:fe} can be applied:
+\begin{equation}
+  T^1 = T^0 + \Delta t F(T^0).
+\end{equation}
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+levels    = TimeLevels(levels = [n, n + 1, n + 2],
+  cycle_map = {n:n + 1, n + 1:n + 2}, last_past_level = n + 1)
+levels_dT = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1},
+  last_past_level = n)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[0]) * dx == dt * F(test, T[0]),
+  dT[0], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (1.0, dT[0])),
+  T[1])
+
+system.add_solve(inner(test, dT[n + 1]) * dx ==
+  dt * F(test, T[n + 1]),
+  dT[n + 1], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n + 1]),
+                                   (3.0 / 2.0, dT[n + 1]),
+                                   (-1.0 / 2.0, dT[n])),
+  T[n + 2])
+\end{lstlisting}
+At the end of the simulation the final model solution is stored in
+\verb=T[N + 1]=.
+
+\subsection{Third order Adams-Bashforth}
+
+For $N = 3$ the $\gamma_i$ are chosen to be:
+\begin{align}
+  \gamma_1 = \frac{23}{12}, & & \gamma_2 = -\frac{4}{3}, & & \gamma_3 = \frac{5}{12},
+\end{align}
+yielding:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t \left[ \frac{23}{12} F(T^n) - \frac{4}{3} F(T^{n - 1}) + \frac{5}{12} F(T^{n - 2}) \right].
+\end{equation}
+This scheme cannot be used to perform the first two model timesteps as $T^{-1}$
+and $T^{-2}$ are not available. In order to ensure global third order accuracy
+in time the method used to perform the first two timesteps, and compute $T^0$
+and $T^1$, must be at least second order accurate. Hence, for example, second
+order Runge-Kutta \ref{sect:rk2} and second order Adams-Bashforth \ref{sect:ab2}
+steps can be applied:
+\begin{subequations}
+  \begin{equation}
+    T^{\frac{1}{2}} = T^0 + \frac{1}{2} \Delta t F(T^0).
+  \end{equation}
+  \begin{equation}
+    T^1 = T^0 + \Delta t F(T^{\frac{1}{2}}).
+  \end{equation}
+  \begin{equation}
+    T^2 = T^1 + \Delta t \left[ \frac{3}{2} F(T^1) - \frac{1}{2} F(T^0) \right].
+  \end{equation}
+\end{subequations}
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+from fractions import Fraction
+hf = Fraction(1, 2)
+
+levels    = TimeLevels(levels = [n, n + hf, n + 1, n + 2, n + 3],
+  cycle_map = {n:n + 1, n + 1:n + 2, n + 2:n + 3},
+  last_past_level = n + 2)
+levels_dT = TimeLevels(levels = [n, n + hf, n + 1, n + 2],
+  cycle_map = {n:n + 1, n + 1:n + 2}, last_past_level = n + 1)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[0]) * dx == dt * F(test, T[0]),
+  dT[0], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (0.5, dT[0])),
+  T[hf])
+
+system.add_solve(inner(test, dT[hf]) * dx == dt * F(test, T[hf]),
+  dT[hf], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (1.0, dT[hf])),
+  T[1])
+
+system.add_solve(inner(test, dT[1]) * dx == dt * F(test, T[1]),
+  dT[1], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[1]),
+                                   (3.0 / 2.0, dT[1]),
+                                   (-1.0 / 2.0, dT[0])),
+  T[2])
+
+system.add_solve(inner(test, dT[n + 2]) * dx ==
+  dt * F(test, T[n + 2]),
+  dT[n + 2], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n + 2]),
+                                   (23.0 / 12.0, dT[n + 2]),
+                                   (-4.0 / 3.0, dT[n + 1]),
+                                   (5.0 / 12.0, dT[n])),
+  T[n + 3])
+\end{lstlisting}
+At the end of the simulation the final model solution is stored in
+\verb=T[N + 2]=.
+
+\section{Leapfrog}
+
+A centred discretisation in time leads to:
+\begin{equation}
+  T^{n + 1} = T^{n - 1} + 2 \Delta t F( T^n ).
+\end{equation}
+This method is also known as the leapfrog scheme, and has second order accuracy.
+This scheme cannot be used to perform the first model timestep as $T^{-1}$ is
+not available. In order to ensure global second order accuracy in time the
+method used to perform the first timestep, and compute $T^1$, must be at least
+first order accurate. Hence, for example, a single forward Euler step
+\ref{sect:fe} can be applied:
+\begin{equation}
+  T^1 = T^0 + \Delta t F(T^0).
+\end{equation}
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+levels    = TimeLevels(levels = [n, n + 1, n + 2],
+  cycle_map = {n:n + 1, n + 1:n + 2}, last_past_level = n + 1)
+levels_dT = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1},
+  last_past_level = n)
+T  = TimeFunction(levels,    space, name = "T")
+dT = TimeFunction(levels_dT, space, name = "dT")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, dT[0]) * dx == dt * F(test, T[0]),
+  dT[0], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[0]),
+                                   (1.0, dT[0])),
+  T[1])
+
+system.add_solve(inner(test, dT[n + 1]) * dx ==
+  dt * F(test, T[n + 1]),
+  dT[n + 1], solver_parameters = solver_parameters)
+system.add_solve(LinearCombination((1.0, T[n]),
+                                   (2.0, dT[n + 1])),
+  T[n + 2])
+\end{lstlisting}
+At the end of the simulation the final model solution is stored in
+\verb=T[N + 1]=.
+
+\section{Adams-Moulton}
+
+Adams-Moulton schemes take the form:
+\begin{equation}\label{eqn:adams_moulton}
+  T^{n + 1} = T^n + \Delta t \sum_{i = 0}^N \gamma_i F( T^{n + 1 - i}),
+\end{equation}
+where the $\gamma_i$ are constants. For $n > 1$ Adams-Moulton schemes therefore
+require the values of $F(T)$ from $N - 1$ previous time levels. Given a choice
+of $N$, the $\gamma_i$ are chosen so as to yield $N$th order accuracy. If
+$N > 2$ then the scheme is not self starting, and care is required when
+initialising in order to maintain the accuracy of the method. Further details of
+Adams-Moulton schemes can be found in \citet[section III.1]{hairer1987}.
+
+In an Adams-Moulton scheme the right-hand-side of equation
+\eqref{eqn:adams_moulton} depends implicitly on the future solution $T^{n + 1}$.
+If $F(T)$ is linear then the implicit term contributes to the left-hand-side
+matrix in a linear solve for $T^{n + 1}$. If $F(T)$ is non-linear then
+a non-linear equation for $T^{n + 1}$ must be solved.
+
+\subsection{Backward Euler}\label{sect:be}
+
+For $N = 0$ the single coefficient $\gamma_0$ is chosen to be:
+\begin{equation}
+  \gamma_0 = 1,
+\end{equation}
+yielding:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t F(T^{n + 1}).
+\end{equation}
+This is the backward Euler scheme, and has first order accuracy.
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, space, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, T[n + 1]) * dx ==
+  inner(test, T[n]) * dx + dt * F(test, T[n + 1]),
+  T[n + 1], solver_parameters = solver_parameters)
+\end{lstlisting}
+Note that \verb+solver_parameters+ should define Newton solver parameters if
+$F(T)$ is non-linear.
+
+\subsection{Implicit trapezium rule}\label{sect:implicit_trapezium}
+
+For $N = 1$ the $\gamma_i$ are chosen to be:
+\begin{align}
+  \gamma_0 = \frac{1}{2}, & & \gamma_1 = \frac{1}{2},
+\end{align}
+yielding:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t \left[ \frac{1}{2} F(T^{n + 1}) + \frac{1}{2} F(T^n) \right].
+\end{equation}
+This is the implicit trapezium rule, and has second order accuracy. If $F(T)$ is
+linear then this scheme is identical to the implicit midpoint rule
+\ref{sect:implicit_midpoint}. This method is also known as the Crank-Nicolson
+scheme \citep[e.g.][section 3.4]{donea2003}, (although, confusingly,
+``Crank-Nicolson'' can refer to the implicit midpoint rule).
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, space, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, T[n + 1]) * dx ==
+  inner(test, T[n]) * dx +
+  dt * (0.5 * F(test, T[n + 1]) + 0.5 * F(test, T[n])),
+  T[n + 1], solver_parameters = solver_parameters)
+\end{lstlisting}
+Note that \verb+solver_parameters+ should define Newton solver parameters if
+$F(T)$ is non-linear.
+
+\section{Implicit midpoint rule}\label{sect:implicit_midpoint}.
+
+The implicit midpoint rule takes the form:
+\begin{equation}
+  T^{n + 1} = T^n + \Delta t F \left( \frac{1}{2} T^{n + 1} + \frac{1}{2} T^n \right).
+\end{equation}
+This has second order accuracy. If $F(T)$ is linear then this scheme is
+identical to the implicit trapezium rule \ref{sect:implicit_trapezium}. This
+method is also known as the Crank-Nicolson scheme
+\citep[e.g.][section 2.22]{mitchell1980} (although, confusingly,
+``Crank-Nicolson'' can refer to the implicit trapezium rule).
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, space, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, T[n + 1]) * dx ==
+  inner(test, T[n]) * dx +
+  dt * F(test, 0.5 * T[n + 1] + 0.5 * T[n]),
+  T[n + 1], solver_parameters = solver_parameters)
+\end{lstlisting}
+Note that \verb+solver_parameters+ should define Newton solver parameters if
+$F(T)$ is non-linear.
+
+\section{Explicit Runge-Kutta}
+
+Runge-Kutta schemes are a class of multi-stage integration schemes. In the
+context of equation \eqref{eqn:t_timestep} an $S$-stage explicit Runge-Kutta
+scheme takes the form:
+\begin{subequations}
+  \begin{equation}
+    F^{n + 1,i} = F \left( T^n + \Delta t \sum_{j = 1}^{i - 1} \gamma_{ij} F^{n + 1,j} \right) \quad \textrm{ for } i \in \left\{ 1, \ldots, S \right\},
+  \end{equation}
+  \begin{equation}
+    T^{n + 1} = T^n + \Delta t \sum_{j = 1}^S b_j F^{n + 1,j}.
+  \end{equation}
+\end{subequations}
+The $\gamma_{ij}$ and $b_j$ are chosen so as to yield desired accuracy and
+stability properties. Further details of explicit Runge-Kutta schemes can be
+found in \citet[section 3.2]{iserles2009}.
+
+\subsection{First order explicit Runge-Kutta}
+
+The one stage first order explicit Runge-Kutta scheme is:
+\begin{subequations}
+  \begin{equation}
+    F^{n + 1,1} = F \left( T^n \right),
+  \end{equation}
+  \begin{equation}
+    T^{n + 1} = T^n + \Delta t F^{n + 1,1}.
+  \end{equation}
+\end{subequations}
+This is the forward Euler scheme \ref{sect:fe}.
+
+\subsection{Second order explicit Runge-Kutta}\label{sect:rk2}
+
+A simple second order explicit Runge-Kutta scheme takes the form:
+\begin{subequations}
+  \begin{equation}
+    F^{n + 1,1} = F \left( T^n \right),
+  \end{equation}
+  \begin{equation}
+    F^{n + 1,2} = F \left( T^n + \Delta t \frac{1}{2} F^{n + 1,1} \right),
+  \end{equation}
+  \begin{equation}
+    T^{n + 1} = T^n + \Delta t F^{n + 1,2}.
+  \end{equation}
+\end{subequations}
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+from fractions import Fraction
+hf = Fraction(1, 2)
+
+levels   = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+levels_F = TimeLevels(levels = [n, n + hf], cycle_map = {},
+  last_past_level = n - hf)
+T    = TimeFunction(levels,   space, name = "T")
+F_s  = TimeFunction(levels_F, space, name = "F_s")
+
+system = TimeSystem()
+
+system.add_solve(T_ic, T[0])
+
+system.add_solve(inner(test, F_s[n]) * dx == F(test, T[n]),
+  F_s[n], solver_parameters = solver_parameters)
+system.add_solve(inner(test, F_s[n + hf]) * dx ==
+  F(test, T[n] + 0.5 * dt * F_s[n]),
+  F_s[n + hf], solver_parameters = solver_parameters)
+  
+system.add_solve(LinearCombination((1.0, T[n]),
+                                   (dt, F_s[n + hf])),
+  T[n + 1])
+\end{lstlisting}
+
+\section{Discontinuous Galerkin Finite Element}
+
+Let a time interval $(0, T)$ be divided into a series of $N$ elements
+$( (n - 1) \Delta t, n \Delta t )$,
+$n \in \left\{ 1, \ldots, N \right\}$. Equip this one dimensional mesh with
+discontinuous Lagrange basis functions, thus defining a discrete function
+space $V^\delta \in L^2$ on $(0, T)$. A temporal discontinuous Galerkin
+discretisation of equation \eqref{eqn:t_timestep} yields:
+\begin{align}\label{eqn:t_timestep_dg}
+  \left[ \phi^\delta T^{\delta,-} \right]_{(n - 1) \Delta t}^{n \Delta t}
+    - \int_{(n - 1) \Delta t}^{n \Delta t} d_t \phi^\delta T^\delta dt
+    = \int_{(n - 1) \Delta t}^{n \Delta t} \phi^\delta F \left( T^\delta \right) dt \nonumber \\
+    \forall \phi^\delta \in V^\delta, \quad \forall n \in \left\{ 1, \ldots, N \right\},
+\end{align}
+where $T^\delta \in V^\delta$ and
+$T^{\delta,-} = \lim_{t \rightarrow \left[ (n - 1) \Delta t \right]^{-}} T^{\delta}$.
+The time derivative has been integrated by parts and an ``upwind flux''
+approximation applied.
+
+\subsection{$P0$}
+
+Consider discontinuous Lagrange elements of degree zero. Introduce basis
+functions:
+\begin{equation}
+  \phi^n = \left\{ \begin{array}{l} 1 \quad \textrm{ if } t \in \left( (n - 1) \Delta t, n \Delta t \right) \\
+                                    0 \quad \textrm{ otherwise}\end{array} \right. ,
+\end{equation}
+and let $T^\delta = \sum_{n = 1}^N T^n \phi^n$. Then equation
+\eqref{eqn:t_timestep_dg} leads to:
+\begin{equation}
+  T^{n + 1} - T^n
+    = \int_{n \Delta t}^{(n + 1) \Delta t} \phi^{n + 1} F \left( T^\delta \right) dt \quad \forall n \in \left\{ 0, \ldots, N - 1 \right\}.
+\end{equation}
+Application of the midpoint rule to the right-hand-side yields:
+\begin{equation}
+  T^{n + 1} - T^n = \Delta t F \left( T^{n + 1} \right).
+\end{equation}
+This is the backward Euler scheme \ref{sect:be}, and has first order accuracy
+in time.
+
+\subsection{$P1_{DG}$}
+
+Consider discontinuous Lagrange elements of degree one. Introduce basis
+functions:
+\begin{subequations}
+  \begin{align}
+    \phi^{n,1} & = \left\{ \begin{array}{l} 1 - \frac{n \Delta t - t}{\Delta t} \\ 0 \end{array}\begin{array}{l} \textrm{ if } t \in \left( (n - 1) \Delta t, n \Delta t \right) \\ \textrm{ otherwise} \end{array} \right. , \\
+    \phi^{n,2} & = \left\{ \begin{array}{l} \frac{t - (n - 1) \Delta t}{\Delta t} \\ 0 \end{array}\begin{array}{l} \textrm{ if } t \in \left( (n - 1) \Delta t, n \Delta t \right) \\ \textrm{ otherwise} \end{array} \right. ,
+  \end{align}
+\end{subequations}
+and let $T^\delta = \sum_{n = 1}^N T^{n,1} \phi^{n,1} + T^{n,2} \phi^{n,2}$.
+Then equation \eqref{eqn:t_timestep_dg} leads to:
+\begin{subequations}
+  \begin{align}
+    -T^{n,2} + \frac{1}{2} \left( T^{n + 1,1} + T^{n + 1,2} \right)
+      = \int_{n \Delta t}^{(n + 1) \Delta t} \phi^{n + 1,1} F \left( T^\delta \right) dt \nonumber \\ \forall n \in \left\{ 0, \ldots, N - 1 \right\},
+  \end{align}
+  \begin{align}
+    T^{n + 1,2} - \frac{1}{2} \left( T^{n + 1,1} + T^{n + 1,2} \right)
+      = \int_{n \Delta t}^{(n + 1) \Delta t} \phi^{n + 1,2} F \left( T^\delta \right) dt \nonumber \\ \forall n \in \left\{ 0, \ldots, N - 1 \right\},
+  \end{align}
+\end{subequations}
+The right-hand-sides can be integrated to second order accuracy in time via
+application of the trapezium rule to yield:
+\begin{subequations}
+  \begin{align}
+    -T^{n,2} + \frac{1}{2} \left( T^{n + 1,1} + T^{n + 1,2} \right)
+      = \Delta t \left[ \frac{1}{3} F \left( T^{n + 1,1} \right) + \frac{1}{6} F \left( T^{n + 1,2} \right) \right] \nonumber \\ \forall n \in \left\{ 0, \ldots, N - 1 \right\},
+  \end{align}
+  \begin{align}
+    T^{n + 1,2} - \frac{1}{2} \left( T^{n + 1,1} + T^{n + 1,2} \right)
+      = \Delta t \left[ \frac{1}{6} F \left( T^{n + 1,1} \right) + \frac{1}{3} F \left( T^{n + 1,2} \right) \right]  \nonumber \\ \forall n \in \left\{ 0, \ldots, N - 1 \right\},
+  \end{align}
+\end{subequations}
+The complete discretisation, including the initial condition, becomes:
+\begin{subequations}
+  \begin{align}
+    \frac{1}{2} T^{n + 1,2} + \frac{1}{2} T^{n + 1,1} - T^{n,2}
+      = \Delta t \left[ \frac{1}{3} F \left( T^{n + 1,1} \right) + \frac{1}{6} F \left( T^{n + 1,2} \right) \right] \nonumber \\ \forall n \in \left\{ 0, \ldots, N - 1 \right\},
+  \end{align}
+  \begin{align}
+    \frac{1}{2} T^{n + 1,2} - \frac{1}{2} T^{n + 1,1}
+      = \Delta t \left[ \frac{1}{6} F \left( T^{n + 1,1} \right) + \frac{1}{3} F \left( T^{n + 1,2} \right) \right]  \nonumber \\ \forall n \in \left\{ 0, \ldots, N - 1 \right\},
+  \end{align}
+  \begin{equation}
+    T^{0,2} = T_0.
+  \end{equation}
+\end{subequations}
+Each timestep involves the implicit solution for two new values, $T^{n + 1,1}$
+and $T^{n + 1,2}$, given one previous value, $T^{n,2}$.
+
+This discretisation can be implemented via:
+\begin{lstlisting}
+spaces = space * space
+tests = TestFunction(spaces)
+test1, test2 = split(tests)
+
+levels = TimeLevels(levels = [n, n + 1], cycle_map = {n:n + 1})
+T = TimeFunction(levels, spaces, name = "T")
+
+system = TimeSystem()
+
+system.add_solve(inner(tests, T[0]) * dx == inner(test2, T_ic) * dx,
+  T[0], solver_parameters = solver_parameters)
+
+m_11 = m_22 = 1.0 / 3.0
+m_12 = m_21 = 1.0 / 6.0
+system.add_solve(
+  inner(test1, 0.5 * T[n + 1][1] + 0.5 * T[n + 1][0] - T[n][1])*dx +
+  inner(test2, 0.5 * T[n + 1][1] - 0.5 * T[n + 1][0]) * dx ==
+  dt*(m_11 * F(test1, T[n + 1][0]) + m_12 * F(test1, T[n + 1][1]) +
+      m_21 * F(test2, T[n + 1][0]) + m_22 * F(test2, T[n + 1][1])),
+  T[n + 1], solver_parameters = {"linear_solver":"lu"})
+\end{lstlisting}
+Note that \verb+solver_parameters+ should define Newton solver parameters if
+$F(T)$ is non-linear.
+
+\bibliography{bibliography}
+\bibliographystyle{plainnat}
+
+\end{document}

=== added directory 'timestepping/python'
=== added directory 'timestepping/python/dolfin_adjoint_timestepping'
=== added file 'timestepping/python/dolfin_adjoint_timestepping/__init__.py'
--- timestepping/python/dolfin_adjoint_timestepping/__init__.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/dolfin_adjoint_timestepping/__init__.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,51 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import dolfin_adjoint
+import dolfin_adjoint_timestepping
+
+import timestepping
+
+from timestepping import *
+from dolfin_adjoint_timestepping import *
+
+__doc__ = \
+"""
+A timestepping abstraction and automated adjoining library. This library
+utilises the FEniCS system for symbolic manipulation and automated code
+generation, and supplements this system with a syntax for the description of
+timestepping finite element models.
+
+This version of the library integrates with dolfin-adjoint.
+"""
+
+__license__ = "LGPL-3"
+
+__version__ = "1.2.0"
+
+__all__ = \
+  dolfin_adjoint_timestepping.__all__ + \
+  [
+    "__doc__",
+    "__license__",
+    "__name__",
+    "__version__",
+    "dolfin_adjoint",
+    "dolfin_adjoint_timestepping",
+    "timestepping"
+  ]
+  
\ No newline at end of file

=== added file 'timestepping/python/dolfin_adjoint_timestepping/dolfin_adjoint_timestepping.py'
--- timestepping/python/dolfin_adjoint_timestepping/dolfin_adjoint_timestepping.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/dolfin_adjoint_timestepping/dolfin_adjoint_timestepping.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,439 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+# Based on code from dolfin-adjoint bzr trunk 640
+# Code first added: 13/05/05
+
+import copy
+import types
+  
+import dolfin
+import dolfin_adjoint
+import libadjoint
+import ufl
+
+import timestepping
+
+__all__ = timestepping.__all__ + dir(dolfin_adjoint)
+for val in copy.copy(__all__):
+  if val.startswith("__") and val.endswith("__"):
+    __all__.remove(val)
+
+from timestepping import *
+from dolfin_adjoint import *
+
+def system_info():
+  """
+  Print system information and assorted library versions.
+  """
+  
+  timestepping.system_info()
+  dolfin.info("dolfin-adjoint version: %s" % dolfin_adjoint.__version__)
+  
+  return
+
+dolfin.parameters["timestepping"]["pre_assembly"]["linear_forms"]["whole_form_optimisation"] = True
+
+# dolfin-adjoint internals expect Constant s and Function s to have an
+# adj_name attribute. Patch __getattr__ to provide it.
+def adj_name__getattr__(self, name):
+  if name == "adj_name":
+    return self.name()
+  else:
+    return object.__getattr__(self, name)
+dolfin.Constant.__getattr__ = adj_name__getattr__
+dolfin.Function.__getattr__ = adj_name__getattr__
+del(adj_name__getattr__)
+
+# constant.get_constant expects to be dealing with dolfin-adjoint Constant s.
+# Patch it so that it can handle any Constant. Modified version of get_constant
+# from constant.py.
+def get_constant(a):
+  if isinstance(a, dolfin.Constant):
+    return a
+  else:
+    return constant_objects[a]
+constant.get_constant.func_code = get_constant.func_code
+del(get_constant)
+
+# Resolve namespace clashes
+def assemble(*args, **kwargs):
+  if isinstance(args[0], QForm):
+    if "form_compiler_parameters" in kwargs:
+      raise InvalidArgumentException("Cannot supply form_compiler_parameters argument when assembling a QForm")
+    return dolfin_adjoint.assemble(form_compiler_parameters = args[0].form_compiler_parameters(), *args, **kwargs)
+  elif isinstance(args[0], tuple(_assemble_classes)):
+    return args[0].assemble(*args[1:], **kwargs)
+  else:
+    return dolfin_adjoint.assemble(*args, **kwargs)
+
+def Constant(value, cell = None, name = "u"):
+  if isinstance(value, tuple):
+    return dolfin.as_vector([dolfin_adjoint.Constant(val, cell = cell, name = "%s_%i" % (name, i)) for i, val in enumerate(value)])
+  else:
+    return dolfin_adjoint.Constant(value, cell = cell, name = name)
+
+class TimeFunction(timestepping.TimeFunction):
+  def __init__(self, tlevels, space, name = "u"):
+    timestepping.TimeFunction.__init__(self, tlevels, space, name = name)
+    # Ensure that all time level Function s are defined at all times
+    cycle_map = self.final_cycle_map()
+    for level in cycle_map:
+      self.__lfns[level].wrap(self.__fns[cycle_map[level]])
+    return
+  
+  def cycle(self, extended = True):
+    """
+    Perform a timestep cycle. The optional extended argument has no effect.
+    """
+    
+    cycle_map = self.cycle_map()
+    for level in cycle_map:
+      # Annotate the timestep variable cycle
+      record = da_annotate_assign(self[cycle_map[level]], self[level])
+      assert(not record)
+      self[level].vector()[:] = self[cycle_map[level]].vector()  
+      
+    return
+
+  def final_cycle(self):
+    """
+    Perform the final cycle.
+    """
+    
+    return
+
+__ForwardModel_timestep_cycle_orig = timestepping.ForwardModel.timestep_cycle
+def ForwardModel_timestep_cycle(self, extended = True):
+  """
+  Perform the timestep cycle. The optional extended argument has no effect.
+  """
+    
+  __ForwardModel_timestep_cycle_orig(self, extended = False)
+  # Signal the end of the timestep
+  adj_inc_timestep()
+  
+  return
+timestepping.ForwardModel.timestep_cycle = ForwardModel_timestep_cycle
+del(ForwardModel_timestep_cycle)
+
+__ForwardModel_finalise_orig = timestepping.ForwardModel.finalise
+def ForwardModel_finalise(self):
+  """
+  Solve final equations and perform the final variable cycle.
+  """
+    
+  __ForwardModel_finalise_orig(self)
+  # Signal the end of the forward model
+  adj_inc_timestep()
+  
+  return
+timestepping.ForwardModel.finalise = ForwardModel_finalise
+del(ForwardModel_finalise)
+
+__AssignmentSolver__init__orig = AssignmentSolver.__init__
+def AssignmentSolver__init__(self, *args, **kwargs):
+  __AssignmentSolver__init__orig(self, *args, **kwargs)
+  solve_orig = self.solve.im_class.solve
+  # Equation annotation based on solve in solving.py
+  def solve(self, *args, **kwargs):
+    # Annotate an assignment solve
+    record = da_annotate_equation_solve(self)
+    annotate = not dolfin.parameters["adjoint"]["stop_annotating"]
+    dolfin.parameters["adjoint"]["stop_annotating"] = True
+    ret = solve_orig(self, *args, **kwargs)
+    dolfin.parameters["adjoint"]["stop_annotating"] = not annotate
+    if record:
+      x = self.x()
+      if isinstance(x, WrappedFunction):
+        x = x.fn()
+      adjglobals.adjointer.record_variable(adjglobals.adj_variables[x], libadjoint.MemoryStorage(adjlinalg.Vector(x)))
+    return ret
+  solve.__doc__ = solve_orig.__doc__
+  self.solve = types.MethodType(solve, self)
+  return
+AssignmentSolver.__init__ = AssignmentSolver__init__
+del(AssignmentSolver__init__)
+
+__EquationSolver__init__orig = EquationSolver.__init__
+def EquationSolver__init__(self, *args, **kwargs):
+  __EquationSolver__init__orig(self, *args, **kwargs)
+  solve_orig = self.solve.im_class.solve
+  # Equation annotation based on solve in solving.py
+  def solve(self, *args, **kwargs):
+    # Annotate an equation solve
+    record = da_annotate_equation_solve(self)
+    # The EquationSolver solve method could do a lot of work, including
+    # further solves. Temporarily disable annotation during the solve ...
+    annotate = not dolfin.parameters["adjoint"]["stop_annotating"]
+    dolfin.parameters["adjoint"]["stop_annotating"] = True
+    ret = solve_orig(self, *args, **kwargs)
+    # ... and then restore the previous annotation status.
+    dolfin.parameters["adjoint"]["stop_annotating"] = not annotate
+    if record:
+      x = self.x()
+      # We still want to wrap functions with WrappedFunction so that we can for
+      # example write separate equations for u[0] and u[n]. However
+      # dolfin-adjoint needs to treat these as the same Function, so
+      # strategically unwrap WrappedFunction s.
+      if isinstance(x, WrappedFunction):
+        x = x.fn()
+      adjglobals.adjointer.record_variable(adjglobals.adj_variables[x], libadjoint.MemoryStorage(adjlinalg.Vector(x)))
+    return ret
+  solve.__doc__ = solve_orig.__doc__
+  self.solve = types.MethodType(solve, self)
+  return
+EquationSolver.__init__ = EquationSolver__init__
+del(EquationSolver__init__)
+
+def unwrap_fns(form):
+  """
+  Return a form with all WrappedFunction s unwrapped.
+  """
+  
+  repl = {}
+  for dep in ufl.algorithms.extract_coefficients(form):
+    if isinstance(dep, WrappedFunction):
+      repl[dep] = dep.fn()
+      
+  return replace(form, repl)  
+
+# Based on the Functional class in functional.py
+class Functional(functional.Functional):    
+  def __init__(self, timeform, verbose = False, name = None):
+    # Unwrap WrappedFunction s in the Functional
+    if isinstance(timeform, ufl.form.Form):
+      timeform = unwrap_fns(timeform)
+    else:
+      for term in timeform.terms:
+        term.form = unwrap_fns(term.form)
+    functional.Functional.__init__(self, timeform, verbose = verbose, name = name)
+    
+    return
+  
+# Based on the ReducedFunctional class in reduced_functional.py
+class ReducedFunctional(reduced_functional.ReducedFunctional):
+  def eval_array(self, m_array):
+    # Clear caches
+    clear_caches()
+    return reduced_functional.ReducedFunctional.eval_array(self, m_array)
+
+# Based on dolfin_adjoint_assign in function.py.
+def da_annotate_assign(y, x):
+  """
+  Annotate an assignment. Returns whether the variable being solved for should
+  be recorded by dolfin-adjoint.
+  """
+  
+  if dolfin.parameters["adjoint"]["stop_annotating"]:
+    # Annotation disabled
+    return False
+  
+  if isinstance(x, WrappedFunction):
+    x = x.fn()
+  if isinstance(y, WrappedFunction):
+    y = y.fn()
+  if not x == y:
+    # ?? What does this do ??
+    if not adjglobals.adjointer.variable_known(adjglobals.adj_variables[x]):
+      adjglobals.adj_variables.forget(x)
+    assign.register_assign(x, y)
+    
+  return False
+  
+# A simple cache for use by dolfin-adjoint callbacks.
+da_matrix_cache = {}
+def clear_caches(*args):
+  """
+  Clear caches. Constant s or Function s can be supplied, indicating that only
+  cached data associated with those coefficients should be cleared.
+  """
+  
+  timestepping.clear_caches(*args)
+  da_matrix_cache.clear()
+  
+  return
+
+# Based on annotate in solving.py
+def da_annotate_equation_solve(solve):
+  """
+  Annotate an equation solve. Returns whether the variable being solved for
+  should be recorded by dolfin-adjoint.
+  """
+  
+  if dolfin.parameters["adjoint"]["stop_annotating"]:
+    # Annotation disabled
+    return False
+  
+  # The Function being solved for
+  x = solve.x()
+  if isinstance(x, WrappedFunction):
+    x_fn = x.fn()
+  else:
+    x_fn = x
+    
+  if isinstance(solve, AssignmentSolver):
+    # Assignment solve case
+    
+    rhs = solve.rhs()
+    if isinstance(rhs, (list, tuple)):
+      if len(rhs) == 1 and rhs[0][0] == 1.0:
+        # This is a direct assignment, so register an assignment
+        return da_annotate_assign(rhs[0][1], x_fn)
+      nrhs = rhs[0][0] * rhs[0][1]
+      for term in rhs[1:]:
+        nrhs += term[0] * term[1]
+      rhs = nrhs;  del(nrhs)
+    if isinstance(rhs, (float, int, ufl.constantvalue.FloatValue, ufl.constantvalue.IntValue, ufl.constantvalue.Zero)):
+      # This is a direct assignment, so register an assignment
+      return da_annotate_assign(Constant(rhs), x_fn)
+    elif isinstance(rhs, (dolfin.Constant, dolfin.Function)):
+      # This is a direct assignment, so register an assignment
+      return da_annotate_assign(rhs, x_fn)
+    # This is a LinearCombination assignment or expression assignment. For now
+    # register this as a Galerkin projection.
+    eq = dolfin.inner(dolfin.TestFunction(x.function_space()), dolfin.TrialFunction(x.function_space())) * dolfin.dx == \
+      dolfin.inner(dolfin.TestFunction(x.function_space()), rhs) * dolfin.dx
+    bcs = []
+    solver_parameters = {"linear_solver":"lu"}
+    adjoint_solver_parameters = solver_parameters
+  else:
+    # Equation solve case
+    assert(isinstance(solve, EquationSolver))
+    
+    eq = solve.eq()
+    bcs = solve.bcs()
+    solver_parameters = solve.solver_parameters()
+    adjoint_solver_parameters = solve.adjoint_solver_parameters()
+  
+  # Unwrap WrappedFunction s in the equation
+  eq.lhs = unwrap_fns(eq.lhs)
+  if not is_zero_rhs(eq.rhs):
+    eq.rhs = unwrap_fns(eq.rhs)
+  
+  if hasattr(x, "_time_level_data"):
+    # This is a time level solve. Set up a Matrix class with some caching
+    # enabled.
+    
+    class DAMatrix(adjlinalg.Matrix):
+      def __init__(self, *args, **kwargs):
+        adjlinalg.Matrix.__init__(self, *args, **kwargs)
+        self.__x = x
+        self.__x_fn = x_fn
+        self.__eq = eq
+        self.__bcs = bcs
+        self.__solver_parameters = solver_parameters
+        self.__adjoint_solver_parameters = adjoint_solver_parameters
+        self.parameters = dolfin.Parameters(**dolfin.parameters["timestepping"]["pre_assembly"])
+        
+        return
+
+      # Based on Matrix.solve in adjlinalg.py
+      def solve(self, var, b):
+        if isinstance(self.data, adjlinalg.IdentityMatrix) or b.data is None:
+          return adjlinalg.Matrix.solve(self, var, b)
+        
+        # Configure the cache
+        if not self.__x in da_matrix_cache:
+          da_matrix_cache[self.__x] = {}
+        cache = da_matrix_cache[self.__x]
+      
+        # Boundary conditions
+        if var.type in ["ADJ_ADJOINT", "ADJ_SOA", "ADJ_TLM"]:
+          bcs = [homogenize(bc) for bc in self.__bcs]
+        else:
+          bcs = self.__bcs
+        static_bcs = n_non_static_bcs(bcs) == 0
+        
+        if ("pa_a", var.type) in cache:
+          # We have cached data for this matrix
+          if cache[("pa_a", var.type)] is None:
+            # The cache is empty
+            static_a = False
+            a = assemble(self.data)
+            apply_a_bcs = True
+          else:              
+            # The cache contains a matrix, let's use it
+            static_a = True
+            a, apply_a_bcs = cache[("pa_a", var.type)]
+        else:
+          if extract_form_data(self.__eq.lhs).rank == 2:
+            assert(not self.__x_fn in ufl.algorithms.extract_coefficients(self.__eq.lhs))
+            if not is_zero_rhs(self.__eq.rhs):
+              assert(not self.__x_fn in ufl.algorithms.extract_coefficients(self.__eq.rhs))
+            # The forward equation is a linear variational problem. Is the
+            # forward LHS matrix static?
+            static_a = is_static_form(self.__eq.lhs)
+          else:
+            # Is the matrix static?
+            static_a = is_static_form(self.data)
+          if static_a:
+            # The matrix is static, so we can cache it
+            if not self.parameters["equations"]["symmetric_boundary_conditions"] and static_bcs:
+              # Cache with boundary conditions
+              a = assembly_cache.assemble(self.data, bcs = bcs)
+              apply_a_bcs = False
+            else:
+              # Cache without boundary conditions
+              a = assembly_cache.assemble(self.data)
+              apply_a_bcs = True
+            # Cache
+            cache[("pa_a", var.type)] = a, apply_a_bcs
+          else:
+            # The matrix is not static, so we cannot cache it. Add an empty
+            # entry to the cache to prevent repeated processing.
+            cache[("pa_a", var.type)] = None
+            a = assemble(self.data)
+            apply_a_bcs = True
+        
+        # Assemble the RHS
+        if isinstance(b.data, ufl.form.Form):
+          b = assemble(b.data)
+        else:
+          b = b.data.vector()
+
+        # Apply boundary conditions
+        if apply_a_bcs:
+          apply_bcs(a, bcs, L = b, symmetric_bcs = self.parameters["equations"]["symmetric_boundary_conditions"])
+        else:
+          enforce_bcs(b, bcs)
+
+        if ("solver", var.type) in cache:
+          # Extract a linear solver from the cache
+          solver = cache[("solver", var.type)]
+        else:
+          # Create a new linear solver and cache it
+          if var.type in ["ADJ_ADJOINT", "ADJ_SOA"]:
+            solver_parameters = self.__adjoint_solver_parameters
+          else:
+            solver_parameters = self.__solver_parameters
+          solver = cache[("solver", var.type)] = solver_cache.solver(self.data, solver_parameters, static = static_a and static_bcs, bcs = bcs, symmetric_bcs = self.parameters["equations"]["symmetric_boundary_conditions"])
+
+        # RHS solution vector
+        x = adjlinalg.Vector(dolfin.Function(self.__x.function_space()))
+
+        # Solve and return
+        solver.set_operator(a)
+        solver.solve(x.data.vector(), b)
+        return x
+  else:
+    # This is not a time level solve. Use the default Matrix type.
+    DAMatrix = adjlinalg.Matrix
+  
+  # Annotate the equation
+  solving.annotate(eq, x_fn, bcs, solver_parameters = solver_parameters, matrix_class = DAMatrix)
+  return True

=== added directory 'timestepping/python/timestepping'
=== added file 'timestepping/python/timestepping/__init__.py'
--- timestepping/python/timestepping/__init__.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/timestepping/__init__.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,112 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from fenics_patches import *
+from parameters import *
+
+import caches
+import checkpointing
+import equation_solvers
+import embedded_cpp
+import exceptions
+import fenics_overrides
+import fenics_utils
+import quadrature
+import pre_assembled_adjoint
+import pre_assembled_equations
+import pre_assembled_forms
+import statics
+import time_functions
+import time_levels
+import time_systems
+import versions
+import vtu_io
+
+from caches import *
+from checkpointing import *
+from equation_solvers import *
+from embedded_cpp import *
+from exceptions import *
+from fenics_overrides import *
+from fenics_utils import *
+from pre_assembled_adjoint import *
+from pre_assembled_equations import *
+from pre_assembled_forms import *
+from quadrature import *
+from statics import *
+from time_functions import *
+from time_levels import *
+from time_systems import *
+from versions import *
+from vtu_io import *
+
+__doc__ = \
+"""
+A timestepping abstraction and automated adjoining library. This library
+utilises the FEniCS system for symbolic manipulation and automated code
+generation, and supplements this system with a syntax for the description of
+timestepping finite element models.
+"""
+
+__license__ = "LGPL-3"
+
+__version__ = "1.2.0"
+
+__all__ = \
+  caches.__all__ + \
+  checkpointing.__all__ + \
+  equation_solvers.__all__ + \
+  embedded_cpp.__all__ + \
+  exceptions.__all__ + \
+  fenics_overrides.__all__ + \
+  fenics_patches.__all__ + \
+  fenics_utils.__all__ + \
+  parameters.__all__ + \
+  pre_assembled_adjoint.__all__ + \
+  pre_assembled_equations.__all__ + \
+  pre_assembled_forms.__all__ + \
+  quadrature.__all__ + \
+  statics.__all__ + \
+  time_functions.__all__ + \
+  time_levels.__all__ + \
+  time_systems.__all__ + \
+  versions.__all__ + \
+  vtu_io.__all__ + \
+  [ 
+    "__doc__",
+    "__license__",
+    "__name__",
+    "__version__",
+    "caches",
+    "checkpointing",
+    "equation_solvers",
+    "embedded_cpp",
+    "exceptions",
+    "fenics_overrides",
+    "fenics_utils",
+    "pre_assembled_adjoint",
+    "pre_assembled_equations",
+    "pre_assembled_forms",
+    "quadrature",
+    "statics",
+    "time_functions",
+    "time_levels",
+    "time_systems",
+    "versions",
+    "vtu_io"
+  ]
+  
\ No newline at end of file

=== added file 'timestepping/python/timestepping/caches.py'
--- timestepping/python/timestepping/caches.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/timestepping/caches.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,266 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from collections import OrderedDict
+import copy
+
+import dolfin
+import ufl
+
+from exceptions import *
+from fenics_overrides import *
+from fenics_utils import *
+
+__all__ = \
+  [
+    "AssemblyCache",
+    "SolverCache",
+    "assembly_cache",
+    "solver_cache",
+    "cache_info",
+    "clear_caches"
+  ]
+
+def cache_info(msg, info = dolfin.info):
+  """
+  Print a message if verbose pre-assembly is enabled.
+  """
+  
+  if dolfin.parameters["timestepping"]["pre_assembly"]["verbose"]:
+    info(msg)
+  return
+    
+class AssemblyCache:
+  """
+  A cache of assembled Form s. The assemble method can be used to assemble a
+  given Form. If an assembled version of the Form exists in the cache, then the
+  cached result is returned. Note that this does not check that the Form
+  dependencies are unchanged between subsequent assemble calls -- that is
+  deemed the responsibility of the caller.
+  """
+  
+  def __init__(self):
+    self.__cache = {}
+
+    return
+
+  def assemble(self, form, bcs = [], symmetric_bcs = False):
+    """
+    Return the result of assembling the supplied Form, optionally with boundary
+    conditions, which are optionally applied so as to yield a symmetric matrix.
+    If an assembled version of the Form exists in the cache, return the cached
+    result. Note that this does not check that the Form dependencies are
+    unchanged between subsequent assemble calls -- that is deemed the
+    responsibility of the caller.
+    """
+    
+    if not isinstance(form, ufl.form.Form):
+      raise InvalidArgumentException("form must be a Form")
+    if not isinstance(bcs, list):
+      raise InvalidArgumentException("bcs must be a list of DirichletBC s")
+    for bc in bcs:
+      if not isinstance(bc, dolfin.cpp.DirichletBC):
+        raise InvalidArgumentException("bcs must be a list of DirichletBC s")
+
+    form_data = extract_form_data(form)
+    if len(bcs) == 0:
+      key = (expand(form), None, None)
+      if not key in self.__cache:
+        cache_info("Assembling form with rank %i" % form_data.rank, dolfin.info_red)
+        self.__cache[key] = assemble(form)
+      else:
+        cache_info("Using cached assembled form with rank %i" % form_data.rank, dolfin.info_green)
+    else:
+      if not form_data.rank == 2:
+        raise InvalidArgumentException("form must be rank 2 when applying boundary conditions")
+
+      key = (expand(form), tuple(bcs), symmetric_bcs)
+      if not key in self.__cache:
+        cache_info("Assembling form with rank 2, with boundary conditions", dolfin.info_red)
+        mat = assemble(form)
+        apply_bcs(mat, bcs, symmetric_bcs = symmetric_bcs)
+        self.__cache[key] = mat
+      else:
+        cache_info("Using cached assembled form with rank 2, with boundary conditions", dolfin.info_green)
+
+    return self.__cache[key]
+
+  def info(self):
+    """
+    Print some cache status information.
+    """
+    
+    counts = [0, 0, 0]
+    for key in self.__cache.keys():
+      counts[extract_form_data(key[0]).rank] += 1
+
+    dolfin.info_blue("Assembly cache status:")
+    for i in range(3):
+      dolfin.info_blue("Pre-assembled rank %i forms: %i" % (i, counts[i]))
+
+    return
+
+  def clear(self, *args):
+    """
+    Clear the cache. If arguments are supplied, clear only the cached assembled
+    Form s which depend upon the supplied Constant s or Function s.
+    """
+    
+    if len(args) == 0:
+      self.__cache = {}
+    else:
+      for dep in args:
+        if not isinstance(dep, (dolfin.Constant, dolfin.Function)):
+          raise InvalidArgumentException("Arguments must be Constant s or Function s")
+
+      for dep in args:
+        for key in copy.copy(self.__cache.keys()):
+          form = key[0]
+          if dep in ufl.algorithms.extract_coefficients(form):
+            del(self.__cache[key])
+
+    return
+
+class SolverCache:
+  """
+  A cache of LUSolver s and KrylovSolver s. The solver method can be used to
+  return an LUSolver or KrylovSolver suitable for solving an equation with the
+  supplied rank 2 Form defining the LHS matrix.
+  """
+  
+  def __init__(self):
+    self.__cache = OrderedDict()
+
+    return
+
+  def __del__(self):
+    for key in self.__cache:
+      del(self.__cache[key])
+
+    return
+
+  def solver(self, form, solver_parameters, static = False, bcs = [], symmetric_bcs = False):
+    """
+    Return an LUSolver or KrylovSolver suitable for solving an equation with the
+    supplied rank 2 Form defining the LHS. If such a solver exists in the cache,
+    return the cached solver. Optionally accepts boundary conditions which
+    are optionally applied so as to yield a symmetric matrix. If static is true
+    then it is assumed that the Form is "static", and solver caching
+    options are enabled. The appropriate value of the static argument is
+    deemed the responsibility of the caller.
+    """
+    
+    if not isinstance(form, ufl.form.Form):
+      raise InvalidArgumentException("form must be a rank 2 Form")
+    elif not extract_form_data(form).rank == 2:
+      raise InvalidArgumentException("form must be a rank 2 Form")
+    if not isinstance(solver_parameters, dict):
+      raise InvalidArgumentException("solver_parameters must be a dictionary")
+    if not isinstance(bcs, list):
+      raise InvalidArgumentException("bcs must be a list of DirichletBC s")
+    for bc in bcs:
+      if not isinstance(bc, dolfin.cpp.DirichletBC):
+        raise InvalidArgumentException("bcs must be a list of DirichletBC s")
+
+    def flatten_parameters(opts):
+      assert(isinstance(opts, dict))
+      fopts = []
+      for key in opts.keys():
+        if isinstance(opts[key], dict):
+          fopts.append(flatten_parameters(opts[key]))
+        else:
+          fopts.append((key, opts[key]))
+      return tuple(fopts)
+
+    if static:
+      if not "linear_solver" in solver_parameters or solver_parameters["linear_solver"] == "lu":
+        solver_parameters = expand_solver_parameters(solver_parameters, default_solver_parameters = {"lu_solver":{"reuse_factorization":True, "same_nonzero_pattern":True}})
+      else:
+        solver_parameters = expand_solver_parameters(solver_parameters, default_solver_parameters = {"krylov_solver":{"preconditioner":{"reuse":True}}})
+    else:
+      solver_parameters = expand_solver_parameters(solver_parameters)
+
+      static_parameters = False
+      if solver_parameters["linear_solver"] == "lu":
+        static_parameters = solver_parameters["lu_solver"]["reuse_factorization"] or solver_parameters["lu_solver"]["same_nonzero_pattern"]
+      else:
+        static_parameters = solver_parameters["krylov_solver"]["preconditioner"]["reuse"]
+      if static_parameters:
+        raise ParameterException("Non-static solve supplied with static solver parameters")
+
+    if static:
+      if len(bcs) == 0:
+        key = (expand(form), flatten_parameters(solver_parameters))
+      else:
+        key = (expand(form), tuple(bcs), symmetric_bcs, flatten_parameters(solver_parameters))
+    else:
+      args = ufl.algorithms.extract_arguments(form)
+      assert(len(args) == 2)
+      test, trial = args
+      if test.count() > trial.count():
+        test, trial = trial, test
+      if len(bcs) == 0:
+        key = ((test, trial), flatten_parameters(solver_parameters))
+      else:
+        key = ((test, trial), tuple(bcs), symmetric_bcs, flatten_parameters(solver_parameters))
+
+    if not key in self.__cache:
+      if static:
+        cache_info("Creating new static linear solver", dolfin.info_red)
+      else:
+        cache_info("Creating new non-static linear solver", dolfin.info_red)
+      self.__cache[key] = LinearSolver(solver_parameters)
+    else:
+      cache_info("Using cached linear solver", dolfin.info_green)
+    return self.__cache[key]
+
+  def clear(self, *args):
+    """
+    Clear the cache. If arguments are supplied, clear only the solvers
+    associated with Form s which depend upon the supplied Constant s or
+    Function s.
+    """
+    
+    if len(args) == 0:
+      for key in self.__cache.keys():
+        del(self.__cache[key])
+    else:
+      for dep in args:
+        if not isinstance(dep, (dolfin.Constant, dolfin.Function)):
+          raise InvalidArgumentException("Arguments must be Constant s or Function s")
+
+      for key in self.__cache:
+        form = key[0]
+        if isinstance(form, ufl.form.Form) and dep in ufl.algorithms.extract_coefficients(form):
+          del(self.__cache[key])
+
+    return
+
+# Default assembly and solver caches.
+assembly_cache = AssemblyCache()
+solver_cache = SolverCache()
+def clear_caches(*args):
+  """
+  Clear the default assembly and solver caches. If arguments are supplied, clear
+  only cached data associated with Form s which depend upon the supplied
+  Constant s or Function s.
+  """
+
+  assembly_cache.clear(*args)
+  solver_cache.clear(*args)
+  
+  return
\ No newline at end of file

=== added file 'timestepping/python/timestepping/checkpointing.py'
--- timestepping/python/timestepping/checkpointing.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/timestepping/checkpointing.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,394 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import copy
+import pickle
+import os
+
+import dolfin
+
+from exceptions import *
+
+__all__ = \
+  [
+    "Checkpointer",
+    "DiskCheckpointer",
+    "MemoryCheckpointer"
+  ]
+
+class Checkpointer:
+  """
+  A template for Constant and Function storage.
+  """
+  
+  def __init__(self):
+    return
+
+  def __pack(self, c):
+    if isinstance(c, dolfin.Constant):
+      return float(c)
+    else:
+      assert(isinstance(c, dolfin.Function))
+      return c.vector().array()
+
+  def __unpack(self, c, c_c):
+    if isinstance(c, dolfin.Constant):
+      c.assign(c_c)
+    else:
+      assert(isinstance(c, dolfin.Function))
+      c.vector().set_local(c_c)
+      c.vector().apply("insert")
+
+    return
+
+  def __verify(self, c, c_c, tolerance = 0.0):
+    if isinstance(c, dolfin.Constant):
+      err = abs(float(c) - c_c)
+      if err > tolerance:
+        raise CheckpointException("Invalid checkpoint data for Constant with value %.6g, error %.6g" % (float(c), err))
+    else:
+      assert(isinstance(c, dolfin.Function))
+      err = abs(c.vector().array() - c_c).max()
+      if err > tolerance:
+        raise CheckpointException("Invalid checkpoint data for Function %s, error %.6g" % (c.name(), err))
+
+    return
+        
+  def __check_cs(self, cs):
+    if not isinstance(cs, (list, set)):
+      raise InvalidArgumentException("cs must be a list of Constant s or Function s")
+    for c in cs:
+      if not isinstance(c, (dolfin.Constant, dolfin.Function)):
+        raise InvalidArgumentException("cs must be a list of Constant s or Function s")
+
+    if not isinstance(cs, set):
+      cs = set(cs)
+    cs = list(cs)
+    def cmp(x, y):
+      return x.id() - y.id()
+    cs.sort(cmp = cmp)
+    
+    return cs
+    
+  def checkpoint(self, key, cs):
+    """
+    Store, with the supplied key, the supplied Constant s and Function s.
+    """
+    
+    raise AbstractMethodException("checkpoint method not overridden")
+
+  def restore(self, key, cs = None):
+    """
+    Restore Constant s and Function s with the given key. If cs is supplied,
+    only restore Constant s and Function s found in cs.
+    """
+    
+    raise AbstractMethodException("restore method not overridden")
+
+  def has_key(key):
+    """
+    Return whether any data is associated with the given key.
+    """
+    
+    raise AbstractMethodException("has_key method not overridden")
+    
+  def verify(self, key, tolerance = 0.0):
+    """
+    Verify data associated with the given key, with the specified tolerance.
+    """
+    
+    raise AbstractMethodException("verify method not overridden")
+
+  def remove(self, key):
+    """
+    Remove data associated with the given key.
+    """
+    
+    raise AbstractMethodException("remove method not overridden")
+
+  def clear(self, keep = []):
+    """
+    Clear all stored data, except for those with keys in keep.
+    """
+    
+    raise AbstractMethodException("clear method not overridden")
+    
+class MemoryCheckpointer(Checkpointer):
+  """
+  Constant and Function storage in memory.
+  """
+  
+  def __init__(self):
+    Checkpointer.__init__(self)
+    
+    self.__cache = {}
+
+    return
+
+  def checkpoint(self, key, cs):
+    """
+    Store, with the supplied key, the supplied Constant s and Function s.
+    """
+    
+    if key in self.__cache:
+      raise CheckpointException("Attempting to overwrite checkpoint with key %s" % str(key))
+    cs = self._Checkpointer__check_cs(cs)
+  
+    c_cs = {}
+    for c in cs:
+      c_cs[c] = self._Checkpointer__pack(c)
+
+    self.__cache[key] = c_cs
+
+    return
+
+  def restore(self, key, cs = None):
+    """
+    Restore Constant s and Function s with the given key. If cs is supplied,
+    only restore Constant s and Function s found in cs.
+    """
+    
+    if not key in self.__cache:
+      raise CheckpointException("Missing checkpoint with key %s" % str(key))
+    if not cs is None:
+      cs = self._Checkpointer__check_cs(cs)
+
+    c_cs = self.__cache[key]
+    if cs is None:
+      cs = c_cs.keys()
+      
+    for c in cs:
+      self._Checkpointer__unpack(c, c_cs[c])
+
+    return
+
+  def has_key(self, key):
+    """
+    Return whether any data is associated with the given key.
+    """
+    
+    return key in self.__cache
+
+  def verify(self, key, tolerance = 0.0):
+    """
+    Verify data associated with the given key, with the specified tolerance.
+    """
+    
+    if not key in self.__cache:
+      raise CheckpointException("Missing checkpoint with key %s" % str(key))
+    if not isinstance(tolerance, float) or tolerance < 0.0:
+      raise InvalidArgumentException("tolerance must be a non-negative float")
+    c_cs = self.__cache[key]
+
+    try:
+      for c in c_cs:
+        self._Checkpointer__verify(c, c_cs[c], tolerance = tolerance)
+      dolfin.info_green("Verified checkpoint with key %s" % str(key))
+    except CheckpointException as e:
+      dolfin.info_red(str(e))
+      raise CheckpointException("Failed to verify checkpoint with key %s" % str(key))
+
+    return
+
+  def remove(self, key):
+    """
+    Remove data associated with the given key.
+    """
+    
+    if not key in self.__cache:
+      raise CheckpointException("Missing checkpoint with key %s" % str(key))
+
+    del(self.__cache[key])
+
+    return
+
+  def clear(self, keep = []):
+    """
+    Clear all stored data, except for those with keys in keep.
+    """
+    
+    if not isinstance(keep, list):
+      raise InvalidArgumentException("keep must be a list")
+    
+    if len(keep) == 0:
+      self.__cache = {}
+    else:
+      for key in copy.copy(self.__cache.keys()):
+        if not key in keep:
+          del(self.__cache[key])
+
+    return
+
+class DiskCheckpointer(Checkpointer):
+  """
+  Constant and Function storage on disk. All keys handled by a DiskCheckpointer
+  are internally cast to strings.
+
+  Constructor arguments:
+    dirname: The directory in which data is to be stored.
+  """
+  
+  def __init__(self, dirname = "checkpoints~"):
+    if not isinstance(dirname, str):
+      raise InvalidArgumentException("dirname must be a string")
+    
+    Checkpointer.__init__(self)
+
+    if dolfin.MPI.process_number() == 0:    
+      if not os.path.exists(dirname):
+        os.mkdir(dirname)
+    dolfin.MPI.barrier()
+    
+    self.__dirname = dirname
+    self.__filenames = {}
+    self.__id_map = {}
+
+    return
+
+  def __filename(self, key):
+    return os.path.join(self.__dirname, "checkpoint_%s_%i" % (str(key), dolfin.MPI.process_number()))
+
+  def checkpoint(self, key, cs):
+    """
+    Store, with the supplied key, the supplied Constant s and Function s. The
+    key is internally cast to a string.
+    """
+    
+    key = str(key)
+    if key in self.__filenames:
+      raise CheckpointException("Attempting to overwrite checkpoint with key %s" % key)
+    cs = self._Checkpointer__check_cs(cs)
+
+    c_cs = {}
+    id_map = {}
+    for c in cs:
+      c_id = c.id()
+      c_cs[c_id] = self._Checkpointer__pack(c)
+      id_map[c_id] = c
+
+    filename = self.__filename(key)
+    handle = open(filename, "wb")
+    pickler = pickle.Pickler(handle, -1)
+    pickler.dump(c_cs)
+
+    self.__filenames[key] = filename
+    self.__id_map[key] = id_map
+
+    return
+
+  def restore(self, key, cs = None):
+    """
+    Restore Constant s and Function s with the given key. If cs is supplied,
+    only restore Constant s and Function s found in cs. The key is internally
+    cast to a string.
+    """
+    
+    key = str(key)
+    if not key in self.__filenames:
+      raise CheckpointException("Missing checkpoint with key %s" % key)
+    if not cs is None:
+      cs = self._Checkpointer__check_cs(cs)
+      cs = [c.id() for c in cs]
+
+    handle = open(self.__filename(key), "rb")
+    pickler = pickle.Unpickler(handle)
+    c_cs = pickler.load()
+    if cs is None:
+      cs = c_cs.keys()
+
+    id_map = self.__id_map[key]
+    for c_id in cs:
+      c = id_map[c_id]
+      self._Checkpointer__unpack(c, c_cs[c_id])
+
+    return
+
+  def has_key(self, key):
+    """
+    Return whether any data is associated with the given key. The key is
+    internally cast to a string.
+    """
+    
+    key = str(key)
+    return key in self.__filenames
+
+  def verify(self, key, tolerance = 0.0):
+    """
+    Verify data associated with the given key, with the specified tolerance. The
+    key is internally cast to a string.
+    """
+    
+    key = str(key)
+    if not key in self.__filenames:
+      raise CheckpointException("Missing checkpoint with key %s" % key)
+    if not isinstance(tolerance, float) or tolerance < 0.0:
+      raise InvalidArgumentException("tolerance must be a non-negative float")
+    handle = open(self.__filename(key), "rb")
+    pickler = pickle.Unpickler(handle)
+    c_cs = pickler.load()
+
+    try:
+      id_map = self.__id_map[key]
+      for c_id in c_cs:
+        c = id_map[c_id]
+        self._Checkpointer__verify(c, c_cs[c_id], tolerance = tolerance)
+      dolfin.info_green("Verified checkpoint with key %s" % key)
+    except CheckpointException as e:
+      dolfin.info_red(str(e))
+      raise CheckpointException("Failed to verify checkpoint with key %s" % key)
+
+    return
+
+  def remove(self, key):
+    """
+    Remove data associated with the given key. The key is internally cast to a
+    string.
+    """
+    
+    key = str(key)
+    if not key in self.__filenames:
+      raise CheckpointException("Missing checkpoint with key %s" % key)
+
+#    os.remove(self.__filenames[key])
+    del(self.__filenames[key])
+    del(self.__id_map[key])
+
+    return
+
+  def clear(self, keep = []):
+    """
+    Clear all stored data, except for those with keys in keep. The keys are
+    internally cast to strings.
+    """
+    
+    if not isinstance(keep, list):
+      raise InvalidArgumentException("keep must be a list")
+
+    if len(keep) == 0:
+#      for key in self.__filenames:
+#        os.remove(self.__filenames[key])
+      self.__filenames = {}
+      self.__id_map = {}
+    else:
+      keep = [str(key) for key in keep]
+      for key in copy.copy(self.__filenames.keys()):
+        if not key in keep:
+ #         os.remove(self.__filenames[key])
+          del(self.__filenames[key])
+          del(self.__id_map[key])
+
+    return
\ No newline at end of file

=== added file 'timestepping/python/timestepping/embedded_cpp.py'
--- timestepping/python/timestepping/embedded_cpp.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/timestepping/embedded_cpp.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,232 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import copy
+import ctypes
+import os
+
+import dolfin
+import instant
+import numpy
+
+from exceptions import *
+from fenics_versions import *
+
+__all__ = \
+  [
+    "EmbeddedCpp",
+    "CellKernel",
+    "double_arr",
+    "int_arr",
+    "long_arr"
+  ]
+
+int_arr, long_arr, double_arr = 2, 3, 4
+
+class EmbeddedCpp:
+  """
+  A wrapper for short sections of embedded C++ code.
+
+  Constructor arguments:
+    code:          C++ code.
+    includes:      Code which can, for example be used to include header files.
+    include_dirs:  Header file directories.
+  Remaining keyword arguments form a list of name:type pairs, with:
+    name:          The name of a variable in the code, which will be passed from
+                   Python.
+    type:          One of int, float, int_arr, long_arr double_arr, Function,
+                   GenericVector, or Mesh, identifying the variable type.
+  """
+
+  if dolfin_version() < (1, 1, 0):
+    __default_includes = """#include "dolfin.h"
+#define la_index uint"""
+  else:
+    __default_includes = """#include "dolfin.h" """
+
+  def __init__(self, code, includes = "", include_dirs = [], **kwargs):
+    if not isinstance(code, str):
+      raise InvalidArgumentException("code must be a string")
+    if not isinstance(includes, str):
+      raise InvalidArgumentException("includes must be a string")
+    if not isinstance(include_dirs, list):
+      raise InvalidArgumentException("include_dirs must be a list of strings")
+    for path in include_dirs:
+      if not isinstance(path, str):
+        raise InvalidArgumentException("include_dirs must be a list of strings")
+    for arg in kwargs.keys():
+      if not isinstance(arg, str):
+        raise InvalidArgumentException("Argument name must be a string")
+    for arg in kwargs.values():
+      if not arg in [int, float, int_arr, double_arr, long_arr, dolfin.Function, dolfin.GenericVector, dolfin.Mesh]:
+        raise InvalidArgumentException("Argument type must be int, float, int_arr, long_arr, double_arr, Function, GenericVector or Mesh")
+
+    self.__code = code
+    self.__includes = """%s
+
+%s""" % (self.__default_includes, includes)
+    self.__include_dirs = copy.copy(include_dirs)
+    self.__args = copy.copy(kwargs)
+    self.__lib = None
+
+    return
+
+  def compile(self):
+    """
+    Compile the code.
+    """
+
+    args = ""
+    cast_code = ""
+    for name in sorted(self.__args.keys()):
+      arg = self.__args[name]
+      if len(args) > 0:
+        args += ", "
+      if arg == int:
+        args += "int %s" % name
+      elif arg == float:
+        args += "double %s" % name
+      elif arg == int_arr:
+        args += "int* %s" % name
+      elif arg == long_arr:
+        args += "long* %s" % name
+      elif arg == double_arr:
+        args += "double* %s" % name
+      else:
+        name_mangle = name
+        while name_mangle in self.__args.keys():
+          name_mangle = "%s_" % name_mangle
+        args += "void* %s" % name_mangle
+        if arg == dolfin.Function:
+          cast_code += "    boost::shared_ptr<Function> %s = (*((boost::shared_ptr<Function>*)%s));\n" % (name, name_mangle)
+        elif arg == dolfin.GenericVector:
+          cast_code += "    boost::shared_ptr<GenericVector> %s = (*((boost::shared_ptr<GenericVector>*)%s));\n" % (name, name_mangle)
+        else:
+          assert(arg == dolfin.Mesh)
+          cast_code += "    boost::shared_ptr<Mesh> %s = (*((boost::shared_ptr<Mesh>*)%s));\n" % (name, name_mangle)
+
+    code = \
+"""%s
+
+using namespace dolfin;
+
+extern "C" {
+  int code(%s){
+%s
+
+%s
+    return 0;
+  }
+}""" % (self.__includes, args, cast_code, self.__code)
+
+    mod = instant.build_module(code = code, cppargs = dolfin.parameters["form_compiler"]["cpp_optimize_flags"], include_dirs = self.__include_dirs)
+    path = os.path.dirname(mod.__file__)
+    name = os.path.split(path)[-1]
+    self.__lib = ctypes.cdll.LoadLibrary(os.path.join(path, "_%s.so" % name))
+    self.__lib.code.restype = int
+
+    return
+
+  def run(self, **kwargs):
+    """
+    Run the code. The keyword arguments form a list of name:variable pairs,
+    with:
+      name:     The name of a variable in the C++ code, which will be passed
+                from Python.
+      variable: The variable to be passed to the C++ code. The type must be
+                consistent with the type passed to the constructor.
+    """
+
+    args = kwargs
+    if not len(args) == len(self.__args) or not tuple(sorted(args.keys())) == tuple(sorted(self.__args.keys())):
+      raise InvalidArgumentException("Invalid argument names")
+    for name in args.keys():
+      arg = args[name]
+      if isinstance(arg, numpy.ndarray):
+        if arg.dtype == "int32":
+          if not self.__args[name] == int_arr:
+            raise InvalidArgumentException("Argument %s is of invalid type" % name)
+        elif arg.dtype == "int64":
+          if not self.__args[name] == long_arr:
+            raise InvalidArgumentException("Argument %s is of invalid type" % name)
+        elif arg.dtype == "float64":
+          if not self.__args[name] == double_arr:
+            raise InvalidArgumentException("Argument %s is of invalid type" % name)
+        else:
+          raise InvalidArgumentException("Argument %s is of invalid type" % name)
+      elif not isinstance(arg, self.__args[name]):
+        raise InvalidArgumentException("Argument %s is of invalid type" % name)
+
+    if self.__lib is None:
+      self.compile()
+
+    largs = []
+    for name in sorted(args.keys()):
+      arg = args[name]
+      if isinstance(arg, int):
+        largs.append(ctypes.c_int(arg))
+      elif isinstance(arg, float):
+        largs.append(ctypes.c_double(arg))
+      elif isinstance(arg, numpy.ndarray):
+        largs.append(arg.ctypes.data)
+      else:
+        assert(isinstance(arg, (dolfin.Function, dolfin.GenericVector, dolfin.Mesh)))
+        largs.append(ctypes.c_void_p(int(arg.this)))
+
+    ret = self.__lib.code(*largs)
+    if not ret == 0:
+      raise StateException("Non-zero return value: %i" % ret)
+
+    return
+
+class CellKernel(EmbeddedCpp):
+  """
+  A wrapper for short sections of embedded C++ code which iterate over cells.
+
+  Constructor arguments:
+    mesh:                 A Mesh.
+    kernel_code:          Code executed for each cell in the mesh.
+    initialisation_code:  Code executed before iterating over the cells.
+    finalisation_code:    Code executed after iterating over the cells.
+    includes:             Code which can, for example, be used to include header
+                          files.
+    include_dirs:         Header file directories.
+  Remaining keyword arguments are as for the EmbeddedCpp constructor. An
+  additional size_t variable cell is defined in the cell iteration, indicating
+  the cell number.
+  """
+
+  def __init__(self, mesh, kernel_code, initialisation_code = "", finalisation_code = "", includes = "", include_dirs = [], **kwargs):
+    if not isinstance(mesh, dolfin.Mesh):
+      raise InvalidArgumentException("mesh must be a Mesh")
+    if not isinstance(kernel_code, str):
+      raise InvalidArgumentException("kernel_code must be a string")
+    if not isinstance(initialisation_code, str):
+      raise InvalidArgumentException("initialisation_code must be a string")
+    if not isinstance(finalisation_code, str):
+      raise InvalidArgumentException("finalisation_code must be a string")
+
+    code = \
+"""%s
+    for(size_t cell = 0;cell < %i;cell++){
+%s
+    }
+%s""" % (initialisation_code, mesh.num_cells(), kernel_code, finalisation_code)
+
+    EmbeddedCpp.__init__(self, code, includes = includes, include_dirs = include_dirs, **kwargs)
+    
+    return
\ No newline at end of file

=== added file 'timestepping/python/timestepping/equation_solvers.py'
--- timestepping/python/timestepping/equation_solvers.py	1970-01-01 00:00:00 +0000
+++ timestepping/python/timestepping/equation_solvers.py	2013-05-15 14:57:23 +0000
@@ -0,0 +1,556 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2011-2012 by Imperial College London
+# Copyright (C) 2013 University of Oxford
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, version 3 of the License
+#
+# 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 Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import copy
+
+import dolfin
+import ufl
+
+from exceptions import *
+from fenics_overrides import *
+from fenics_utils import *
+
+__all__ = \
+  [
+    "AssignmentSolver",
+    "EquationSolver",
+    "LinearCombination"
+  ]
+
+class LinearCombination:
+  """
+  A linear combination.
+
+  Constructor arguments: An arbitrary number of tuples, each of which contains
+  two elements:
+    (alpha, x)
+  where alpha is one of:
+      1. A float.
+    or:
+      2. An arbitrary Expr.
+  and x is a Function. The supplied alphas are not checked to ensure that this
+  really is a true linear combination.
+  """
+  
+  def __init__(self, *args):
+    for arg in args:
+      if not isinstance(arg, tuple) or not len(arg) == 2 \
+        or not isinstance(arg[0], (float, ufl.expr.Expr)) or not isinstance(arg[1], dolfin.Function):
+        raise InvalidArgumentException("Require tuples of (float or Expr, Function) as arguments")
+
+    alpha = [arg[0] for arg in args]
+    for i, lalpha in enumerate(alpha):
+      if isinstance(lalpha, float):
+        alpha[i] = ufl.constantvalue.FloatValue(lalpha)
+    y = [arg[1] for arg in args]
+
+    use_exp = False
+    for lalpha in alpha:
+      if not isinstance(lalpha, (ufl.constantvalue.FloatValue, dolfin.Constant)):
+        use_exp = True
+        break
+    if use_exp:
+      fl = ufl.constantvalue.Zero()
+      for lalpha, ly in zip(alpha, y):
+        fl += lalpha * ly
+    else:
+      fl = [(lalpha, ly) for lalpha, ly in zip(alpha, y)]
+
+    self.__alpha = alpha
+    self.__y = y
+    self.__fl = fl
+
+    return
+
+  def dependencies(self, non_symbolic = False):
+    """
+    Return all dependencies associated with the LinearCombination. The
+    optional non_symbolic argument has no effect.
+    """
+    
+    deps = copy.copy(self.__y)
+    for alpha in self.__alpha:
+      deps += ufl.algorithms.extract_coefficients(alpha)
+
+    return deps
+
+  def nonlinear_dependencies(self):
+    """
+    Return all non-linear dependencies associated with the LinearCombination.
+    """
+    
+    if isinstance(self.__fl, ufl.expr.Expr):
+      nl_deps = []
+      for dep in self.dependencies():
+        if isinstance(dep, dolfin.Function):
+          nl_deps += ufl.algorithms.extract_coefficients(differentiate_expr(self.__fl, dep))
+        elif not isinstance(dep, (dolfin.Constant, dolfin.Expression)):
+          raise DependencyException("Invalid dependency")
+      return nl_deps
+    else:
+      return []
+
+  def flatten(self):
+    """
+    Return a flattened version of the LinearCombination. If all alphas are
+    floats or Constant s then return a list of (alpha, x) tuples. Otherwise
+    return an Expr.
+    """
+    
+    return self.__fl
+
+  def solve(self, x):
+    """
+    Assign the given Function x to be equal to the value of the
+    LinearCombination.
+    """
+