#include "utils.h"
// Maximum number of source terms
#define NSource 10
// Required from solver
// fields updated by time-integration
extern scalar * evolving;
// how to compute updates
double (* update) (scalar * evolving, scalar * updates, double dtmax) = NULL;
// User-provided functions
// gradient
double (* gradient) (double, double, double) = minmod2;
double t = 0., dt = 0.;
static void advance_generic (scalar * output, scalar * input, scalar * updates,
double dt){
if (input != output)
trash (output);
foreach() {
scalar o, i, u;
for (o,i,u in output,input,updates)
o[] = i[] + dt*u[];
}
boundary (output);
}
static void (* advance) (scalar * output, scalar * input, scalar * updates,
double dt) = advance_generic;
// how to compute source terms
void (* updatesource[NSource] ) (scalar * evolving,
scalar * sources,
double dtmax,
int numbersource);
// Fixing minor bug when a pointer is calling NULL
static void fnull(){
}
// Because of the additivity of source terms, we need to reset them at each step
static void resetsource(scalar * evolving,
scalar * sources,
double dtmax,
int numbersource){
scalar dsh = sources[0];
vector dshu = { sources[1], sources[2] };
foreach(){
dsh[] = 0;
foreach_dimension() dshu.x[] = 0;
}
numbersource++;
updatesource[numbersource](evolving, sources,
dtmax, numbersource);
}
event defaults (i = 0)
{
// limiting
for (scalar s in all)
s.gradient = gradient;
// default values
foreach()
for (scalar s in all)
s[] = 0.;
boundary (all);
}
void run()
{
t = 0.;
init_grid (N);
// main loop
timer start = timer_start();
int i = 0; long tnc = 0;
while (events (i, t, true)) {
// list of updates
scalar * updates = list_clone (evolving);
scalar * sources = list_clone (evolving);
dt = dtnext (t, update (evolving, updates, DT));
if (gradient != zero) {
/* 2nd-order time-integration */
scalar * predictor = list_clone (evolving);
/* predictor */
advance (predictor, evolving, updates, dt/2.);
updatesource[0](predictor, sources, dt/2.,0);
advance (predictor, predictor, sources, dt/2.);
/* corrector */
update (predictor, updates, dt);
delete (predictor);
free (predictor);
}
advance (evolving, evolving, updates, dt);
updatesource[0] (evolving, sources, dt,0);
advance (evolving, evolving, sources, dt);
delete (updates);
free (updates);
delete(sources);
free(sources);
long nc = 0;
foreach (reduction(+:nc))
nc++;
tnc += nc;
i++; t = tnext;
}
timer_print (start, i, tnc);
free_grid();
}
// Overloading the first updatesource function with resetsource
int numbersource;
event initres(i = 0 ){
numbersource = 0;
updatesource[numbersource]=resetsource;
numbersource++;
updatesource[numbersource] = fnull;
}