Arman Akbarian
UNIVERSITY OF BRITISH COLUMBIA
PHYSICS & ASTRONOMY DEPT.

//=============================================================================
// application interface functions for wave example
//=============================================================================
// PAMR example to solve 2D wave equation
// See http://laplace.physics.ubc.ca/People/arman/file/wave2d_pamr.tar.gz for source

#include <stdlib.h>
#include <stdio.h>
#include <pamr.h>
#include <amrd.h>
#include <math.h>
#include <mpi.h>
#include "wave.h"

#include <sys/timeb.h>
void elapsed_time(void);

//=============================================================================
// id parameters (time symmetric gaussian)
//=============================================================================

//Initialization parameters:
real amp, tdam, xc, xwid, yc, ywid;
//=====================================

//=============================================================================
// some convenient, "local" global variables
//=============================================================================

//The pointers to actual grid function data:
real *f_t_n, *f_t_np1, *f_n, *f_np1;
//=======================================

//The pointers to grid:
real *x, *y;
//===========================

//Parameters that describe the grid:
int shape[2]; // number of grid points in each dim, basically shape[0]=Nx, shape[1]=Ny
int ghost_width[4]; //side of ghost cells, 2*dim array, at x_min, x_max, y_min and so on
int Nx,Ny;
int phys_bdy[4]; // 1 or 0, wether or not boundary is physical, has size = dim*2
int lb_flag,rb_flag,tb_flag,bb_flag; //convenient variable for phys boundary indication
int size; //size=Nx*Ny*...
int dim; //the spatial dimension of the grid (only 1,2 or 3)
int g_rank; // To store the MPI rank of the execution

real base_bbox[4], bbox[4], dx, dy, dt;
// bbox stores the coordinate bounding box [x1,x2,y1...]
// base_bbox is the coordinate of the entire domain
int g_L; //to store the grid level (in AMR heirarchy)

//To store the grid function numbers, see set_gfns subroutine
int f_t_n_gfn, f_t_np1_gfn, f_n_gfn, f_np1_gfn;
//===========

//=============================================================================
// call after variables have been defined
//=============================================================================
void set_gfns(void)
{
if ((f_t_np1_gfn = PAMR_get_gfn("f_t",PAMR_AMRH,1))<0) AMRD_stop("set_gnfs error",0);
if ((f_t_n_gfn   = PAMR_get_gfn("f_t",PAMR_AMRH,2))<0) AMRD_stop("set_gnfs error",0);

if ((f_np1_gfn   = PAMR_get_gfn("f",PAMR_AMRH,1))<0) AMRD_stop("set_gnfs error",0);
if ((f_n_gfn     = PAMR_get_gfn("f",PAMR_AMRH,2))<0) AMRD_stop("set_gnfs error",0);

// get_gfn returns the grid function number:
// PAMR_get_gfn(A,B,C)
// A >> grid function name
// B >> in which heirarchy, two options: PAMR_AMRH: AMR heirarchy, PAMR_MGH: multigrid h.
// C >> time level if it is in AMRH, else 0, latest time ordered first.
}

//=============================================================================
// call with valid iter to set up globals:
//=============================================================================
void ldptr(void)
{
real dx0[2];
real *x0[2]; //array of size dim, pointing to local perimeter coordinate array
real *gfs[PAMR_MAX_GFNS]; //array of size ngfs, pointers to local data array
static int first=1;

if (first)
{
first=0;
set_gfns();
PAMR_get_global_bbox(base_bbox); //AAK??
}
//These are to get the information of the grid
PAMR_get_g_dim(&dim); //saves the dimension of the grid in dim
PAMR_get_g_rank(&g_rank); // gets the MPI rank
PAMR_get_g_shape(shape); //gets the shape (Nx,Ny...)
PAMR_get_g_bbox(bbox); //gets the boundary (x1,x2,y1,y2...)
PAMR_get_g_ghost_width(ghost_width); //gets the ghots width at boundaries
PAMR_get_g_level(&g_L); // gets the level number of the current grid
PAMR_get_dxdt(g_L,dx0,&dt); //clear
dx=dx0[0]; //convenient variables??
dy=dx0[1];

//AAK: This is for finding the real boundary
if ((bbox[0]-base_bbox[0])<dx/2) phys_bdy[0]=1; else phys_bdy[0]=0;
if ((base_bbox[1]-bbox[1])<dx/2) phys_bdy[1]=1; else phys_bdy[1]=0;
Nx=shape[0];

if ((bbox[2]-base_bbox[2])<dy/2) phys_bdy[2]=1; else phys_bdy[2]=0;
if ((base_bbox[3]-bbox[3])<dy/2) phys_bdy[3]=1; else phys_bdy[3]=0;
Ny=shape[1];
size=Nx*Ny;

//setting the convenient variables
lb_flag = phys_bdy[0];
rb_flag = phys_bdy[1];
bb_flag = phys_bdy[2];
tb_flag = phys_bdy[3];
//=========================
PAMR_get_g_x(x0); //gets the pointer to the coordinate array

x=x0[0]; //the pointers to the array of actual grid data (size Nx*Ny..)
y=x0[1];

PAMR_get_g_gfs(gfs); //gets the pointer to local data of grid function

//AAK:store the pointer in convenient variables
f_n     = gfs[f_n_gfn-1]; //AAK?? what is -1, numbering starts from  1?
f_np1   = gfs[f_np1_gfn-1];
f_t_n   = gfs[f_t_n_gfn-1];
f_t_np1 = gfs[f_t_np1_gfn-1];

}

//=============================================================================
// utility routines
//=============================================================================
void const_f(real *f, real c)
{
int i;

for (i=0; i<Nx*Ny; i++) f[i]=c;
}

void zero(real *f)
{
const_f(f,0);
}

//=============================================================================
// Routines required by amrd:
//=============================================================================

//=============================================================================
// Returns 0 to use default mechanism, or is expected to calculate
// the correct initial hierarchy and return 1:
//=============================================================================
int wave_id(void)
{
if( my_rank == 0 ) elapsed_time();
return 0;
}

//=============================================================================
// Sets custom parameters, variables, etc. Split up into two segments,
// one called before the pamr context is initialized and standard
// parameters are read, and the other afterwards
//=============================================================================
void wave_var_pre_init(char *pfile)
{
return;
}

//AAK: one of the hook functions passed to amrd, to read the parameters from file
void wave_var_post_init(char *pfile)
{
int i,j;
char buf[64]; //AAK??

if (my_rank==0)
{
system("date > date.dat");
printf("===================================================================\n");
}

amp=tdam=xc=yc=0;
xwid=ywid=1;
//AAK: the file is in RNPL format
// last argument is the size of the array, 1 for just a real number
AMRD_real_param(pfile,"amp",&amp,1);
AMRD_real_param(pfile,"tdam",&tdam,1);
AMRD_real_param(pfile,"xc",&xc,1);
AMRD_real_param(pfile,"yc",&yc,1);
AMRD_real_param(pfile,"xwid",&xwid,1);
AMRD_real_param(pfile,"ywid",&ywid,1);

if (my_rank==0) printf("===================================================================\n");
return;
}

//=============================================================================
// Sets all t=n variables to their 'zero' values:
//=============================================================================
void wave_AMRH_var_clear(void)
{
ldptr();

//   zero(phi_n);

zero(f_n);
zero(f_t_n);
return;
}

//=============================================================================
//
// currently, we only allow for time-symmetric initial data
//=============================================================================
void wave_free_data(void)
{
ldptr();
//This initializes the grid, all the calculations are encapsulated in the routine.
initializer0_(f_t_n,f_n,&Nx,&Ny,x,y,&amp,&tdam,&xc,&xwid,&yc,&ywid);

return;
}

//=============================================================================
// Initial constraint data --- called after each MG iteration.
//=============================================================================
void wave_t0_cnst_data(void)
{
return;
}

//=============================================================================
// Calculations prior to saving info to disk.
//
// NOTE: at this point, the time sequence is: n,nm1,np1
//=============================================================================
void wave_pre_io_calc(void)
{
return;
}

//=============================================================================
// Returns some norm of the residual for the evolution variables,
// called after an evolution iteration.
// We're using an explicit scheme to solve for phi, hence return 0
//=============================================================================
real wave_evo_residual(void)
{
return 0;
}

//=============================================================================
// Returns some norm of the residual for the MG variables, *AND*
// stores the point-wise residual in "f_res" for each MG variable "f" (for
// computing new RHS's)
//=============================================================================
real wave_MG_residual(void)
{
return 0;
}

//=============================================================================
// Performs 1 iteration of the evolution equations
//=============================================================================
void wave_evolve(int iter)
{
ldptr();
//This updates the grid, all the calculations are encapsulated in the routine
#ifdef __L
update0__(f_t_np1,f_t_n,f_np1,f_n,&Nx,&Ny,&dt,&dx,&dy,&bb_flag,&lb_flag,&rb_flag,&tb_flag);
#else
update0_(f_t_np1,f_t_n,f_np1,f_n,&Nx,&Ny,&dt,&dx,&dy,&bb_flag,&lb_flag,&rb_flag,&tb_flag);
#endif

return;
}

//=============================================================================
//=============================================================================
{
}

//=============================================================================
void wave_fill_bh_bboxes(real *bbox, int *num, int max_num)
{
}

//=============================================================================
void wave_post_tstep(int L)
{
return;
}

//=============================================================================
// Performs 1 relaxation sweep of the MG equations, and returns an estimate
// of the norm of the residual.
//=============================================================================
real wave_MG_relax(void)
{
return 0;
}

//=============================================================================
// Computes the differential operator L acting on the current grid,
// in a region specified by the grid function "mask". Stores the result
// in "f_lop" for each MG variable "f"
//=============================================================================
void wave_L_op(void)
{
return;
}

//=============================================================================
// Called after calculating the TRE for all variables
//=============================================================================
void wave_scale_tre(void)
{
return;
}

//=============================================================================
// post-regrid initialization of constant functions
//=============================================================================
void wave_post_regrid(void)
{
}

//=============================================================================
int main(int argc, char **argv)
{
amrd(argc,argv,&wave_id,&wave_var_pre_init,
&wave_var_post_init, &wave_AMRH_var_clear,
&wave_free_data, &wave_t0_cnst_data,
&wave_evo_residual, &wave_MG_residual,
&wave_evolve, &wave_MG_relax, &wave_L_op,
&wave_pre_io_calc, &wave_scale_tre,
&wave_post_regrid, &wave_post_tstep,
if (my_rank==0) elapsed_time();
}

//=============================================================================
// Maintains/reports elapsed wall-clock time.
//=============================================================================
void elapsed_time(void) {
static int    first = 1;
struct        timeb t;
static double msinit;
double        mscurr, mselapsed;

ftime(&t);
mscurr = 1000.0 * t.time + t.millitm;
if( first ) {
msinit = mscurr;
first = 0;
}
mselapsed = mscurr - msinit;
printf("elapsed_time: Seconds since initial call: %12.3f\n",
mselapsed / 1000.0);
}

last update: Wed May 11, 2016