Subversion Repositories channelflow

[/] [trunk/] [programs/] [diffop.cpp] - Rev 414

Compare with Previous | Blame | View Log

#include <fstream>
#include <iomanip>
#include <string>
#include <iostream>
#include "channelflow/flowfield.h"
#include "channelflow/diffops.h"
#include "channelflow/symmetry.h"
#include "channelflow/utilfuncs.h"

using namespace std;
using namespace channelflow;

int main(int argc, char* argv[]) {

  string purpose("apply a differential operation to a given field");

  ArgList args(argc, argv, purpose);

  //const bool channel  = args.getflag("-c", "--channel", "channelflow instead of plane Couette");
  //const bool couette  = args.getflag("-pcf", "--planecouette", "channelflow instead of plane Couette");

  const bool ddx = args.getflag("-ddx", "--ddx",         "apply d/dx");
  const bool ddy = args.getflag("-ddy", "--ddy",         "apply d/dy");
  const bool ddz = args.getflag("-ddz", "--ddz",         "apply d/dz");
  const bool grd = args.getflag("-grad", "--gradient",   "apply gradient");
  const bool lpl = args.getflag("-lapl", "--laplacian",  "apply laplacian");
  const bool crl = args.getflag("-curl", "--curl",       "apply curl");
  const bool dvv = args.getflag("-div",  "--divergence", "compute divergence");
  const bool nnl = args.getflag("-nonl", "--nonlinearity","compute Navier-Stokes nonlinearity");
  const string meth = args.getstr("-meth", "--nonlmethod","conv", "method for computing nonl: conv|skew|div|rot");
  const bool nrg = args.getflag("-e",    "--energy",      "compute energy operator");
  const bool qcr = args.getflag("-Q",    "--Qcriterion",  "compute Q criterion");
  const bool nrm = args.getflag("-norm", "--norm",        "compute pointwise vector norm of field");

  const string Uname  = args.getstr("-U",    "--Ubase",  "", "U(y) base flow");
  const string iname  = args.getstr(2, "<flowfield>", "input field");
  const string oname  = args.getstr(1, "<flowfield>", "output field");

  args.check();
  args.save();

  FlowField u(iname);
  u.makeSpectral();

  if (Uname == "parabolic") {
    ChebyCoeff U(u.Ny(), u.a(), u.b(), Spectral);
    U[0] =  0.5;
    U[2] = -0.5;
    u += U;
  }
  else if (Uname == "linear") {
    ChebyCoeff U(u.Ny(), u.a(), u.b(), Spectral);
    U[1] = 1.0;
    u += U;
  }
  else if (Uname.length() > 0) {
    ChebyCoeff U(Uname);
    U.makeSpectral();
    u += U;
  }

  FlowField v;

  // Note: could also use v = grad(u) etc, but grad(u,v) is more efficient
  
  if (ddx)
    xdiff(u,v);

  else if (ddy)
    ydiff(u,v);

  else if (ddz)
    zdiff(u,v);

  else if (grd)
    grad(u,v);

  else if (lpl)
    lapl(u,v);

  else if (crl)
    curl(u,v);

  else if (dvv)
    div(u,v);

  else if (nnl) {
    // Have already added Ubase to u, so don't need to here (and we're assuming Wbase==0). 
    // But navierstokeNL fuinc sig requires both. So make up zero Ubase and Wbase
    ChebyCoeff zero(u.Ny(), u.a(), u.b(), Spectral); 
    FlowField tmp(u);
    NonlinearMethod m = s2nonlmethod(meth);
    navierstokesNL(u,zero,zero,v,tmp, m);
    //convectionNL(u,v,tmp);
  }
  else if (qcr)
    Qcriterion(u,v);

  else if (nrg)
    energy(u,v);

  else if (nrm)
    norm(u,v);
  
  else 
    cferror("diffop: please choose a differential operator, rerun with -h option");
    
  v.makeSpectral();
  v.save(oname);
}

Compare with Previous | Blame