Changeset 8274


Ignore:
Timestamp:
Dec 8, 2011, 11:38:18 AM (12 years ago)
Author:
paul
Message:

Fixed loss of mass in sqpipe and rearranged code a little

Location:
trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/parabolic_canal.py

    r8265 r8274  
    9494    pylab.ion()
    9595
    96     x = domain.get_centroids()
    97     z = domain.get_quantity('elevation', 'centroids')
    98     w = domain.get_quantity('stage', 'centroids')
    99     h = domain.get_quantity('height', 'centroids')
    100     v = domain.get_quantity('velocity', 'centroids')
    101     t = domain.get_quantity('top', 'centroids')
    102     s = domain.state
    103     m = domain.get_mass()
    104     M = m.sum() * numpy.ones_like(x)
     96    x, z, w, h, v, t, s, m, M = get_quantities(domain)
    10597
    106     plot1 = pylab.subplot(511)
    107     zplot, = pylab.plot(x, z)
    108     wplot, = pylab.plot(x, w)
    109     ztplot, = pylab.plot(x, z+t)   
    110     plot1.set_ylim([-1,50])
    111     pylab.xlabel('Position')
    112     pylab.ylabel('Stage')
    113 
    114     plot2 = pylab.subplot(512)
    115     hplot, = pylab.plot(x, h)
    116     tplot, = pylab.plot(x, t)
    117     plot2.set_ylim([-1,20])
    118     pylab.xlabel('Position')
    119     pylab.ylabel('Height')
    120 
    121     plot3 = pylab.subplot(513)
    122     vplot, = pylab.plot(x, v)
    123     plot3.set_ylim([-30,30])
    124     pylab.xlabel('Position')
    125     pylab.ylabel('Velocity')
    126 
    127     plot4 = pylab.subplot(514)
    128     splot, = pylab.plot(x, s)
    129     plot4.set_ylim([-1,2])
    130     pylab.xlabel('Position')
    131     pylab.ylabel('State')
    132 
    133     plot5 = pylab.subplot(515)
    134     Mplot, = pylab.plot(x, m)
    135     mplot, = pylab.plot(x, M)
    136     plot5.set_ylim([-1,45000])
    137     pylab.xlabel('Position')
    138     pylab.ylabel('Mass')
    139 
     98    zplot, wplot, ztplot, hplot, tplot, vplot, splot, Mplot, mplot = make_plots(x, z, w, h, v, t, s, m, M)
    14099
    141100    for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
    142         z = domain.get_quantity('elevation', 'centroids')
    143         w = domain.get_quantity('stage', 'centroids')
    144         h = domain.get_quantity('height', 'centroids')
    145         t = domain.get_quantity('top', 'centroids')
    146         v = domain.get_quantity('velocity', 'centroids')
    147         s = domain.state
    148         m = domain.get_mass()
    149         M = m.sum() * numpy.ones_like(x)
     101        x, z, w, h, v, t, s, m, M = get_quantities(domain)
    150102       
    151103        zplot.set_ydata(z)
     
    166118    import pylab
    167119
    168     x = domain.get_centroids()
    169     z = domain.get_quantity('elevation', 'centroids')
    170     w = domain.get_quantity('stage', 'centroids')
    171     h = domain.get_quantity('height', 'centroids')
    172     t = domain.get_quantity('top', 'centroids')
    173     u = domain.get_quantity('velocity', 'centroids')
     120    x, z, w, h, v, t, s, m, M = get_quantities(domain)
    174121
    175122    pylab.ioff()
    176123    pylab.hold(False)
    177124
    178     plot1 = pylab.subplot(311)
    179     pylab.plot(x, z, x, w, x, z+t)
    180     plot1.set_ylim([-1,40])
    181     pylab.xlabel('Position')
    182     pylab.ylabel('Stage')
    183 
    184     plot2 = pylab.subplot(312)
    185     pylab.plot(x, h, x, t)
    186     plot2.set_ylim([-1,10])
    187     pylab.xlabel('Position')
    188     pylab.ylabel('Height')
    189 
    190     plot3 = pylab.subplot(313)
    191     pylab.plot(x, u)
    192     plot3.set_ylim([-15,15])
    193     pylab.xlabel('Position')
    194     pylab.ylabel('Velocity')
    195 
     125    make_plots(x, z, w, h, v, t, s, m, M)
     126   
    196127    pylab.show()
    197128
     
    205136        f.write("%s %s %s\n" % (x[i], z[i], w[i]))
    206137    f.close()
     138
     139def get_quantities(domain):
     140    x = domain.get_centroids()
     141    z = domain.get_quantity('elevation', 'centroids')
     142    w = domain.get_quantity('stage', 'centroids')
     143    h = domain.get_quantity('height', 'centroids')
     144    v = domain.get_quantity('velocity', 'centroids')
     145    t = domain.get_quantity('top', 'centroids')
     146    s = domain.state
     147    m = domain.get_mass()
     148    M = m.sum() * numpy.ones_like(x)
     149
     150    return x, z, w, h, v, t, s, m, M
     151
     152def make_plots(x, z, w, h, v, t, s, m, M):
     153    import pylab
     154   
     155    plot1 = pylab.subplot(321)
     156    zplot, = pylab.plot(x, z)
     157    wplot, = pylab.plot(x, w)
     158    ztplot, = pylab.plot(x, z+t)   
     159    plot1.set_ylim([-1,50])
     160    pylab.xlabel('Position')
     161    pylab.ylabel('Stage')
     162
     163    plot2 = pylab.subplot(322)
     164    hplot, = pylab.plot(x, h)
     165    tplot, = pylab.plot(x, t)
     166    plot2.set_ylim([-1,20])
     167    pylab.xlabel('Position')
     168    pylab.ylabel('Height')
     169
     170    plot3 = pylab.subplot(323)
     171    vplot, = pylab.plot(x, v)
     172    plot3.set_ylim([-30,30])
     173    pylab.xlabel('Position')
     174    pylab.ylabel('Velocity')
     175
     176    plot4 = pylab.subplot(324)
     177    splot, = pylab.plot(x, s)
     178    plot4.set_ylim([-1,2])
     179    pylab.xlabel('Position')
     180    pylab.ylabel('State')
     181
     182    plot5 = pylab.subplot(325)
     183    mplot, = pylab.plot(x, m)
     184    plot5.set_ylim([-1,1000])
     185    pylab.xlabel('Position')
     186    pylab.ylabel('Mass')
     187
     188    plot6 = pylab.subplot(326) Mplot, = pylab.plot(x, M)
     189    plot6.set_ylim([-1,45000]) pylab.xlabel('Position')
     190    pylab.ylabel('Total Mass')
     191
     192    return zplot, wplot, ztplot, hplot, tplot, vplot, splot, Mplot, mplot
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe.h

    r8238 r8274  
    6969  double b;
    7070  double t;
     71  long state;
    7172};
    7273
     
    7980
    8081// Methods
     82double* sqpipe_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
     83double sqpipe_quantity_sound_speed (struct quantity *q, struct params *p);
     84double sqpipe_quantity_get_conserved (struct quantity *q, int k, double normal);
     85
     86// Helper methods
    8187double* sqpipe_quantity_free_surface_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
    8288double* sqpipe_quantity_pressurised_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
     
    8591
    8692// Implementation
     93double* sqpipe_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
     94  if (((struct sqpipe_quantity *)q)->state == 0) {
     95    return sqpipe_quantity_free_surface_flux_formula(q, normal, p, quantityflux);
     96  } else {
     97    return sqpipe_quantity_pressurised_flux_formula(q, normal, p, quantityflux);
     98  }
     99}
     100
     101double sqpipe_quantity_sound_speed (struct quantity *q, struct params *p) {
     102  if (((struct sqpipe_quantity *)q)->state == 0) {
     103    return sqpipe_quantity_free_surface_sound_speed(q, p);
     104  } else {
     105    return sqpipe_quantity_pressurised_sound_speed(q, p);
     106  }
     107}
     108
     109double sqpipe_quantity_get_conserved (struct quantity *q, int k, double normal) {
     110  struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
     111  double c;
     112 
     113  switch (k) {
     114  case 0:
     115    c = p->a;
     116    break;
     117  case 1:
     118    // This should be normal^2 p->d and normal is +/- 1
     119    c = p->d;
     120    break;
     121  default:
     122    c = 0;
     123  }
     124  return c;
     125}
     126
     127
    87128// Flux formula taking into account the normal
    88129// Note u, d should be normal*u, normal*d but since normal is +/- 1
     
    121162void sqpipe_quantity_init(struct sqpipe_quantity *q) {
    122163  static struct quantity_vtable vtable = {
    123     NULL,
    124     NULL,
    125     NULL,
     164    &sqpipe_quantity_flux_formula,
     165    &sqpipe_quantity_sound_speed,
     166    &sqpipe_quantity_get_conserved,
    126167    &sqpipe_quantity_destroy
    127168  };
     
    136177  q->b = 0.0;
    137178  q->t = 0.0;
     179  q->state = 0;
    138180}
    139181
     
    237279struct sqpipe_cell {
    238280  struct cell super;
     281 
     282  long state;
    239283};
    240284
     
    248292// Methods
    249293double sqpipe_cell_extreme_sound_speeds (struct cell *c, struct params *p);
     294double* sqpipe_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
     295
     296// Helper methods
    250297double* sqpipe_cell_free_surface_forcing_terms (struct cell *c, struct params *p, double *cellforce);
    251298double* sqpipe_cell_pressurised_forcing_terms (struct cell *c, struct params *p, double *cellforce);
     
    270317}
    271318
     319double* sqpipe_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
     320  if (((struct sqpipe_cell *)c)->state == 0) {
     321    return sqpipe_cell_free_surface_forcing_terms(c, p, cellforce);   
     322  } else {
     323    return sqpipe_cell_pressurised_forcing_terms(c, p, cellforce);
     324  }
     325}
     326
    272327double* sqpipe_cell_free_surface_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
    273328  struct sqpipe_quantity *qin = (struct sqpipe_quantity *) c->edges[0]->qout;
     
    297352    &cell_flux_function_generic,
    298353    &sqpipe_cell_extreme_sound_speeds,
    299     NULL,
     354        &sqpipe_cell_forcing_terms,
    300355    &sqpipe_cell_destroy
    301356  };
     
    309364
    310365  c->super.num_edges = 2;
     366  ((struct sqpipe_cell *)c)->state = 0;
    311367}
    312368
     
    365421// Methods
    366422struct quantity* sqpipe_quantities_get_quantity(struct quantities *qs, int i, struct quantity *q);
     423void sqpipe_quantities_update (struct quantities *qs, double *flux, int k);
    367424
    368425// Implementation
     
    383440}
    384441
     442void sqpipe_quantities_update (struct quantities *qs, double *flux, int k) {
     443  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
     444  ps->a[k] = flux[0];
     445  ps->d[k] = flux[1];
     446}
     447
     448
    385449void sqpipe_quantities_init(struct sqpipe_quantities *q) {
    386450  static struct quantities_vtable vtable = {
    387451    &sqpipe_quantities_get_quantity,
    388     NULL,
     452    &sqpipe_quantities_update,
    389453    &sqpipe_quantities_destroy
    390454  };
     
    426490struct sqpipe_domain {
    427491  struct domain super;
     492  long *state;
    428493};
    429494
    430495// Construction
    431 void sqpipe_domain_init(struct sqpipe_domain *q);
     496void sqpipe_domain_init(struct sqpipe_domain *d);
    432497struct sqpipe_domain* sqpipe_domain_new();
    433498
    434499// Destruction
    435 struct domain* sqpipe_domain_destroy(struct domain *D);
     500struct domain* sqpipe_domain_destroy(struct domain *d);
    436501
    437502// Methods
    438 void sqpipe_domain_init(struct sqpipe_domain *D);
     503struct cell* sqpipe_domain_get_cell(struct domain *D, int k, struct cell *c);
    439504
    440505// Implementation
     506struct cell* sqpipe_domain_get_cell(struct domain *D, int k, struct cell *c) {
     507  struct sqpipe_cell *sc = (struct sqpipe_cell *) c;
     508  struct sqpipe_domain *sD = (struct sqpipe_domain *) D;
     509  int i, ki, n;
     510  long state_in, state_out, state;
     511 
     512  // First get a generic cell
     513  c = domain_get_cell_generic(D, k, c);
     514
     515  // Set state
     516  // The state should be set to state_in for qin and state_out for qout
     517  // Temporarily, I set the state of the edge, aribtrarily choosing
     518  // pressurised for transition points
     519  // This will be fixed when the flux is computed at the edge rather than
     520  // at the quantity
     521  sc->state = sD->state[k];
     522  for (i = 0; i<c->num_edges; i++) {
     523        // The inner state is that of this cell
     524        state_in = sc->state;
     525        // The outer state is that of the corresponding neighbour cell
     526        ki = k*(c->num_edges) + i;
     527        n = D->neighbours[ki];
     528        state_out = sD->state[n];
     529       
     530        // If both in and out state are the same, choose that state
     531        // Otherwise, choose pressurised
     532        if (state_in == state_out) {
     533                state = state_in;
     534        } else {
     535                state = 1;
     536        }
     537    ((struct sqpipe_quantity *)(c->edges[i]->qin))->state = state;     
     538    ((struct sqpipe_quantity *)(c->edges[i]->qout))->state = state;
     539  }
     540 
     541  return c;
     542}
     543
     544
    441545void sqpipe_domain_init(struct sqpipe_domain *d) {
    442546  static struct domain_vtable vtable = {
    443547    NULL,
    444548    NULL,
    445     &domain_compute_fluxes_generic,
    446     &domain_get_cell_generic,
     549    &domain_compute_fluxes_cells_generic,
     550    &sqpipe_domain_get_cell,
    447551    &sqpipe_domain_destroy
    448552  };
     
    463567  ((struct domain *)d)->max_speed_array = NULL;
    464568
     569  d->state = NULL;
    465570  // Can't set params to NULL here since sqpipe_domain_new() has allocated
    466571  // memory
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_comp_flux.c

    r8254 r8274  
    2525  struct domain *D = (struct domain *) sqpipe_default_domain_new();
    2626  double timestep;
    27   PyArrayObject *state;
    2827  PyObject *domain;
    2928
     
    3433  }
    3534
     35  // Generic sqpipe domain
    3636  D = (struct domain *) sqpipe_domain_python_get((struct sqpipe_domain *)D, domain, timestep);
    3737
    3838  // We need the state
    39   state = get_consecutive_array(domain, "state");
    40   ((struct sqpipe_default_domain *)D)->state = (long *) state->data;
     39  //state = get_consecutive_array(domain, "state");
     40  //((struct sqpipe_default_domain *)D)->state = (long *) state->data;
    4141
     42  // Compute the fluxes and return the new timestep
    4243  timestep = domain_compute_fluxes(D);
    4344  domain_destroy(D);
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_default.h

    r8254 r8274  
    77struct sqpipe_default_quantity {
    88  struct sqpipe_quantity super;
    9   long state;
    109};
    1110
     
    1716
    1817// Methods
    19 double* sqpipe_default_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
    20 double sqpipe_default_quantity_sound_speed (struct quantity *q, struct params *p);
    21 double sqpipe_default_quantity_get_conserved (struct quantity *q, int k, double normal);
    2218
    2319// Implementation
    24 double* sqpipe_default_quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
    25   if (((struct sqpipe_default_quantity *)q)->state == 0) {
    26     return sqpipe_quantity_free_surface_flux_formula(q, normal, p, quantityflux);
    27   } else {
    28     return sqpipe_quantity_pressurised_flux_formula(q, normal, p, quantityflux);
    29   }
    30 }
    31 
    32 double sqpipe_default_quantity_sound_speed (struct quantity *q, struct params *p) {
    33   if (((struct sqpipe_default_quantity *)q)->state == 0) {
    34     return sqpipe_quantity_free_surface_sound_speed(q, p);
    35   } else {
    36     return sqpipe_quantity_pressurised_sound_speed(q, p);
    37   }
    38 }
    39 
    40 double sqpipe_default_quantity_get_conserved (struct quantity *q, int k, double normal) {
    41   struct sqpipe_quantity *p = (struct sqpipe_quantity*) q;
    42   double c;
    43  
    44   switch (k) {
    45   case 0:
    46     c = p->a;
    47     break;
    48   case 1:
    49     // This should be normal^2 p->d and normal is +/- 1
    50     c = p->d;
    51     break;
    52   default:
    53     c = 0;
    54   }
    55   return c;
    56 }
    57 
    5820void sqpipe_default_quantity_init (struct sqpipe_default_quantity *q) {
    5921  struct sqpipe_quantity *p = (struct sqpipe_quantity *)q;
    6022
    61   sqpipe_quantity_init(p);
    62 
    63   static struct quantity_vtable vtable = {
    64     &sqpipe_default_quantity_flux_formula,
    65     &sqpipe_default_quantity_sound_speed,
    66     &sqpipe_default_quantity_get_conserved,
    67     &sqpipe_quantity_destroy
    68   };
    69   p->super.vtable = &vtable;
    70 
    71   ((struct sqpipe_default_quantity *)p)->state = 0;
     23  sqpipe_quantity_init(p); 
    7224}
    7325
     
    9850
    9951// Implementation
    100 
    10152void sqpipe_default_edge_init (struct sqpipe_default_edge *q) {
    10253  sqpipe_edge_init((struct sqpipe_edge *)q);
     
    12273struct sqpipe_default_cell {
    12374  struct sqpipe_cell super;
    124 
    125   long state;
    12675};
    12776
     
    13382
    13483// Methods
    135 double* sqpipe_default_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
    13684
    13785// Implementation
    138 double* sqpipe_default_cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
    139   if (((struct sqpipe_default_cell *)c)->state == 0) {
    140     return sqpipe_cell_free_surface_forcing_terms(c, p, cellforce);   
    141   } else {
    142     return sqpipe_cell_pressurised_forcing_terms(c, p, cellforce);
    143   }
    144 }
    145 
    146 void sqpipe_default_cell_init (struct sqpipe_default_cell *q) {
    147   struct sqpipe_cell *p = (struct sqpipe_cell *) q;
     86void sqpipe_default_cell_init (struct sqpipe_default_cell *c) {
     87  struct sqpipe_cell *p = (struct sqpipe_cell *) c;
    14888
    14989  sqpipe_cell_init(p);
    150 
    151   static struct cell_vtable vtable = {
    152     &cell_flux_function_generic,
    153     &sqpipe_cell_extreme_sound_speeds,
    154     &sqpipe_default_cell_forcing_terms,
    155     &sqpipe_cell_destroy
    156   };
    157   p->super.vtable = &vtable;
    158 
    159   ((struct sqpipe_default_cell *)p)->state = 0;
    16090}
    16191
     
    193123
    194124// Methods
    195 void sqpipe_default_quantities_update (struct quantities *qs, double *flux, int k);
    196 
    197125
    198126// Implementation
    199 void sqpipe_default_quantities_update (struct quantities *qs, double *flux, int k) {
    200   struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
    201   ps->a[k] = flux[0];
    202   ps->d[k] = flux[1];
    203 }
    204 
    205127void sqpipe_default_quantities_init (struct sqpipe_default_quantities *q) {
    206128  struct sqpipe_quantities *p = (struct sqpipe_quantities *)q;
    207129
    208130  sqpipe_quantities_init(p);
    209 
    210   static struct quantities_vtable vtable = {
    211     &sqpipe_quantities_get_quantity,
    212     &sqpipe_default_quantities_update,
    213     &sqpipe_quantities_destroy
    214   };
    215   p->super.vtable = &vtable;
    216131}
    217132
     
    232147struct sqpipe_default_domain {
    233148  struct sqpipe_domain super;
    234   long *state;
    235149};
    236150
     
    244158struct edge* sqpipe_default_domain_new_edge(struct domain *d);
    245159struct cell* sqpipe_default_domain_new_cell(struct domain *d);
    246 struct cell* sqpipe_default_domain_get_cell(struct domain *D, int k, struct cell *c);
    247160
    248161// Implementation
     
    258171}
    259172
    260 struct cell* sqpipe_default_domain_get_cell(struct domain *D, int k, struct cell *c) {
    261   struct sqpipe_default_cell *sc = (struct sqpipe_default_cell *) c;
    262   struct sqpipe_default_domain *sD = (struct sqpipe_default_domain *) D;
    263   int i;
    264  
    265   c = domain_get_cell_generic(D, k, c);
    266 
    267   // Set state
    268   sc->state = sD->state[k];
    269   for (i = 0; i<c->num_edges; i++) {
    270     ((struct sqpipe_default_quantity *)(c->edges[i]->qin))->state = sc->state;
    271     ((struct sqpipe_default_quantity *)(c->edges[i]->qout))->state = sc->state;
    272   }
    273  
    274   return c;
    275 }
    276 
    277173void sqpipe_default_domain_init (struct sqpipe_default_domain *d) {
    278174  struct sqpipe_domain *p = (struct sqpipe_domain *)d;
     
    283179    &sqpipe_default_domain_new_edge,
    284180    &sqpipe_default_domain_new_cell,
    285     &domain_compute_fluxes_generic,
    286     &sqpipe_default_domain_get_cell,
     181    &domain_compute_fluxes_cells_generic,
     182    &sqpipe_domain_get_cell,
    287183    &sqpipe_domain_destroy
    288184  };
    289185  p->super.vtable = &vtable;
    290 
    291   ((struct sqpipe_default_domain *)p)->state = NULL;
    292186}
    293187
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py

    r8265 r8274  
    270270        """
    271271
    272         # Equivalent area A is defined by A = (\rho/\rho_0) A_{pipe}
    273         # A_{pipe} = width * top
    274         # A = width * height
    275         # \rho_0 = 1
    276         # Therefore \rho = height/top for pressurised states
    277         # For free surface \rho = 1
    278 
    279         # Coarse approximation of mass M in a cell, assume rho is constant
    280         # M = \Delta x * \rho * A
    281         #   = \Delta_x * A * height/top (pressurised)
    282         #   = \Delta_x * A (un-pressurised)
    283 
    284         area      = self.quantities['area']       
    285         height    = self.quantities['height']       
    286         top       = self.quantities['top']
    287        
    288         #Arrays   
    289         a   = area.centroid_values       
    290         h   = height.centroid_values
    291         t   = top.centroid_values
    292        
    293         M = numpy.where(self.state == 0, self.areas * a, self.areas * a * h/t)
    294 
    295         return M
     272        # The equivalent area is the conserved mass quantity
     273        # It is equal to \rho/\rho_0 A where A is the cross sectional
     274        # area of the fluid where $\rho_0 is some fixed reference density
     275        area = self.quantities['area']
     276        a = area.centroid_values
     277       
     278        return a * self.areas
     279       
    296280
    297281# Auxillary methods
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_python.h

    r8238 r8274  
    44struct sqpipe_domain* sqpipe_domain_python_get(struct sqpipe_domain *sD, PyObject *domain, double timestep) {
    55  PyObject *quantities;
     6  PyArrayObject *state;
    67  struct domain *D = (struct domain *)sD;
    78  struct sqpipe_params *sp = (struct sqpipe_params *) D->params;
    8 
     9 
    910  // Get a generic domain
    1011  D = get_python_domain(D, domain, timestep);
     
    1516  D->boundary_values = sqpipe_quantities_python_get(D->boundary_values, quantities, "boundary_values");
    1617  D->explicit_update = sqpipe_quantities_python_get(D->explicit_update, quantities, "explicit_update");
     18
     19  // Get cell level quantities
     20  state = get_consecutive_array(domain, "state");
     21  sD->state = (long *) state->data;
    1722
    1823  // Get the specific parameters for square pipes
  • trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/test_sqpipe_domain.py

    r8263 r8274  
    55
    66domain = anuga_1d.sqpipe.parabolic_canal.get_domain(dom)
    7 finaltime = 10000.0
     7finaltime = 1000.0
    88yieldstep = 10.0
     9#m0 = (domain.get_mass()).sum()
     10#for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime):
     11#       m = (domain.get_mass()).sum()
     12#       if abs(m - m0) > 0.5:
     13#               print m0 - m
     14#               
     15#       m0 = m
     16       
     17#outfile="out"
     18#anuga_1d.sqpipe.parabolic_canal.write_domain(domain, outfile)
     19
    920anuga_1d.sqpipe.parabolic_canal.animate_domain(domain, yieldstep, finaltime)
     21print "finished"
    1022anuga_1d.sqpipe.parabolic_canal.plot_domain(domain)
Note: See TracChangeset for help on using the changeset viewer.