source: trunk/anuga_work/development/2010-projects/anuga_1d/utilities/solver.h @ 8276

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

Should have been included in revision 8274

File size: 12.3 KB
Line 
1/**********************************************************
2Generic params class
3**********************************************************/
4struct params;
5
6// vtable
7struct params_vtable { 
8  struct params* (*destroy) (struct params *q);
9};
10
11// structure
12struct params {
13  struct params_vtable *vtable;
14 
15  double epsilon;
16  double cfl;
17};
18
19// Public API
20
21// Destruction
22struct params* params_destroy(struct params *q);
23
24// method declarations
25
26// method implementations
27struct params* params_destroy(struct params *q) {
28  return q->vtable->destroy(q);
29}
30/**********************************************************/
31
32
33/**********************************************************
34Generic quantity class
35**********************************************************/
36struct quantity;
37
38// vtable
39struct quantity_vtable {
40  double* (*flux_formula) (struct quantity *q, double normal, struct params *p, double *quantityflux);
41  double (*sound_speed) (struct quantity *q, struct params *p);
42  double (*get_conserved) (struct quantity *q, int k, double normal);
43  struct quantity* (*destroy) (struct quantity *q);
44};
45
46// structure
47struct quantity {
48  struct quantity_vtable *vtable;
49};
50
51// Public API
52// Destruction
53struct quantity* quantity_destroy(struct quantity *q);
54
55// method declarations
56double* quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux);
57double quantity_sound_speed (struct quantity *q, struct params *p);
58double quantity_get_conserved (struct quantity *q, int k, double normal);
59
60// method implementations
61double* quantity_flux_formula (struct quantity *q, double normal, struct params *p, double *quantityflux) {
62  return q->vtable->flux_formula(q, normal, p, quantityflux);
63}
64
65double quantity_sound_speed (struct quantity *q, struct params *p) {
66  return q->vtable->sound_speed(q, p);
67}
68
69double quantity_get_conserved (struct quantity *q, int k, double normal) {
70  return q->vtable->get_conserved(q, k, normal);
71}
72
73struct quantity* quantity_destroy(struct quantity *q) {
74  return q->vtable->destroy(q);
75}
76/**********************************************************/
77
78
79/**********************************************************
80Generic edge class
81**********************************************************/
82struct edge;
83
84// vtable
85struct edge_vtable{
86  double* (*flux_function) (struct edge *e, struct params *p, double *edgeflux, int num_eqns);
87  double (*extreme_sound_speeds) (struct edge *e, struct params *p);
88  struct edge* (*destroy) (struct edge *e);
89};
90
91// structure
92struct edge {
93  struct edge_vtable *vtable;
94  struct quantity *qin, *qout;
95  double normal;
96  double smin, smax;
97};
98
99// Public API
100// Destruction
101struct edge* edge_destroy(struct edge *e);
102
103// method declarations
104double* edge_flux_function (struct edge *e, struct params *p, double *edgeflux, int num_eqns);
105double edge_extreme_sound_speeds (struct edge *e, struct params *p);
106
107// method implementations
108double* edge_flux_function (struct edge *e, struct params *p, double *edgeflux, int num_eqns) {
109  return e->vtable->flux_function(e, p, edgeflux, num_eqns);
110}
111
112double edge_extreme_sound_speeds (struct edge *e, struct params *p) {
113  return e->vtable->extreme_sound_speeds(e, p);
114}
115
116struct edge* edge_destroy(struct edge *e) {
117  return e->vtable->destroy(e);
118}
119/**********************************************************/
120
121
122/**********************************************************
123Generic cell class
124**********************************************************/
125struct cell;
126
127// vtable
128struct cell_vtable{
129  double* (*flux_function) (struct cell *c, struct params *p, double *cellflux, int num_eqns);
130  double (*extreme_sound_speeds) (struct cell *c, struct params *p);
131  double* (*forcing_terms) (struct cell *c, struct params *p, double *cellforce);
132  struct cell* (*destroy) (struct cell *c);
133};
134
135// structure
136struct cell {
137  struct cell_vtable *vtable;
138  struct edge **edges;
139  int num_edges;
140};
141
142// Public API
143
144// Destruction
145struct cell* cell_destroy(struct cell *c);
146
147// method declarations
148double* cell_flux_function (struct cell *c, struct params *p, double *cellflux, int num_eqns);
149double cell_extreme_sound_speeds (struct cell *c, struct params *p);
150double* cell_forcing_terms (struct cell *c, struct params *p, double *cellforce);
151
152// method implementations
153double* cell_flux_function (struct cell *c, struct params *p, double *cellflux, int num_eqns) {
154  return c->vtable->flux_function(c, p, cellflux, num_eqns);
155}
156
157double cell_extreme_sound_speeds (struct cell *c, struct params *p) {
158  return c->vtable->extreme_sound_speeds(c, p);
159}
160
161double* cell_forcing_terms (struct cell *c, struct params *p, double *cellforce) {
162  return c->vtable->forcing_terms(c, p, cellforce);
163}
164
165struct cell* cell_destroy(struct cell *c) {
166  return c->vtable->destroy(c);
167}
168/**********************************************************/
169
170
171/**********************************************************
172Generic quantities class
173**********************************************************/
174struct quantities;
175
176// vtable
177struct quantities_vtable{
178  struct quantity* (*get_quantity) (struct quantities *qs, int i, struct quantity* q);
179  void (*update) (struct quantities *qs, double *flux, int k);
180  struct quantities* (*destroy) (struct quantities *qs);
181};
182
183// structure
184struct quantities {
185  struct quantities_vtable *vtable;
186};
187
188// Public API
189
190// Destruction
191struct quantities* quantities_destroy(struct quantities *qs);
192
193// method declarations
194struct quantity* quantities_get_quantity(struct quantities *qs, int i, struct quantity *q); 
195void quantities_update (struct quantities *qs, double *flux, int k);
196
197// method implementations
198struct quantity* quantities_get_quantity(struct quantities *qs, int i, struct quantity *q) {
199  return qs->vtable->get_quantity(qs, i, q);
200}
201
202void quantities_update (struct quantities *qs, double *flux, int k) {
203  return qs->vtable->update(qs, flux, k);
204}
205
206struct quantities* quantities_destroy(struct quantities *qs) {
207  return qs->vtable->destroy(qs);
208}
209/**********************************************************/
210
211/**********************************************************
212Generic domain class
213**********************************************************/
214struct domain;
215
216// vtable
217struct domain_vtable {
218  struct edge* (*new_edge) (struct domain *D);
219  struct cell* (*new_cell) (struct domain *D);
220  double (*compute_fluxes) (struct domain *D);
221  struct cell* (*get_cell) (struct domain *d, int i, struct cell *c); 
222  struct domain* (*destroy) (struct domain *D);
223};
224
225// structure
226struct domain {
227  struct domain_vtable *vtable;
228
229  struct quantities *vertex_values;
230  struct quantities *boundary_values;
231  struct quantities *explicit_update;
232  long* neighbours;
233  long* neighbour_vertices;
234  double* normals;
235  double* areas;
236  int number_of_elements;
237  int number_of_equations;
238  double *max_speed_array;
239  double timestep;
240  struct params *params;
241};
242
243// Public API
244// Destruction
245struct domain* domain_destroy(struct domain* D);
246
247// method declarations
248struct edge* domain_new_edge(struct domain *qs);
249struct cell* domain_new_cell(struct domain *qs);
250double domain_compute_fluxes (struct domain *D);
251struct cell* domain_get_cell(struct domain *d, int i, struct cell *c); 
252
253// method implementations
254struct edge* domain_new_edge(struct domain *qs) {
255  return qs->vtable->new_edge(qs);
256}
257
258struct cell* domain_new_cell(struct domain *qs) {
259  return qs->vtable->new_cell(qs);
260}
261
262double domain_compute_fluxes (struct domain *D) {
263  return D->vtable->compute_fluxes(D);
264}
265
266struct cell* domain_get_cell(struct domain *D, int i, struct cell *c) {
267  return D->vtable->get_cell(D, i, c);
268}
269
270struct domain* domain_destroy(struct domain* D) {
271  return D->vtable->destroy(D);
272}
273/**********************************************************/
274
275
276/**********************************************************
277Generic methods
278**********************************************************/
279// Declarations
280double domain_compute_fluxes_cells_generic(struct domain *D);
281double domain_compute_fluxes_edges_generic(struct domain *D);
282struct cell* sqpipe_domain_get_cell(struct domain *D, int i, struct cell *c);
283double* edge_flux_function_godunov (struct edge *c, struct params *p, double *edgeflux, int num_eqns);
284double* cell_flux_function_generic (struct cell *c, struct params *p, double *cellflux, int num_eqns);
285
286// Implementations
287
288// Computational function for flux computation
289// This iterates over cells.
290double domain_compute_fluxes_cells_generic (struct domain *D) {
291  double flux[D->number_of_equations];
292  double max_speed;
293  int k, i;
294  struct cell *c;
295   
296  // pre-create a cell object to use in the loop
297  // so we need only allocate memory once
298  c = domain_new_cell(D);
299
300  for (k=0; k<D->number_of_elements; k++) {
301    c = domain_get_cell(D, k, c);
302   
303    // Compute sound speeds (this updates c) and update timestep
304    max_speed = cell_extreme_sound_speeds(c, D->params);
305    if (max_speed > D->params->epsilon) {                                   
306      D->timestep = min(D->timestep, 0.5*D->params->cfl * D->areas[k]/max_speed); 
307    }
308
309    // Compute flux
310    cell_flux_function(c, D->params, flux, D->number_of_equations);   
311    for (i=0; i<D->number_of_equations; i++) {
312      flux[i] /= D->areas[k];
313    }
314    quantities_update(D->explicit_update, flux, k);
315
316    //Keep track of maximal speeds
317    D->max_speed_array[k]=max_speed;
318  }
319 
320  cell_destroy(c);
321
322  return D->timestep;   
323}
324
325// Computational function for flux computation
326// This iterates over edges.
327double domain_compute_fluxes_edges_generic (struct domain *D) {
328  double flux[D->number_of_equations];
329  double max_speed;
330  int k, i;
331  // pre-create a cell object to use in the loop
332  // so we need only allocate memory once
333  struct cell *c = domain_new_cell(D);
334
335  for (k=0; k<D->number_of_elements; k++) {
336    c = domain_get_cell(D, k, c);
337   
338    // Compute sound speeds (this updates c) and update timestep
339    max_speed = cell_extreme_sound_speeds(c, D->params);
340    if (max_speed > D->params->epsilon) {                                   
341      D->timestep = min(D->timestep, 0.5*D->params->cfl * D->areas[k]/max_speed); 
342    }
343
344    // Compute flux and store in explicit_update
345    cell_flux_function(c, D->params, flux, D->number_of_equations);   
346    for (i=0; i<D->number_of_equations; i++) {
347      flux[i] /= D->areas[k];
348    }
349    quantities_update(D->explicit_update, flux, k);
350
351    //Keep track of maximal speeds
352    D->max_speed_array[k]=max_speed;
353  }
354 
355  cell_destroy(c);
356
357  return D->timestep;   
358}
359
360
361struct cell* domain_get_cell_generic(struct domain *D, int k, struct cell *c) {
362  int i, ki, n, m, nm;
363
364  for (i=0; i<c->num_edges; i++) {
365    ki = k*(c->num_edges) + i;
366    c->edges[i]->normal = D->normals[ki];     
367    n = D->neighbours[ki];
368
369    c->edges[i]->qin = quantities_get_quantity(D->vertex_values, ki, c->edges[i]->qin);
370
371    if (n<0) {
372      m = -n-1;
373      c->edges[i]->qout = quantities_get_quantity(D->boundary_values, m, c->edges[i]->qout);
374    } else {
375      m = D->neighbour_vertices[ki];
376      nm = n*2+m;
377      c->edges[i]->qout = quantities_get_quantity(D->vertex_values, nm, c->edges[i]->qout);
378    }
379  } 
380
381  return c;
382}
383
384
385// Generic Godunov flux computation method
386double* edge_flux_function_godunov (struct edge *e, struct params *p, double *edgeflux, int num_eqns) {
387  int i;
388  double flux_in[num_eqns], flux_out[num_eqns];
389  double denom;
390
391  // Flux contribution from inside of  edge
392  quantity_flux_formula(e->qin, e->normal, p, flux_in);
393
394  // Flux contribution from outside of edge
395  quantity_flux_formula(e->qout, e->normal, p, flux_out);
396
397  // Flux computation
398  denom = e->smax - e->smin;
399  if (denom < p->epsilon) {
400    for (i=0; i<num_eqns; i++) edgeflux[i] = 0.0;
401  } else {
402    for (i=0; i<num_eqns; i++) {
403      edgeflux[i] = e->smax*flux_in[i] - e->smin*flux_out[i];
404      edgeflux[i] += e->smax*e->smin * (quantity_get_conserved(e->qout, i, e->normal) - quantity_get_conserved(e->qin, i, e->normal));
405      edgeflux[i] /= denom;
406    }
407  }
408
409  return edgeflux;
410}
411
412double* cell_flux_function_generic (struct cell *c, struct params *p, double *cellflux, int num_eqns) {
413  int i, j;
414  double cellforce[num_eqns];
415  double flux[num_eqns];
416
417  // Compute forcing terms first
418  cell_forcing_terms(c, p, cellforce);
419  for (i=0; i<num_eqns; i++) {
420    cellflux[i] = -cellforce[i];
421  }
422
423  // Compute flux through each edge and append to flux
424  for (j=0; j<c->num_edges; j++) {
425    edge_flux_function(c->edges[j], p, flux, num_eqns);
426    for (i=0; i<num_eqns; i++) {
427      cellflux[i] -= flux[i];
428    }
429  }
430
431  return cellflux;
432}
433/**********************************************************/
Note: See TracBrowser for help on using the repository browser.