file: Nektar2d/src/drive.c


/*---------------------------------------------------------------------------*
 *                        RCS Information                                    *
 *                                                                           *
 * $Source: 
 * $Revision:
 * $Date: 
 * $Author:
 * $State: 
 *---------------------------------------------------------------------------*/
// tcew done
#include <stdio.h>
#include <time.h>
#include "nektar.h"

static int    Je;            /* Externals */
static double dt, Re;

void MakeF   (Domain *omega, ACTION act, SLVTYPE slvtype);        //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?); ACTION: nekcomp.h(?), nektar.h(?)
static void StartUp (Domain *);        //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)

void switch_hj(Element_List *EL, double **hj);
void Bsystem_mem_free(Bsystem *Ubsys, Element_List *U);
void P_solve(Element_List *U, Element_List *Uf,Bndry *Ubc,Bsystem *Ubsys,SolveType Stype); 
void V_solve(Element_List *U, Element_List *Uf,Bndry *Ubc,Bsystem *Ubsys, double *us, SolveType Stype,int step);
void S_solve(Domain *Omega);        //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)

#if 1
#define Timing(s) \
 fprintf(stdout,"%s Took %g seconds\n",s,(clock()-st)/cps); \        //OVERLOAD CALL: cps: nekcomp.h(?), Solve_Stokes.c(?)
 st = clock();
#else
#define Timing(s) \
 /* Nothing */
#endif

int Nek_Fullcg = 0;

main(int argc, char *argv[]){
  Domain      *Omega;             /* Solution Domain            */        //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?); Domain: nekcomp.h(?), nektar.h(?)
  int          step, nsteps;      /* Discrete time (step)       */
  double       time = 0.0;        /* Continuous time            */
  ACTION       WaveProp;          /* Convective form            */        //OVERLOAD CALL: ACTION: nekcomp.h(?), nektar.h(?)
  SLVTYPE      SolveType;         /* Solver type                */
  double       begin_clock;
  double       st, cps = (double)CLOCKS_PER_SEC;        //OVERLOAD CALL: cps: nekcomp.h(?), Solve_Stokes.c(?)

#ifdef DEBUG
/*  mallopt(M_DEBUG,1); */
  init_debug();
#endif
  
  Omega     = PreProcess(argc,argv);


  Nek_Fullcg = option("FullCG");
  time      = dparam("STARTIME"); 
  nsteps    = iparam("NSTEPS");
  WaveProp  = (ACTION) iparam("EQTYPE");        //OVERLOAD CALL: ACTION: nekcomp.h(?), nektar.h(?)
  SolveType = (SLVTYPE) iparam("SLVTYPE");
  Je        = iparam("INTYPE");
  dt        = dparam("DELT");
  step      = 0;
  Re        = 1./dparam("KINVIS");   

  StartUp (Omega);

#if defined(WANNIER)||defined(WOMERSLEY)
  dparam_set("TIME",time);
  Omega->U->Terror("u"); 
  Omega->V->Terror("v"); 
#else
  Omega->U->Terror(Omega->soln[0]); 
  Omega->V->Terror(Omega->soln[1]);
  Omega->P->Terror(Omega->soln[2]);
#endif


  /*------------------------------------------------------------------*
   *                    Main Integration Loop                         *
   *------------------------------------------------------------------*/
  begin_clock = (double) clock();
  switch(SolveType){
  case Splitting:
    while (step < nsteps) {
      
      set_order(((step+1 < Je)? step+1: Je));
      
      st = clock();
      MakeF (Omega, Prep, SolveType);
      Timing("Prep.......");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      MakeF (Omega, WaveProp, SolveType);
      Timing("WavePropo..");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      Integrate_SS(Je, dt, Omega->U, Omega->Uf, Omega->u, Omega->uf);
      Integrate_SS(Je, dt, Omega->V, Omega->Vf, Omega->v, Omega->vf);
      Timing("Integrate..");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      MakeF (Omega, Pressure, SolveType);
      Timing("Pressure...");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      P_solve (Omega->P,Omega->Pf,Omega->Pbc,Omega->Pressure_sys,Helm);
      Timing("P_solve....");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      MakeF (Omega, Viscous, SolveType);
      Timing("Viscous....");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      V_solve(Omega->U,Omega->Uf,Omega->Ubc,Omega->Usys,
	      Omega->us[0],Helm,step);
      V_solve(Omega->V,Omega->Vf,Omega->Vbc,Omega->Vsys,
	      Omega->vs[0],Helm,step);
      Timing("V_solve....");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
    
    
      MakeF (Omega, Post, SolveType);
      Timing("Post.......");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
    
#ifdef WOMERSLEY
      dparam_set("TIME",time+dt);
#endif
      Analyser(Omega, ++step, (time += dt));
      Timing("Analyser...");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
    }
    break;
  case StokesSlv:
    while (step < nsteps) {
    
      set_order_CNAB(((step+1 < Je)? step+1: Je));
    
      st = clock();
      MakeF (Omega, Prep, SolveType);
      Timing("Prep.......");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      MakeF (Omega, WaveProp, SolveType);
      Timing("WavePropo..");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
      if((WaveProp != StokesS)&&(WaveProp != Oseen)){
	Integrate_CNAB(Je, dt, Omega->U, Omega->Uf, Omega->u, Omega->uf);
	Integrate_CNAB(Je, dt, Omega->V, Omega->Vf, Omega->v, Omega->vf);
	Timing("Integrate..");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
	
	MakeF(Omega, StokesForce, SolveType);
	Timing("StokesFce..");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      }

      S_solve(Omega);
      Timing("StokesSlv..");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
    
      if((WaveProp != StokesS)&&(WaveProp != Oseen)){
	MakeF (Omega, Post, SolveType);
	Timing("Post.......");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      }

#ifdef WOMERSLEY
      dparam_set("TIME",time+dt);
#endif
      Analyser(Omega, ++step, (time += dt));
      Timing("Analyser...");        //OVERLOAD CALL: Timing: Solve_Stokes.c(?), nekcomp.h(?), Solve_Stokes.c(?), nekcomp.h(?), drive.c(?), drive.c(?)
      
    }
    break;
  }
  printf("User time of solver (seconds):  %lf \n",
	 (clock()-begin_clock)/(double)CLOCKS_PER_SEC);

  PostProcess(Omega, step, time);
  
  exit(0);
}
  


void P_solve(Element_List *U, Element_List *Uf, 
	   Bndry *Ubc, Bsystem *Ubsys,SolveType Stype)
{
  if(!Nek_Fullcg){
    SetBCs (U,Ubc,Ubsys);        //OVERLOAD CALL: SetBCs: Solve.c(?), Solve.c(?)
    Solve  (U,Uf,Ubc,Ubsys,Stype);
  }
  else{
    Solve_CG(U,Uf,Ubc,Ubsys,Stype);
  }

}

void V_solve(Element_List *U, Element_List *Uf, 
	   Bndry *Ubc, Bsystem *Ubsys, double *us, SolveType Stype, int step){
  if(step && step < Je){
    int k;
    for(k=0;k<U->nel;++k)
      Ubsys->lambda[k].d = Re*getgamma()/dt;        //OVERLOAD CALL: getgamma: Integrate.c(?), Integrate.c(?)
    
    if(Ubc&&Ubc->DirRHS){
      free(Ubc->DirRHS);
      Ubc->DirRHS = (double*) 0;
    }

    if(U->fhead->type == 'u' || option("REFLECT")){
      Bsystem_mem_free(Ubsys, U);
      GenMat(U,Ubc,Ubsys,Ubsys->lambda,Helm);
    }
    DirBCs(U,Ubc,Ubsys,Helm);
  }

  if(option("tvarying")){
    Bndry *Bc;
    for(Bc=Ubc;Bc;Bc=Bc->next)
      Bc->elmt->update_bndry(Bc);        //OVERLOAD CALL: update_bndry: Boundary.c(Tri), Boundary.c(Quad), Boundary.c(Tet), Boundary.c(Pyr), Boundary.c(Prism), Boundary.c(Hex), Boundary.c(Element)
  }

  dcopy(U->hjtot, us, 1, U->base_hj, 1);
  if(!Nek_Fullcg){
    SetBCs (U,Ubc,Ubsys);        //OVERLOAD CALL: SetBCs: Solve.c(?), Solve.c(?)
    Solve  (U,Uf,Ubc,Ubsys,Stype);
  }
  else{
    Solve_CG(U,Uf,Ubc,Ubsys,Stype);
  }

}

void S_solve(Domain *Omega){     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)
  Element_List *V[DIM+1], *Vf[DIM+1];
  Bndry   *Vbc[DIM];
  Bsystem *Vbsys[DIM+1];
  
  V[0] = Omega->U;
  V[1] = Omega->V;
  V[2] = Omega->P;

  Vf[0] = Omega->Uf;
  Vf[1] = Omega->Vf;
  Vf[2] = Omega->Pf;

  Vbsys[0] = Omega->Usys;
  Vbsys[1] = Omega->Vsys;
  Vbsys[2] = Omega->Pressure_sys;

  Vbc[0] = Omega->Ubc;
  Vbc[1] = Omega->Vbc;

  dcopy(Omega->U->hjtot, Omega->us[0], 1, Omega->U->base_hj, 1);
  dcopy(Omega->V->hjtot, Omega->vs[0], 1, Omega->V->base_hj, 1);

  SetBCs (Omega->U,Omega->Ubc,Omega->Usys);        //OVERLOAD CALL: SetBCs: Solve.c(?), Solve.c(?)
  SetBCs (Omega->V,Omega->Vbc,Omega->Vsys);        //OVERLOAD CALL: SetBCs: Solve.c(?), Solve.c(?)

  Solve_Stokes  (V,Vf,Vbc,Vbsys);
}

void MakeF(Domain *omega, ACTION act, SLVTYPE slvtype){     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?); ACTION: nekcomp.h(?), nektar.h(?)

  Element_List  *U    =  omega->U,  *V    =  omega->V, 
                *Uf   =  omega->Uf, *Vf   =  omega->Vf,
                *P    =  omega->P,  *Pf   =  omega->Pf;

  switch (act) {

  case Prep: /* put fields in physical space for waveprop */
    U->Trans(U,J_to_Q);        //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
    V->Trans(V,J_to_Q);        //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
    break;
    
  case StokesForce:
    /* at present theta factor is dealt with in the H_Integrate.C */
    /*dsmul(U->htot,1/dt,U->base_h,1,Uf->base_h,1);
    dsmul(V->htot,1/dt,V->base_h,1,Vf->base_h,1); */
    U->Iprod(Uf);        //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
    V->Iprod(Vf);        //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
    dsmul(Uf->hjtot,1/dt,Uf->base_hj,1,Uf->base_hj,1);
    dsmul(Vf->hjtot,1/dt,Vf->base_hj,1,Vf->base_hj,1);
    break;

  case Rotational:
    VxOmega (omega);
    goto AddForcing; 

  case Convective:
    VdgradV (omega); 
    goto AddForcing; 

  case Oseen: case Stokes: case StokesS:
    StokesBC (omega);
    

  AddForcing: {
    double fx = dparam("FFX");
    double fy = dparam("FFY");

#ifdef WOMERSLEY
      fx = -dparam("WOMA")*sin(dparam("WOMW")*dparam("TIME"));
      dsadd(U->htot,fx, Uf->base_h, 1, Uf->base_h, 1);
#else
    if(fx)
      dsadd(U->htot, fx, Uf->base_h, 1, Uf->base_h, 1);
    if(fy)
      dsadd(U->htot, fy, Vf->base_h, 1, Vf->base_h, 1);

    if(omega->ForceFuncs){
      dvadd(U->htot, omega->ForceFuncs[0], 1, Uf->base_h, 1, Uf->base_h, 1);
      dvadd(V->htot, omega->ForceFuncs[1], 1, Vf->base_h, 1, Vf->base_h, 1);
    }
#endif

    if(slvtype == Splitting) SetPBCs(omega);
    break;
  }

  case Pressure: { 
    double    dtinv = 1./dt;
    double *nul= (double*)0;
    static double *ux = (double*)0;

    if(!ux)
      ux = dvector(0,2*U->htot-1);

    double *vy = ux+U->htot;

    U->Grad_d(  ux, nul, nul, 'x');        //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
    V->Grad_d(nul,   vy, nul, 'y');        //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
    dvadd(U->htot, ux, 1, vy, 1, Pf->base_h, 1);
    Pf->Set_state('p');

    Pf->Iprod(Pf);        //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
    dscal(Pf->hjtot, dtinv, Pf->base_hj, 1);

    break;
  }
    
  case Viscous: {
    double dtinv = 1/dt;

    P->Trans(P, J_to_Q);         //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
    P->Set_state('p');
    P->Grad_d (Uf->base_h,Vf->base_h,(double*)0,'a');        //OVERLOAD CALL: Grad_d: Element_List.c(Element_List), Gradient.c(Tri), Gradient.c(Quad), Gradient.c(Tet), Gradient.c(Pyr), Gradient.c(Prism), Gradient.c(Hex)
    P->Set_state('t');

    daxpy(U->htot, -dtinv, U->base_h, 1, Uf->base_h, 1);
    Uf->Set_state('p');    
    Uf->Iprod(Uf);            //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
    dscal(U->hjtot, Re, Uf->base_hj, 1);
    dcopy(U->hjtot, *omega->us, 1, U->base_hj, 1);

    daxpy(V->htot, -dtinv, V->base_h, 1, Vf->base_h, 1);
    Vf->Set_state('p');    
    Vf->Iprod(Vf);        //OVERLOAD CALL: Iprod: InnerProd.c(Element), Element_List.c(Element_List), InnerProd.c(Tri), InnerProd.c(Quad), InnerProd.c(Tet), InnerProd.c(Pyr), InnerProd.c(Prism), InnerProd.c(Hex)
    dscal(V->hjtot, Re, Vf->base_hj, 1);
    dcopy(U->hjtot, *omega->vs, 1, V->base_hj, 1);
    
    break;
  }
  case Post: {
    int Nm = U->hjtot;

    if((slvtype == StokesSlv)){
      double theta = dparam("THETA");

      dsvtvp(Nm,-theta,*omega->us,1,U->base_hj,1,U->base_hj,1);
      dscal (Nm,1.0/(1.0-theta),U->base_hj,1);
      dsvtvp(Nm,-theta,*omega->vs,1,V->base_hj,1,V->base_hj,1);
      dscal (Nm,1.0/(1.0-theta),V->base_hj,1);
    }

    dcopy(Nm, U->base_hj, 1, *omega->us, 1);
    dcopy(Nm, V->base_hj, 1, *omega->vs, 1);
    break;
  }
  default:
    error_msg(MakeF--unknown type of action);        //OVERLOAD CALL: error_msg: gen_utils.h(?), hotel.h(?)
    break;
  }

  return;
}

/* Do initial time step assuming for case where startup field is in physical *
 * space or copy V field to Vs if it is a restart                           */

static void StartUp(Domain *omega){     //OVERLOAD CALL: Domain: nekcomp.h(?), nektar.h(?)

  if(omega->U->fhead->state == 'p'){
    omega->U->Trans(omega->U,Q_to_J);        //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
    SetBCs (omega->U,omega->Ubc,omega->Usys);        //OVERLOAD CALL: SetBCs: Solve.c(?), Solve.c(?)
    omega->V->Trans(omega->V,Q_to_J);        //OVERLOAD CALL: Trans: Element_List.c(Element_List), Transforms.c(Element), Transforms.c(Tri), Transforms.c(Quad), Transforms.c(Tet), Transforms.c(Pyr), Transforms.c(Prism), Transforms.c(Hex)
    SetBCs (omega->U,omega->Ubc,omega->Usys);        //OVERLOAD CALL: SetBCs: Solve.c(?), Solve.c(?)
  }

  dcopy(omega->U->hjtot, omega->U->base_hj, 1, *omega->us, 1);
  dcopy(omega->V->hjtot, omega->V->base_hj, 1, *omega->vs, 1);
  omega->U->Set_state('t');
  omega->V->Set_state('t');
  omega->P->Set_state('t');
}

void switch_hj(Element_List *EL, double **hj){
  double *d = EL->base_hj;
  EL->Mem_shift(EL->base_h, *hj);        //OVERLOAD CALL: Mem_shift: Element_List.c(Element_List), Memory.c(Tri), Memory.c(Quad), Memory.c(Tet), Memory.c(Pyr), Memory.c(Prism), Memory.c(Hex), Memory.c(Element)
  *hj = d;
}







Back to Source File Index


C++ to HTML Conversion by ctoohtml