source: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe.h @ 8274

Last change on this file since 8274 was 8274, checked in by paul, 12 years ago

Fixed loss of mass in sqpipe and rearranged code a little

File size: 16.9 KB
Line 
1/**********************************************************
2Pipe params class
3**********************************************************/
4
5struct sqpipe_params {
6  struct params super;
7
8  double g;
9  double h0;
10  double bulk_modulus;
11};
12
13// Construction
14void sqpipe_params_init(struct sqpipe_params *q);
15struct sqpipe_params* sqpipe_params_new();
16
17// Destruction
18struct params* sqpipe_params_destroy(struct params *q);
19
20// Methods
21
22// Implementation
23void sqpipe_params_init(struct sqpipe_params *q) {
24  static struct params_vtable vtable = {
25    &sqpipe_params_destroy
26  };
27  q->super.vtable = &vtable;
28 
29  ((struct params *)q)->cfl = 0.0; 
30  ((struct params *)q)->epsilon = 0.0;
31  q->g = 0.0;
32  q->h0 = 0.0;
33}
34
35struct sqpipe_params* sqpipe_params_new() {
36  struct sqpipe_params *p = malloc(sizeof(struct sqpipe_params));
37
38  sqpipe_params_init(p);
39
40  return p; 
41}
42
43struct params* sqpipe_params_destroy(struct params *q) {
44  free(q);
45  q = NULL;
46 
47  return q;
48}
49/**********************************************************/
50
51
52/**********************************************************
53Pipe quantity class
54This is a virtual class and should never be created.
55The inherting class needs to implement
56    init, new, flux_formula0, flux_formula1, sound_speed,
57    get_conserved
58**********************************************************/
59
60struct sqpipe_quantity {
61  struct quantity super;
62
63  double a;
64  double d;
65  double w;
66  double h;
67  double u;
68  double z; 
69  double b;
70  double t;
71  long state;
72};
73
74// Construction
75void sqpipe_quantity_init(struct sqpipe_quantity *q);
76struct sqpipe_quantity* sqpipe_quantity_new();
77
78// Destruction
79struct quantity* sqpipe_quantity_destroy(struct quantity *q);
80
81// 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
87double* sqpipe_quantity_free_surface_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
88double* sqpipe_quantity_pressurised_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
89double sqpipe_quantity_free_surface_sound_speed (struct quantity *q, struct params *p);
90double sqpipe_quantity_pressurised_sound_speed (struct quantity *q, struct params *p);
91
92// 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
128// Flux formula taking into account the normal
129// Note u, d should be normal*u, normal*d but since normal is +/- 1
130// (normal*u) * (normal*d) = u*d
131double* sqpipe_quantity_free_surface_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
132  struct sqpipe_quantity *pq = (struct sqpipe_quantity*) q;
133  struct sqpipe_params *sp = (struct sqpipe_params *) p;
134
135  quantityflux[0] = (normal * pq->u) * quantity_get_conserved(q, 0, normal);
136  quantityflux[1] = normal * (pq->u * quantity_get_conserved(q, 1, normal) + 0.5 * sp->g * pq->h *pq->h);
137
138  return quantityflux;
139}
140
141double* sqpipe_quantity_pressurised_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
142  struct sqpipe_quantity *pq = (struct sqpipe_quantity*) q;
143  struct sqpipe_params *sp = (struct sqpipe_params *) p;
144 
145  quantityflux[0] = (normal * pq->u) * pq->a; 
146  quantityflux[1] = normal * (pq->u * pq->d + sp->bulk_modulus*(pq->h - pq->t) + 0.5 * sp->g * pq->t * pq->t);
147
148  return quantityflux;
149}
150
151double sqpipe_quantity_free_surface_sound_speed (struct quantity *q, struct params *p) {
152  struct sqpipe_quantity *sq = (struct sqpipe_quantity*) q;
153  struct sqpipe_params *sp = (struct sqpipe_params*) p; 
154
155  return sqrt(sp->g * sq->h);
156}
157
158double sqpipe_quantity_pressurised_sound_speed (struct quantity *q, struct params *p) {
159  return ((struct sqpipe_params *)p)->bulk_modulus;
160}
161
162void sqpipe_quantity_init(struct sqpipe_quantity *q) {
163  static struct quantity_vtable vtable = {
164    &sqpipe_quantity_flux_formula,
165    &sqpipe_quantity_sound_speed,
166    &sqpipe_quantity_get_conserved,
167    &sqpipe_quantity_destroy
168  };
169  q->super.vtable = &vtable;
170 
171  q->a = 0.0;
172  q->d = 0.0;
173  q->w = 0.0;
174  q->h = 0.0;
175  q->u = 0.0;
176  q->z = 0.0; 
177  q->b = 0.0;
178  q->t = 0.0;
179  q->state = 0;
180}
181
182struct sqpipe_quantity* sqpipe_quantity_new() {
183  struct sqpipe_quantity *p = malloc(sizeof(struct sqpipe_quantity));
184
185  return p; 
186}
187
188struct quantity* sqpipe_quantity_destroy(struct quantity *q) {
189  free(q);
190  q = NULL;
191 
192  return q;
193}
194/**********************************************************/
195
196
197/**********************************************************
198Pipe edge class
199This is a virtual class and should never be created.
200The inherting class needs to implement
201    init, new, flux_function
202**********************************************************/
203struct sqpipe_edge {
204  struct edge super;
205};
206
207// Construction
208void sqpipe_edge_init(struct sqpipe_edge *e);
209struct sqpipe_edge* sqpipe_edge_new();
210
211// Destruction
212struct edge* sqpipe_edge_destroy(struct edge *e);
213
214// Methods
215double sqpipe_edge_extreme_sound_speeds (struct edge *e, struct params *p);
216
217// Implementation
218double sqpipe_edge_extreme_sound_speeds (struct edge *e, struct params *p) {
219  double soundspeed_left, soundspeed_right, max_speed;
220  struct sqpipe_quantity *qin = (struct sqpipe_quantity*) e->qin;
221  struct sqpipe_quantity *qout = (struct sqpipe_quantity*) e->qout;
222
223  soundspeed_left = quantity_sound_speed(e->qin, p);
224  soundspeed_right = quantity_sound_speed(e->qout, p);
225
226  e->smax = max(e->normal * qin->u + soundspeed_left, e->normal * qout->u + soundspeed_right);
227  if (e->smax < 0.0) e->smax = 0.0;
228       
229  e->smin = min(e->normal * qin->u - soundspeed_left, e->normal * qout->u - soundspeed_right);
230  if (e->smin > 0.0) e->smin = 0.0;
231
232  if (e->smax - e->smin < p->epsilon) {
233    max_speed = 0.0;
234  } else {
235    max_speed = max(fabs(e->smin), fabs(e->smax));
236  }
237
238  return max_speed;
239}
240
241void sqpipe_edge_init(struct sqpipe_edge *e) {
242  static struct edge_vtable vtable = {
243    &edge_flux_function_godunov,
244    &sqpipe_edge_extreme_sound_speeds,
245    &sqpipe_edge_destroy
246  };
247  e->super.vtable = &vtable;
248
249  ((struct edge *)e)->qin = NULL;
250  ((struct edge *)e)->qout = NULL;
251  ((struct edge *)e)->normal = 0.0;
252  ((struct edge *)e)->smax = 0.0;
253  ((struct edge *)e)->smin = 0.0;
254}
255
256struct sqpipe_edge* sqpipe_edge_new() {
257  struct sqpipe_edge *p = malloc(sizeof(struct sqpipe_edge));
258
259  return p; 
260}
261
262struct edge* sqpipe_edge_destroy(struct edge *e) { 
263  e->qin = quantity_destroy(e->qin);
264  e->qout = quantity_destroy(e->qout);
265 
266  free(e);
267  e = NULL;
268 
269  return e;
270}
271/**********************************************************/
272
273/**********************************************************
274Pipe cell class
275This is a virtual class and should never be created.
276The inherting class needs to implement
277    init, new, flux_function, forcing_terms
278**********************************************************/
279struct sqpipe_cell {
280  struct cell super;
281 
282  long state;
283};
284
285// Construction
286void sqpipe_cell_init(struct sqpipe_cell *e);
287struct sqpipe_cell* sqpipe_cell_new();
288
289// Destruction
290struct cell* sqpipe_cell_destroy(struct cell *c);
291
292// Methods
293double 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
297double* sqpipe_cell_free_surface_forcing_terms (struct cell *c, struct params *p, double *cellforce);
298double* sqpipe_cell_pressurised_forcing_terms (struct cell *c, struct params *p, double *cellforce);
299
300// Implementation
301double sqpipe_cell_extreme_sound_speeds (struct cell *c, struct params *p) {
302  int i;
303  // The compiler doesn't know c->num_edges is positive hence whether the
304  // loop will run. This throws the warning:
305  // ‘max_speed’ may be used uninitialized in this function
306  // and the same for tmp_max_speed
307  // We initiallise them to avoid the warning since we know nothing could
308  // possibly go wrong...
309  double max_speed = 0.0, tmp_max_speed = 0.0;
310
311  for (i=0; i<c->num_edges; i++) {
312    tmp_max_speed = edge_extreme_sound_speeds(c->edges[i], p);
313    max_speed = max(max_speed, tmp_max_speed);
314  }
315
316  return max_speed;
317}
318
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
327double* sqpipe_cell_free_surface_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
328  struct sqpipe_quantity *qin = (struct sqpipe_quantity *) c->edges[0]->qout;
329  struct sqpipe_quantity *qout = (struct sqpipe_quantity *) c->edges[1]->qin;
330  struct sqpipe_params *sp = (struct sqpipe_params *) p;
331
332  cellforce[0] = 0.0; 
333  cellforce[1] = 0.5 * sp->g * (qout->h + qin->h) * (qout->z - qin->z);
334 
335  return cellforce;
336}
337
338double* sqpipe_cell_pressurised_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
339  struct sqpipe_quantity *qin = (struct sqpipe_quantity *) c->edges[0]->qout;
340  struct sqpipe_quantity *qout = (struct sqpipe_quantity *) c->edges[1]->qin;
341  struct sqpipe_params *sp = (struct sqpipe_params *) p;
342
343  cellforce[0] = 0.0;
344  cellforce[1] = sp->bulk_modulus * (sp->g * (0.5*(qout->h + qin->h) - 0.5*(qout->t + qin->t)) * (qout->b - qin->b) + sp->g * (0.5*(qout->h + qin->h) - 0.5*(qout->t + qin->t)) * 0.5*(qout->b + qin->b) * (qout->t - qin->t) * 1/(0.5*(qout->t + qin->t)));
345
346  return cellforce;
347}
348
349
350void sqpipe_cell_init(struct sqpipe_cell *c) {
351  static struct cell_vtable vtable = {
352    &cell_flux_function_generic, 
353    &sqpipe_cell_extreme_sound_speeds,
354        &sqpipe_cell_forcing_terms,
355    &sqpipe_cell_destroy
356  };
357  c->super.vtable = &vtable;
358
359  // Can't set edges to NULL here since sqpipe_cell_new() has allocated memory
360  // If a cell is created without using cell_new(), then edges will be an
361  // unitialised pointer.
362  // Can't allocate memory here since init should not allocate memory
363  // The moral: Use cell_new() to create cells!
364
365  c->super.num_edges = 2;
366  ((struct sqpipe_cell *)c)->state = 0;
367}
368
369struct sqpipe_cell* sqpipe_cell_new() {
370  struct sqpipe_cell *c = malloc(sizeof(struct sqpipe_cell));
371
372  // Allocate memory for a length 2 array of pointers to edges
373  ((struct cell *)c)->edges = (struct edge  **) malloc(2 * sizeof(struct sqpipe_edge *));
374 
375  return c;
376}
377
378struct cell* sqpipe_cell_destroy(struct cell *c) { 
379  int i;
380
381  for (i=0; i<c->num_edges; i++) {
382    c->edges[i] = edge_destroy(c->edges[i]);
383  } 
384
385  free(c->edges);
386  c->edges = NULL;
387  free(c);
388  c = NULL;
389 
390  return c;
391}
392/**********************************************************/
393
394
395/**********************************************************
396Pipe quantities class
397This is a virtual class and should never be created.
398The inherting class needs to implement:
399    init, new, get_quantities
400**********************************************************/
401struct sqpipe_quantities {
402  struct quantities super;
403
404  double *a;
405  double *d;
406  double *w;
407  double *h;
408  double *u;
409  double *z; 
410  double *b;
411  double *t;
412};
413
414// Construction
415void sqpipe_quantities_init(struct sqpipe_quantities *q);
416struct sqpipe_quantities* sqpipe_quantities_new();
417
418// Destruction
419struct quantities* sqpipe_quantities_destroy(struct quantities *qs);
420
421// Methods
422struct quantity* sqpipe_quantities_get_quantity(struct quantities *qs, int i, struct quantity *q);
423void sqpipe_quantities_update (struct quantities *qs, double *flux, int k);
424
425// Implementation
426struct quantity* sqpipe_quantities_get_quantity(struct quantities *qs, int i, struct quantity *q) { 
427  struct sqpipe_quantities *ps = (struct sqpipe_quantities *) qs;
428  struct sqpipe_quantity *p = (struct sqpipe_quantity *) q;
429
430  p->a = ps->a[i];
431  p->d = ps->d[i];
432  p->w = ps->w[i];
433  p->h = ps->h[i];
434  p->u = ps->u[i]; 
435  p->z = ps->z[i];
436  p->b = ps->b[i];
437  p->t = ps->t[i];
438
439  return q;
440}
441
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
449void sqpipe_quantities_init(struct sqpipe_quantities *q) {
450  static struct quantities_vtable vtable = {
451    &sqpipe_quantities_get_quantity,
452    &sqpipe_quantities_update,
453    &sqpipe_quantities_destroy
454  };
455  q->super.vtable = &vtable;
456
457
458  q->a = NULL;
459  q->d = NULL;
460  q->w = NULL;
461  q->h = NULL;
462  q->u = NULL;
463  q->z = NULL; 
464  q->b = NULL;
465  q->t = NULL;
466}
467
468struct sqpipe_quantities* sqpipe_quantities_new() {
469  struct sqpipe_quantities *p = malloc(sizeof(struct sqpipe_quantities));
470
471  return p; 
472}
473
474
475struct quantities* sqpipe_quantities_destroy(struct quantities *qs) {
476  free(qs);
477  qs = NULL;
478 
479  return qs;
480}
481/**********************************************************/
482
483
484/**********************************************************
485Pipe domain class
486This is a virtual class and should never be created.
487The inheriting class needs to implement:
488  init, new, compute_fluxes
489**********************************************************/
490struct sqpipe_domain {
491  struct domain super;
492  long *state;
493};
494
495// Construction
496void sqpipe_domain_init(struct sqpipe_domain *d);
497struct sqpipe_domain* sqpipe_domain_new();
498
499// Destruction
500struct domain* sqpipe_domain_destroy(struct domain *d);
501
502// Methods
503struct cell* sqpipe_domain_get_cell(struct domain *D, int k, struct cell *c);
504
505// 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
545void sqpipe_domain_init(struct sqpipe_domain *d) {
546  static struct domain_vtable vtable = {
547    NULL,
548    NULL,
549    &domain_compute_fluxes_cells_generic,
550    &sqpipe_domain_get_cell,
551    &sqpipe_domain_destroy
552  };
553  d->super.vtable = &vtable;
554
555  ((struct domain *)d)->number_of_equations = 2;
556
557  ((struct domain *)d)->vertex_values = NULL;
558  ((struct domain *)d)->boundary_values = NULL;
559  ((struct domain *)d)->explicit_update = NULL;
560  ((struct domain *)d)->neighbours = NULL;
561  ((struct domain *)d)->neighbour_vertices = NULL;
562  ((struct domain *)d)->normals = NULL;
563  ((struct domain *)d)->areas = NULL;
564  ((struct domain *)d)->number_of_elements = 0; 
565
566  ((struct domain *)d)->timestep = 0.0;
567  ((struct domain *)d)->max_speed_array = NULL;
568
569  d->state = NULL;
570  // Can't set params to NULL here since sqpipe_domain_new() has allocated
571  // memory
572  // If a domain is created without using domain_new(), then params will be an
573  // unitialised pointer.
574
575  // Can't allocate memory here since init should not allocate memory
576  // The moral: Use domain_new() to create domains! 
577}
578
579struct sqpipe_domain* sqpipe_domain_new() {
580  struct sqpipe_domain *p = malloc(sizeof(struct sqpipe_domain));
581
582  ((struct domain *)p)->params = (struct params *) sqpipe_params_new();
583
584  return p;
585}
586
587
588struct domain* sqpipe_domain_destroy(struct domain *D) {
589  D->vertex_values = quantities_destroy(D->vertex_values);
590  D->boundary_values = quantities_destroy(D->boundary_values);
591  D->explicit_update = quantities_destroy(D->explicit_update);
592  D->params = params_destroy(D->params);
593  free(D);
594  D = NULL;
595
596  return D;
597}
598/**********************************************************/
599
600
Note: See TracBrowser for help on using the repository browser.