file: Nektar2d/src/analyser.c


/*---------------------------------------------------------------------------*
 *                        RCS Information                                    *
 *                                                                           *
 * $Source:
 * $Revision:
 * $Date:
 * $Author:
 * $State: 
 *---------------------------------------------------------------------------*/

#include <stdio.h>
#include <string.h>
#include "nektar.h"

/*
 * Run-time Analyzer ... called every time step
 */

static int  History       (Domain *omega, double time);        //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
static void hisHeader     (Domain *omega);        //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
static void addfields     (Element_List *V[]);
static void replacefields (Element_List *V[], int nsteps);
static int check_number=0;

int init = 1, verbose, iostep, hisstep, nsteps, timeavg;
static double last_time=0.0;

void Analyser (Domain *omega, int step, double time)     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
{
  FILE     *fp = 0;
  int      nfields = DIM+1;

  int i;
  Element_List  **V;
  char      fname[FILENAME_MAX];
  double    step_length;
  Element       *E, *F;


  dparam_set("t", time);

  V = (Element_List**) malloc(nfields*sizeof(Element_List*));
  if(init){
    verbose   = option("verbose");
    iostep    = iparam("IOSTEP");
    hisstep   = option("hisstep");
    if(!hisstep) hisstep = iparam("HISSTEP");
    nsteps    = iparam("NSTEPS");
    timeavg   = option("timeavg");
    if (omega->his_list) hisHeader (omega);
    init = 0;
  }
  /* .......... General Output ......... */
  step_length = (clock()-last_time)/(double)CLOCKS_PER_SEC;
  last_time = clock();
  fprintf(stdout, "Time step = %d, Time = %g Cpu-Time = %g\n", 
	 step, time, step_length);
  fflush(stdout);

  if (step == 0){                         /* Everything else is for step > 0 */
    if (omega->his_list)                  /* Do initial history point        */        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
      History (omega, time);
    return;
  }
  if ((step % hisstep) == 0){
    if (omega->his_list){
      History (omega, time);
      fflush(omega->his_file);
    }
    
#ifdef FORCES
    if (omega->fce_file){
      forces(omega,time);        //OVERLOAD CALL: forces: forces.c(?), forces.c(?)
      fflush(omega->fce_file);
    }   
#endif
 
    /* flush stdout at step as well */
    fflush(stdout);
  }
  

  /* ..........  Field Files   ......... */
  
  V[0] = omega->U;
  V[1] = omega->V;
  V[2] = omega->Uf;
  
  if (step % iostep == 0 && step < nsteps) {         
    int *size = ivector(0, MAXFACETS-1);
    
    /* project back pressure onto continous basis */
#ifndef PCONTBASE
    if((SLVTYPE) iparam("SLVTYPE") == StokesSlv)
      for(E=omega->P->fhead;E;E=E->next){
	E->Obwd(E->vert->hj,*E->h,E->dgL);        //OVERLOAD CALL: Obwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
	E->Jfwd(E);        //OVERLOAD CALL: Jfwd: Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex), Transforms.c(Element)
      }
#endif
    
    for(E=omega->P->fhead;E;E=E->next){
      for(i=0;i<E->Nedges;++i)
	size[i] = E->edge[i].l;
      size[i] = E->face[0].l;
      F = V[2]->flist[E->id];
      F->Copy_field(E->vert[0].hj, size);        //OVERLOAD CALL: Copy_field: Fieldfiles.c(Tri), Fieldfiles.c(Quad), Fieldfiles.c(Tet), Fieldfiles.c(Pyr), Fieldfiles.c(Prism), Fieldfiles.c(Hex), Fieldfiles.c(Element)
    }
    V[2]->Set_state('t');
    free(size);
    
    if (option ("checkpt")) {
      if(option("SLICES")&&check_number<100){
	sprintf (fname, "%s_%d.chk",  omega->name,check_number);
	++check_number;
      }
      else
	sprintf (fname, "%s.chk", omega->name);
      fp = fopen(fname,"w");
      Writefield (fp, omega->name, step, time, nfields, V);

      fclose(fp);
    }
  }

  if (step % hisstep == 0) { 
    if(verbose){
      V[2] = omega->P;
#if defined(WANNIER) || defined(WOMERSLEY)
      V[0]->Terror("u");
      V[1]->Terror("v");
#else
      if(omega->soln)
	for(i=0;i<nfields;++i)
	  V[i]->Terror(omega->soln[i]);
      else
	V[0]->Terror("0.0");
#endif
    }
  }
  
  if(timeavg){
    addfields (V);
    if(step == nsteps) replacefields(V,nsteps);
  }

  V[0]->Set_state('t');
  V[1]->Set_state('t');
  //  V[2]->Set_state('t');
  
  return;		
}
    

/* ------------------------------------------------------------------------- *
 * History() -- Process history points                                       *
 *                                                                           *
 * This function processes the history point list and transforms data points *        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
 * if necessary to physical space.                                           *
 * ------------------------------------------------------------------------- */
static int gatherPts(HisPoint *hp, Element_List *V[DIM+1],         //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?)
		     double *vbuf[HP_MAX]);        //OVERLOAD CALL: HP_MAX: nekcomp.h(?), nektar.h(?)

static int History (Domain *omega, double time){     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
  FILE     *fp = omega->his_file;
  HisPoint *hp = omega->his_list;        //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?)
  Element_List  *V [DIM+2];
  double   *vbuf[HP_MAX];        //OVERLOAD CALL: HP_MAX: nekcomp.h(?), nektar.h(?)
  register  int n, cnt;

  if (!hp) return 0;

  V[0] = omega->U; V[1] = omega->V;  V[2] = omega->P;

  gatherPts (hp, V, vbuf);
  cnt     = 0;
  
  do 
    { fprintf (fp, "%lf ", time);
      for (n = 0; n < strlen(hp->flags); n++)
	fprintf (fp, "%#13.6g ", vbuf[cnt][n]);
      fprintf (fp, ":%d\n", cnt+1);
      free (vbuf[cnt++]);    } 
  while 
    (hp = hp->next);
  
  return cnt;
}

/* Collect history points from the Elements */
static int vFe[][2] = {{0,1},{1,2},{2,0},{0,3},{1,3},{2,3}};
static double *modecenter(Element *E, Edge *e);

static int gatherPts (HisPoint *hp, Element_List **V, double *vbuf[HP_MAX])     //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?); HP_MAX: nekcomp.h(?), nektar.h(?)
{
  register  int i, j, n, pos;
  int cnt;

  cnt = DIM+1;

  for (i = 0; hp ; ++i, hp = hp->next) {
    vbuf[i] = dvector (0, (int) strlen (hp->flags)-1);
    for  (n = pos = 0; n < cnt; ++n) {
      if (strchr (hp->flags, V[n]->fhead->type)) {
	switch (hp->mode) {
	case TransVert:
	  if(V[n]->fhead->state == 't')
	    vbuf[i][pos++] = V[n]->flist[hp->id]->vert[hp->i].hj[0];
	  break;
	case TransEdge:
	  if(V[n]->fhead->state == 't'){
	    Edge   *e    = V[n]->flist[hp->id]->edge+hp->i;
	    double *mode = modecenter(V[n]->flist[hp->id],e);
	    
	    vbuf[i][pos] =0.5*(V[n]->flist[hp->id]->vert[vFe[hp->i][0]].hj[0]
		      +V[n]->flist[hp->id]->vert[vFe[hp->i][1]].hj[0]);

	    for(j = 0; j < e->l; ++j) vbuf[i][pos] += e->hj[j]*mode[j];
	    pos++;
	  }
	  break;
	case Physical:
	  break;
	default:
	  error_msg (History -- unknown history point mode);        //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?); point: Coords.c(?), Curvi.c(?)
	    break;
	}
      }
    }
  }
  
  return i;
}

/* find the center of modes  and store */
static double *modecenter(Element *E, Edge *edg){
  static double *mode;
  static int Lmode;

  if(!(mode&&(edg->l <= Lmode))){
    int i,qa = E->qa;
    Mode *e = E->getbasis()->edge[0];        //OVERLOAD CALL: getbasis: Basis.c(Tri), Basis.c(Quad), Basis.c(Tet), Basis.c(Pyr), Basis.c(Prism), Basis.c(Hex), Basis.c(Element)
    
    if(!mode) free(mode);
    mode = dvector(0,edg->l);

    if(qa%2 == 0){ /* if even spacing interpolate to center point */        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
      double **im;
      getim(qa,qa+1,&im,a2a);
      
      for(i = 0; i < edg->l; ++i)
	mode[i] = ddot(qa,im[qa/2],1,e[i].a,1);
    }      
    else           /* else use center value which is always at center point */        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
      for(i = 0; i < edg->l; ++i)
	mode[i] = e[i].a[qa/2];
  }

  return mode;
}

/* Write the header for the history point file */        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)

static void hisHeader (Domain *omega)     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
{
  FILE      *fp = omega->his_file;
  HisPoint  *hp = omega->his_list;        //OVERLOAD CALL: HisPoint: nektar.h(?), nekcomp.h(?)
  Element_List   *U  = omega->U;
  int        n  = 1;

  if (!fp) return;

  fputs ("# Nektar history point file\n"        //OVERLOAD CALL: point: Coords.c(?), Curvi.c(?)
	 "# \n"
	 "# History points:\n", fp);

  do 
    {
      if(hp->mode == TransVert){
	fprintf (fp, "#  %3d: x = %#6.2lf, y = %#6.2lf, ", n++, 
		 U->flist[hp->id]->vert[hp->i].x,
		 U->flist[hp->id]->vert[hp->i].y);
      }
      else if(hp->mode == TransEdge){
	fprintf (fp, "#  %3d: x = %#6.2lf, y = %#6.2lf, ", n++,
		 0.5*(U->flist[hp->id]->vert[vFe[hp->i][0]].x + 
		      U->flist[hp->id]->vert[vFe[hp->i][1]].x),
		 0.5*(U->flist[hp->id]->vert[vFe[hp->i][0]].y + 
		      U->flist[hp->id]->vert[vFe[hp->i][1]].y));
      }	

      fprintf (fp, "fields = %4s, [%2d %3d]\n",
	       hp->flags,  hp->i+1, hp->id+1);  
    }
  while
    (hp = hp->next);
  fputs ("#\n", fp);
  return;
}

static double **avg;

static void addfields (Element_List *V[]){
  register int i;
  Element  *E;
  double   *s;
  int       nf = DIM+1;

  if(!avg){
    int ntot = 0;
    /* count transform storage */
    for(E = V[0]->fhead; E; E = E->next)  ntot += E->Nmodes;
    
    avg = dmatrix(0,DIM,0,ntot-1);
    dzero (nf*ntot,avg[0],1);
  }

  for(i = 0; i < nf; ++i)
    for(E = V[i]->fhead, s = avg[i]; E; E = E->next){
      dvadd (E->Nmodes,E->vert->hj,1,s,1,s,1);
      s += E->Nmodes;
    }
}

static void replacefields (Element_List *V[], int nsteps){
  register int i;
  Element  *E;
  double   *s,fac;
  int       nf = DIM+1;

  fac = 1.0/(double)nsteps;
  
  for(i = 0; i < nf; ++i)
    for(E = V[i]->fhead, s = avg[i]; E; E = E->next){
      dsmul(E->Nmodes,fac,s,1,E->vert->hj,1);
      s += E->Nmodes;
    }

  free_dmatrix(avg,0,0);
}



Back to Source File Index


C++ to HTML Conversion by ctoohtml