sandbox/vheusinkveld/myfunctions/equidistant/equidistant_time.c
Equidistant averaging in time
Testing the equidistant averaging functions for convergence in time. First we define a function, cubed in space and squared in time. Moving at a speed of 1/tau. De time integration of this function is define afterwards. The choice for a more general function would be better, but this is just to see if the functions actually do what they say. Convergence in space and time will be very specific to the problem in question.
#define func(s,d,f) (cube(s-sq(t/tau)) + cube(d) + cube(f))
#define Ifunc(s,d,f) (cube(d) + cube(f) + cube(s) - sq(s)*sq(t)/sq(tau) + 3*s*pow(t, 4)/5./pow(tau,4) - pow(t,6)/(7*pow(tau,6)))
double tau;
double ave_err; // Diagnosing total error
#include "grid/octree.h" // 3D
#include "run.h"
int diagi = 0;
int diaglvl = 4; // Level at which we want equidistant diagnostisation
int minlevel = 4;
int maxlevel = 4; // set equal to diag level such that spacial error = 0
scalar b[];
double *equifield = NULL; // The diagnostics field
Used functions which are define at the end of the document
int equi_diag (scalar s, int lvl, int diagii);
void equi_output (scalar s, FILE * fp, int lvl, int diagii);
int main(){
FILE * fp3 = fopen("output_time", "w"); // Specify te error output file
L0=1.;
int pp = 0;
Loop over different DT to see if average converges
for(DT = 0.1; DT > 0.001; DT = DT/2.){
init_grid(1<<maxlevel);
tau = 0.01; // 'fast' moving function
run();
if(pid() == 0) {
if(pp==0){
fprintf(fp3, "err\tdt\tL0\n");
}
fprintf(fp3, "%g\t%g\t%g\n", ave_err, DT, L0); // Write away error
}
free(equifield); // we dont want memory leaks
equifield = NULL; // reset the field pointer
diagi=0;
pp++;
}
}
Ability to test for non uniform grid
Set values of scalar b based on defined function
Keeping track of the field b (averaging), here we cheat a bit since the error determination is defined in the function itself.
Write away output of the averaging function
event end(t=1){
FILE * fp2 = fopen("output", "w");
equi_output(b, fp2, diaglvl, diagi);
fclose(fp2);
}
Results
The average error seems to nicely converge in a linear fasion.
set xr [0.:0.11]
set yr [0.:0.3]
set xlabel 'DT'
set ylabel 'Error
plot 'output_time' using 2:1 title 'DT vs err'
Used functions
NOTE: Very much in progress.
These functions are geared towards equidistant diagnostics. The level for which diagnostics are done can be specified. Values will be linearly interpolated if necesarry.
So far this function works with MPI and only works in 3D.
Short description of what they do:
//void equi_diag(scalar s){ adds values to the equifield array }
//void equi_output(scalar s, FILE * fp){ outputs equifield to file fp in the format: place = xnn + yn + zn }
TODO implement 2D TODO walking average since values in equifield get big, or divide before summing?
int equi_diag (scalar s, int lvl, int diagii){
int len = 1; // TODO implement multiple scalar field averaging
int n = 1<<lvl;
if(!equifield) {
equifield = calloc(len*n*n*n, sizeof(double)); // Init with zeros
}
double dDelta = 0.9999999*L0/n;
restriction({s}); // Make sure that coarse level values exist
boundary({s}); // Update boundaries
// Iterate over diag level of leaf cell
foreach_level_or_leaf(lvl) {
if(level < lvl){
int pn = 1<<(lvl-level);
for(int i = 0; i < pn; i++){
for(int j = 0; j < pn; j++){
for(int k = 0; k < pn; k++){
double ctrans = dDelta/2. - Delta/2.; // Translation to corner of cell
double xx = x + ctrans + i*dDelta; // Coordinates in equidistant grid
double yy = y + ctrans + j*dDelta;
double zz = z + ctrans + k*dDelta;
int ii = round((xx - dDelta/2.)/dDelta); // Iteration in equidistant grid
int jj = round((yy - dDelta/2.)/dDelta);
int kk = round((zz - dDelta/2.)/dDelta);
int place = ii*n*n + jj*n + kk*len; // Location in equidistant grid
double temp = interpolate(s, xx, yy, zz); // Linear interpolation
equifield[place] += temp;
}
}
}
} else {
int ii = round((x - dDelta/2.)/dDelta); // x adaptive is x equidistant
int jj = round((y - dDelta/2.)/dDelta); // Derive place in equidistant grid
int kk = round((z - dDelta/2.)/dDelta);
int place = ii*n*n + jj*n + kk*len;
equifield[place] += s[];
}
}
return diagii + 1;
}
void equi_output (scalar s, FILE * fp, int lvl, int diagii){
int len = 1;
int n = 1<<lvl;
double dDelta = 0.9999999*L0/n;
if(pid() == 0){
// Make sure all thread values are summed up using MPI
@if _MPI
MPI_Reduce (MPI_IN_PLACE, equifield, len*n*n*n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
@endif
double temperr = 0;
int len = 1;
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
for(int k = 0; k < n; k++){
double x = (i + 0.5)*dDelta; // Coordinates in equidistant grid
double y = (j + 0.5)*dDelta;
double z = (k + 0.5)*dDelta;
int place = i*n*n + j*n + k*len;
double temp = (equifield[place]/((double)diagii)); // Divide by number of additions
double val = Ifunc(x,y,z);
fprintf(fp, "%g\t", temp);
temperr += fabs((temp - val)/val);
}
}
}
ave_err = temperr/(cube(n));
}
@if _MPI
else // slave
MPI_Reduce (equifield, NULL, len*n*n*n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
@endif
}