Re-organise time level handling

Store the current Cactus time (and not a fake Carpet time) in the th
"time hiearchy".  This removes the now redundant "leveltimes" data
structure in Carpet.

Add past time levels to th, so that it can store the time for past
time levels instead of assuming the time step size is constant.  This
allows changing the time step size during evolution.

Share the time hierarchy between all maps, instead of having one time
hierarchy per map.

Simplify the time level cycling and time stepping code used during
evolution.

Improve structure of the code that loops over time levels for certain
schedule bins.  Introduce a new Carpet variable "timelevel", similar
to "reflevel".

This also makes it possible to avoid time interpolation for the past
time levels during regridding.  The past time levels of the fine grid
then remain aligned (in time) with the past time levels of the coarse
grid.  This is controlled by a new parameter
"time_interpolation_during_regridding", which defaults to "yes" for
backwards compatibility.

Simplify the three time level initialisation.  Instead of initialising
all three time levels by taking altogether three time steps (forwards
and backwards), initialise only one past time level by taking one time
step backwards.  The remaining time level is initialised during the
first time step of the evolution, which begins by cycling time levels,
which drops the non-initialised last time level anyway.

Update Carpet and the mode handling correspondingly.

Update the CarpetIOHDF5 checkpoint format correspondingly.

Update CarpetInterp, CarpetReduce, and CarpetRegrid2 correspondingly.

Update CarpetJacobi and CarpetMG correspondingly.

6 files changed, 208 insertions, 89 deletions

diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index 3bfdcd585..a2b366ea2 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -354,6 +354,7 @@ template size_t memoryof (vector<vector<ibbox> > const & v);
 template size_t memoryof (vector<vector<dh::dboxes> > const & v);
 template size_t memoryof (vector<vector<dh::fast_dboxes> > const & v);
 template size_t memoryof (vector<vector<region_t> > const & v);
+template size_t memoryof (vector<vector<vector<CCTK_REAL> > > const & v);
 template size_t memoryof (vector<vector<vector<dh::dboxes> > > const & v);
 template size_t memoryof (vector<vector<vector<dh::fast_dboxes> > > const & v);
 template size_t memoryof (vector<vector<vector<region_t> > > const & v);
@@ -394,6 +395,7 @@ template ostream& output (ostream& os, const vector<CCTK_REAL>& v);
 template ostream& output (ostream& os, const vector<ibbox>& v);
 template ostream& output (ostream& os, const vector<rbbox>& v);
 template ostream& output (ostream& os, const vector<ivect>& v);
+template ostream& output (ostream& os, const vector<rvect>& v);
 template ostream& output (ostream& os, const vector<bbvect>& v);
 template ostream& output (ostream& os, const vector<i2vect>& v);
 template ostream& output (ostream& os, const vector<dh::dboxes> & v);
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 7d947fae5..8a7c643da 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -38,7 +38,7 @@ ggf::ggf (const int varindex_, const operator_type transport_operator_,
     vectorlength(vectorlength_), vectorindex(vectorindex_),
     vectorleader(vectorleader_)
 {
-  assert (&t.h == &d.h);
+  // assert (&t.h == &d.h);
   
   assert (vectorlength >= 1);
   assert (vectorindex >= 0 and vectorindex < vectorlength);
@@ -176,18 +176,6 @@ void ggf::recompose_fill (comm_state & state, int const rl,
     
     assert (d.fast_boxes.AT(ml).AT(rl).do_init);
     
-    vector <int> tls;
-    if (do_prolongate and rl > 0 and
-        transport_operator != op_none and transport_operator != op_sync and
-        transport_operator != op_restrict)
-    {
-      int const numtl = timelevels (ml, rl);
-      tls.resize (numtl);
-      for (int tl = 0; tl < numtl; ++ tl) {
-        tls.AT(tl) = tl;
-      }
-    }
-    
     // Initialise from the same level of the old hierarchy, where
     // possible
     if (rl < (int)oldstorage.AT(ml).size()) {
@@ -207,12 +195,18 @@ void ggf::recompose_fill (comm_state & state, int const rl,
         if (transport_operator != op_none and transport_operator != op_sync and
             transport_operator != op_restrict)
         {
+          int const numtl = timelevels (ml, rl);
+          vector <int> tls (numtl);
+          for (int tl = 0; tl < numtl; ++ tl) {
+            tls.AT(tl) = tl;
+          }
+
           for (int tl = 0; tl < timelevels (ml, rl); ++tl) {
             transfer_from_all (state,
                                tl, rl, ml,
                                & dh::fast_dboxes::fast_old2new_ref_prol_sendrecv,
                                tls, rl - 1, ml,
-                               t.time (tl, rl, ml));
+                               t.get_time (ml, rl, tl));
           } // for tl
         } // if transport_operator
       } // if rl
@@ -277,6 +271,22 @@ void ggf::cycle_all (int const rl, int const ml) {
   }
 }
 
+// Uncycle the time levels by rotating the data sets
+void ggf::uncycle_all (int const rl, int const ml) {
+  assert (rl>=0 and rl<h.reflevels());
+  assert (ml>=0 and ml<h.mglevels());
+  int const ntl = timelevels(ml,rl);
+  assert (ntl > 0);
+  for (int lc=0; lc<(int)storage.AT(ml).AT(rl).size(); ++lc) {
+    fdata & fdatas = storage.AT(ml).AT(rl).AT(lc);
+    gdata * const tmpdata = fdatas.AT(0);
+    for (int tl=0; tl<ntl-1; ++tl) {
+      fdatas.AT(tl) = fdatas.AT(tl+1);
+    }
+    fdatas.AT(ntl-1) = tmpdata;
+  }
+}
+
 // Flip the time levels by exchanging the data sets
 void ggf::flip_all (int const rl, int const ml) {
   assert (rl>=0 and rl<h.reflevels());
@@ -285,9 +295,9 @@ void ggf::flip_all (int const rl, int const ml) {
   assert (ntl > 0);
   for (int lc=0; lc<(int)storage.AT(ml).AT(rl).size(); ++lc) {
     fdata & fdatas = storage.AT(ml).AT(rl).AT(lc);
-    for (int tl=0; tl<ntl/2; ++tl) {
-      const int tl1 =         tl;
-      const int tl2 = ntl-1 - tl;
+    for (int tl=1; tl<(ntl+1)/2; ++tl) {
+      const int tl1 =       tl;
+      const int tl2 = ntl - tl;
       assert (tl1 < tl2);
       gdata * const tmpdata = fdatas.AT(tl1);
       fdatas.AT(tl1) = fdatas.AT(tl2);
@@ -383,15 +393,15 @@ ref_bnd_prolongate_all (comm_state & state,
 void
 ggf::
 mg_restrict_all (comm_state & state,
-             int const tl, int const rl, int const ml,
+                 int const tl, int const rl, int const ml,
                  CCTK_REAL const time)
 {
   static Timer timer ("mg_restrict_all");
   timer.start ();
   // Require same times
   static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml-1))
-	  <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml))));
+  assert (abs(t.get_time(ml,rl,0) - t.get_time(ml-1,rl,0))
+	  <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0))));
   vector<int> const tl2s(1,tl);
   transfer_from_all (state,
                      tl  ,rl,ml,
@@ -414,8 +424,8 @@ mg_prolongate_all (comm_state & state,
   timer.start ();
   // Require same times
   static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml+1))
-	  <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml))));
+  assert (abs(t.get_time(ml,rl,0) - t.get_time(ml+1,rl,0))
+	  <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0))));
   vector<int> const tl2s(1,tl);
   transfer_from_all (state,
                      tl  ,rl,ml,
@@ -431,22 +441,19 @@ mg_prolongate_all (comm_state & state,
 void
 ggf::
 ref_restrict_all (comm_state & state,
-                  int const tl, int const rl, int const ml,
-                  CCTK_REAL const time)
+                  int const tl, int const rl, int const ml)
 {
   // Require same times
   static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature");
-  assert (abs(t.get_time(rl,ml) - t.get_time(rl+1,ml))
-	  <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml))));
+  assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <=
+          1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl))));
   if (transport_operator == op_none or transport_operator == op_sync) return;
   static Timer timer ("ref_restrict_all");
   timer.start ();
-  vector<int> const tl2s(1,tl);
   transfer_from_all (state,
-                     tl  ,rl  ,ml,
+                     tl,rl  ,ml,
                      & dh::fast_dboxes::fast_ref_rest_sendrecv,
-                     tl2s,rl+1,ml,
-                     time);
+                     tl,rl+1,ml);
   timer.stop (0);
 }
 
@@ -522,7 +529,7 @@ transfer_from_all (comm_state & state,
     for (size_t j=0; j<i; ++j) {
       assert (tl2s.AT(i) != tl2s.AT(j));
     }
-    times.AT(i) = t.time(tl2s.AT(i),rl2,ml2);
+    times.AT(i) = t.get_time(ml2,rl2,tl2s.AT(i));
   }
   
   // Interpolation orders
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 43c5c2e93..5538f3941 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -105,6 +105,9 @@ public:
 
   // Cycle the time levels by rotating the data sets
   void cycle_all (int rl, int ml);
+
+  // Uncycle the time levels by rotating the data sets
+  void uncycle_all (int rl, int ml);
   
   // Flip the time levels by exchanging the data sets
   void flip_all (int rl, int ml);
@@ -135,7 +138,7 @@ public:
 
   // Restrict a refinement level
   void ref_restrict_all (comm_state& state,
-                         int tl, int rl, int ml, CCTK_REAL time);
+                         int tl, int rl, int ml);
 
   // Prolongate a refinement level
   void ref_prolongate_all (comm_state& state,
@@ -171,7 +174,7 @@ protected:
   {
     vector <int> tl2s(1);
     tl2s.AT(0) = tl2;
-    CCTK_REAL const time = t.time (tl2,rl2,ml2);
+    CCTK_REAL const time = t.get_time (ml2,rl2,tl2);
     transfer_from_all (state,
                        tl1, rl1, ml1,
                        sendrecvs,
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 21cb399a6..05c1b55ea 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -549,7 +549,7 @@ gh::
 output (ostream & os)
   const
 {
-  os << "gh:"
+  os << "gh:{"
      << "reffacts=" << reffacts << ",refcentering=" << refcent << ","
      << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
      << "baseextents=" << baseextents << ","
@@ -565,6 +565,6 @@ output (ostream & os)
       os << *d;
     }
   }
-  os << "}";
+  os << "}}";
   return os;
 }
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index dd441e455..bed473245 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -19,8 +19,10 @@ list<th*> th::allth;
 
 
 // Constructors
-th::th (gh& h_, const vector<int> & reffacts_, const CCTK_REAL basedelta)
-  : h(h_), reffacts(reffacts_), delta(basedelta)
+th::th (gh& h_, int const timelevels_, vector<int> const& reffacts_,
+        bool const time_interpolation_during_regridding_)
+  : h(h_), timelevels(timelevels_), reffacts(reffacts_),
+    time_interpolation_during_regridding (time_interpolation_during_regridding_)
 {
   assert (reffacts.size() >= 1);
   assert (reffacts.front() == 1);
@@ -42,6 +44,9 @@ th::~th ()
 // Modifiers
 void th::regrid ()
 {
+  CCTK_REAL const basetime  = 0.0;
+  CCTK_REAL const basedelta = 1.0;
+  
   const int old_mglevels = times.size();
   times.resize(h.mglevels());
   deltas.resize(h.mglevels());
@@ -50,17 +55,64 @@ void th::regrid ()
     times.AT(ml).resize(h.reflevels());
     deltas.AT(ml).resize(h.reflevels());
     for (int rl=0; rl<h.reflevels(); ++rl) {
+      if (ml==0) {
+	deltas.AT(ml).AT(rl) = basedelta / reffacts.AT(rl);
+      } else {
+	deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact;
+      }
       if (old_mglevels==0) {
-        times.AT(ml).AT(rl) = 0;
+        times.AT(ml).AT(rl).resize(timelevels);
+        for (int tl=0; tl<timelevels; ++tl) {
+          times.AT(ml).AT(rl).AT(tl) = basetime - tl * deltas.AT(ml).AT(rl);
+        }
       } else if (rl < old_reflevels) {
         // do nothing
+      } else if (rl == 0) {
+        assert (0);
+        times.AT(ml).AT(rl).resize(timelevels);
+        for (int tl=0; tl<timelevels; ++tl) {
+          times.AT(ml).AT(rl).AT(tl) = basetime - tl * deltas.AT(ml).AT(rl);
+        }
       } else {
-        times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1);
+        if (time_interpolation_during_regridding) {
+#warning "We probably don't want to do this, but it is nice for compatibility"
+          times.AT(ml).AT(rl).resize(timelevels);
+          for (int tl=0; tl<timelevels; ++tl) {
+            // linear interpolation between the two surrounding coarse
+            // grid times
+            assert (reffacts.AT(rl) % reffacts.AT(rl-1) == 0);
+            int const rf = reffacts.AT(rl) / reffacts.AT(rl-1);
+            int const ctl = tl / rf;
+            int const mtl = tl % rf;
+            if (mtl == 0) {
+              times.AT(ml).AT(rl).AT(tl) = times.AT(ml).AT(rl-1).AT(ctl);
+            } else {
+              assert (ctl+1<timelevels);
+              CCTK_REAL const alpha = (CCTK_REAL)mtl / rf;
+              assert (alpha>0 and alpha<1);
+              times.AT(ml).AT(rl).AT(tl) =
+                (1-alpha) * times.AT(ml).AT(rl-1).AT(ctl  ) +
+                (  alpha) * times.AT(ml).AT(rl-1).AT(ctl+1);
+            }
+          }
+        } else {
+          times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1);
+        }
       }
-      if (ml==0) {
-	deltas.AT(ml).AT(rl) = delta / reffacts.AT(rl);
-      } else {
-	deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact;
+    }
+  }
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      for (int tl=1; tl<timelevels; ++tl) {
+        assert (times.AT(ml).AT(rl).AT(tl) < times.AT(ml).AT(rl).AT(tl-1));
+      }
+    }
+  }
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      // assert (isfinite(deltas.AT(ml).AT(rl)));
+      for (int tl=0; tl<timelevels; ++tl) {
+        // assert (isfinite(times.AT(ml).AT(rl).AT(tl)));
       }
     }
   }
@@ -72,6 +124,37 @@ void th::regrid_free ()
 
 
 
+void th::advance_time (int const ml, int const rl)
+{
+  for (int tl=timelevels-1; tl>0; --tl) {
+    set_time (ml,rl,tl, get_time(ml,rl,tl-1));
+  }
+  set_time (ml,rl,0, get_time(ml,rl,0) + get_delta(ml,rl));
+}
+
+void th::retreat_time (int const ml, int const rl)
+{
+  CCTK_REAL const t = get_time(ml,rl,0);
+  for (int tl=0; tl<timelevels-1; ++tl) {
+    set_time (ml,rl,tl, get_time(ml,rl,tl+1));
+  }
+  set_time (ml,rl,timelevels-1, t);
+}
+
+void th::flip_timelevels (int const ml, int const rl)
+{
+  for (int tl=1; tl<(timelevels+1)/2; ++tl) {
+    int const tl2 = timelevels - tl;
+    CCTK_REAL const t  = get_time(ml,rl,tl );
+    CCTK_REAL const t2 = get_time(ml,rl,tl2);
+    set_time (ml,rl,tl ,t2);
+    set_time (ml,rl,tl2,t );
+  }
+  set_delta (ml,rl, -get_delta(ml,rl));
+}
+
+
+
 // Memory usage
 size_t
 th::
@@ -80,7 +163,6 @@ memory ()
 {
   return
     memoryof (reffacts) +
-    memoryof (delta) +
     memoryof (times) +
     memoryof (deltas);
 }
@@ -100,19 +182,33 @@ allmemory ()
 
 
 
+// Input
+istream& th::input (istream& is)
+{
+  skipws (is);
+  consume (is, "th:{");
+  consume (is, "timelevels=");
+  is >> timelevels;
+  consume (is, ",");
+  consume (is, "reffacts=");
+  is >> reffacts;
+  consume (is, ",");
+  consume (is, "times=");
+  is >> times;
+  consume (is, ",");
+  consume (is, "deltas=");
+  is >> deltas;
+  consume (is, "}");
+  return is;
+}
+
 // Output
-void th::output (ostream& os) const
+ostream& th::output (ostream& os) const
 {
-  os << "th:"
-     << "times={";
-  const char * sep = "";
-  for (int ml=0; ml<h.mglevels(); ++ml) {
-    for (int rl=0; rl<h.reflevels(); ++rl) {
-      os << sep
-         << ml << ":" << rl << ":"
-	 << times.AT(ml).AT(rl) << "(" << deltas.AT(ml).AT(rl) << ")";
-      sep = ",";
-    }
-  }
-  os << "}";
+  os << "th:{"
+     << "timelevels=" << timelevels << ","
+     << "reffacts=" << reffacts << ","
+     << "times=" << times << ","
+     << "deltas=" << deltas << "}";
+  return os;
 }
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index f1a7b72b7..c1453cf4d 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -2,10 +2,11 @@
 #define TH_HH
 
 #include <cassert>
+#include <cmath>
 #include <iostream>
 #include <vector>
 
-#include "cctk.h"
+#include <cctk.h>
 
 #include "defs.hh"
 #include "gh.hh"
@@ -17,6 +18,8 @@ using namespace std;
 // Forward declaration
 class th;
 
+// Input
+istream& operator>> (istream& is, th& t);
 // Output
 ostream& operator<< (ostream& os, const th& t);
 
@@ -34,18 +37,22 @@ public:				// should be readonly
   gh& h;                        // hierarchy
   gh::th_handle gh_handle;
   
+  int timelevels;               // const
+  
 private:
   
-  const vector<int> reffacts;
-
-  CCTK_REAL delta;		// time step
-  vector<vector<CCTK_REAL> > times; // current times [ml][rl]
-  vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl]
+  vector<int> reffacts;         // const
+  
+  bool const time_interpolation_during_regridding;
+  
+  vector<vector<vector<CCTK_REAL> > > times; // current times [ml][rl][tl]
+  vector<vector<CCTK_REAL> > deltas;         // time steps [ml][rl]
   
 public:
   
   // Constructors
-  th (gh& h, const vector<int> & reffacts, const CCTK_REAL basedelta);
+  th (gh& h, int timelevels, vector<int> const& reffacts,
+      bool time_interpolation_during_regridding);
   
   // Destructors
   ~th ();
@@ -55,50 +62,52 @@ public:
   void regrid_free ();
   
   // Time management
-  CCTK_REAL get_time (const int rl, const int ml) const CCTK_ATTRIBUTE_PURE
+  void set_time (int const ml, int const rl, int const tl, CCTK_REAL const& t)
   {
-    assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    return times.AT(ml).AT(rl);
-  }
-  
-  void set_time (const int rl, const int ml, const CCTK_REAL t)
-  {
     assert (rl>=0 and rl<h.reflevels());
-    assert (ml>=0 and ml<h.mglevels());
-    times.AT(ml).AT(rl) = t;
-  }
-  
-  void advance_time (const int rl, const int ml)
-  {
-    set_time(rl,ml, get_time(rl,ml) + get_delta(rl,ml));
+    assert (tl>=0 and tl<timelevels);
+    // assert (isfinite(t));
+    times.AT(ml).AT(rl).AT(tl) = t;
   }
   
-  CCTK_REAL get_delta (const int rl, const int ml) const CCTK_ATTRIBUTE_PURE
+  CCTK_REAL get_time (int const ml, int const rl, int const tl)
+    const CCTK_ATTRIBUTE_PURE
   {
-    assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    return deltas.AT(ml).AT(rl);
+    assert (rl>=0 and rl<h.reflevels());
+    assert (tl>=0 and tl<timelevels);
+    CCTK_REAL const t = times.AT(ml).AT(rl).AT(tl);
+    // assert (isfinite(t));
+    return t;
   }
   
-  void set_delta (const int rl, const int ml, const CCTK_REAL dt)
+  void set_delta (int const ml, int const rl, CCTK_REAL const& dt)
   {
-    assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
+    assert (rl>=0 and rl<h.reflevels());
+    // assert (isfinite(dt));
     deltas.AT(ml).AT(rl) = dt;
   }
   
-  CCTK_REAL time (const int tl, const int rl, const int ml) const CCTK_ATTRIBUTE_PURE
+  CCTK_REAL get_delta (int const ml, int const rl) const CCTK_ATTRIBUTE_PURE
   {
-    assert (rl>=0 and rl<h.reflevels());
     assert (ml>=0 and ml<h.mglevels());
-    return get_time(rl, ml) - tl * get_delta(rl, ml);
+    assert (rl>=0 and rl<h.reflevels());
+    CCTK_REAL const dt = deltas.AT(ml).AT(rl);
+    // assert (isfinite(dt));
+    return dt;
   }
   
+  void advance_time (int const ml, int const rl);
+  void retreat_time (int const ml, int const rl);
+  void flip_timelevels (int const ml, int const rl);
+  
   // Output
   size_t memory () const CCTK_ATTRIBUTE_PURE;
   static size_t allmemory () CCTK_ATTRIBUTE_PURE;
-  void output (ostream& os) const;
+  istream& input (istream& is);
+  ostream& output (ostream& os) const;
 };
 
 
@@ -108,9 +117,11 @@ inline size_t memoryof (th const & t)
 {
   return t.memory ();
 }
+inline istream& operator>> (istream& is, th& t) {
+  return t.input(is);
+}
 inline ostream& operator<< (ostream& os, const th& t) {
-  t.output(os);
-  return os;
+  return t.output(os);
 }