Subversion Repositories channelflow

[/] [trunk/] [programs/] [fieldprops.cpp] - Rev 447

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"
#include "channelflow/turbstats.h"
#include "channelflow/periodicfunc.h"

using namespace std;
using namespace channelflow;

Real energy(FlowField& u, int i, bool normalize=true);

// returns S,A, for symmetric,antisymmetric to eps accuracy
// returns s,a, for symmetric,antisymmetric to sqrt(eps) accuracy
char symmetric(const FieldSymmetry& s, const FlowField& u, Real& serr, Real& aerr, Real eps=1e-7);

Real project(const FlowField& u, const FieldSymmetry& s, int sign);

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

  string purpose("prints information about a FlowField");

  ArgList args(argc, argv, purpose);
  const bool geom     = args.getflag("-g",  "--geometry", "show geometrical properties");
  const bool mean     = args.getflag("-m",  "--mean",     "show mean properties");
  const bool norm     = args.getflag("-n",  "--norm",     "show norm properties");
  const bool spec     = args.getflag("-sp", "--spectral", "show spectral properties");
  const bool symm     = args.getflag("-sy", "--symmetry", "show symmetry properties");
  const bool dyn      = args.getflag("-d",  "--dynamic",  "show dynamical properties");
  const bool nrgy     = args.getflag("-e",  "--energy",   "show energy properties");
  const bool wall     = args.getflag("-w",  "--wall",     "show wall-unit properties");
  const bool local    = args.getflag("-l",  "--local",    "show z-localization");

  const Real Reynolds = args.getreal("-R", "--Reynolds", 400.0, "Reynolds number == (ub-ua)*(b-a)/4nu ");
  /***/ Real nu       = args.getreal("-nu", "--viscosity", 0.0, "kinematic viscosity (takes precedence over Reynolds, if nonzero");
  const string Ustr   = args.getstr("-U", "--Ubase", "",  "zero|laminar|parabolic|filename base flow U(y)"); 

  const string meanstr= args.getstr("-mc", "--meanconstraint", "gradp", "gradp|bulkv : hold one fixed");
  const Real dPdx     = args.getreal("-P", "--dPdx",   0.0, "value for fixed pressure gradient");
  //const Real Ubulk    = args.getreal("-U", "--Ubulk",  0.0, "value for fixed bulk velocity");
  const Real ua       = args.getreal("-ua", "--ulower",   -1.0, "lower wall speed");
  const Real ub       = args.getreal("-ub", "--uupper",    1.0, "upper wall speed");

  const Real eps      = args.getreal("-eps", "--eps", 1e-7, "spectral noise floor");

  const Real dt       = args.getreal("-dt", "--dt", 0.03125, "integration time for du/dt estimate");

  const Real T        = args.getreal("-T", "--T", dt, "integration time for du/dt estimate");

  const string nlmeth = args.getstr("-nl", "--nonlinearity", "rot", "method for  u grad u [rot or skw]");

  const int  digits   = args.getint("-dg", "--digits", 6, "digits of output for reals");
  const string uname  = args.getstr(1, "<flowfield>", "input field");
  args.check();

  bool all = (!geom && !mean && !norm && !spec && !symm && !dyn && !nrgy && !wall && !local) ? true : false;

  FlowField u(uname);
  u.makeSpectral();
  const Real unorm = L2Norm(u);
  cout << setprecision(digits);

  ChebyCoeff U(u.Ny(), u.a(), u.b());
  if (Ustr.length() == 0 || Ustr == "linear" || Ustr == "planecouette" || Ustr == "pcf") 
    U[1] = 1;
  else if (Ustr == "parabolic" || Ustr == "channel") {
    U[0] = 0.5;
    U[2] = -0.5;
  }
  else if (Ustr != "zero") {
    U = ChebyCoeff(Ustr);
    U.makeSpectral();
  }

  cout << "Ubase mean == " << U.mean() << endl;
  //U.save("U");

  if (nu == 0) 
    nu = abs(0.25*(ub-ua)*(u.b()-u.a())/Reynolds);
  
  if (all || geom) {
    cout << "--------------Geometry------------------" << endl;
    cout << "Nx == " << u.Nx() << endl;
    cout << "Ny == " << u.Ny() << endl;
    cout << "Nz == " << u.Nz() << endl;
    cout << "Nd == " << u.Nd() << endl;
    cout << "Lx == " << u.Lx() << endl;
    cout << "Ly == " << u.Ly() << endl;
    cout << "Lz == " << u.Lz() << endl;
    cout << " a == " << u.a() << endl;
    cout << " b == " << u.b() << endl;
    cout << "lx == " << u.Lx()/(2*pi) << endl;
    cout << "lz == " << u.Lz()/(2*pi) << endl;
    cout << "alpha == " << 2*pi/u.Lx() << endl;
    cout << "gamma == " << 2*pi/u.Lz() << endl;
    cout << "Lx/Lz == " << u.Lx()/u.Lz() << endl;
    cout << "Lx*Lz == " << u.Lx()*u.Lz() << endl;
    cout << "diag  == " << pythag(u.Lx(),u.Lz()) << " == sqrt(Lx^2 + Lz^2)" << endl;
    cout << "vol   == " << u.Lx()*u.Ly()*u.Lz()  << " == Lx*Ly*Lz" << endl;
    cout << "nu    == " << u.nu()  << endl;
    cout << endl;
  }
  if (all || mean) {
    u += U;
    for (int i=0; i<u.Nd(); ++i) {
      cout << "mean u[" << i << "] == " <<  Re(u.profile(0,0,i)).mean() << endl;
    }
    u -= U;
  }
  if (all || symm) {
    cout << "--------------Symmetry-------------------\n";
    FieldSymmetry s1( 1, 1,-1, 0.5, 0.0); // Wally-oriented symmetries, GHC09 S group
    FieldSymmetry s2(-1,-1, 1, 0.5, 0.5);
    FieldSymmetry s3(-1,-1,-1, 0.0, 0.5);

    FieldSymmetry sxtxz(-1,-1, 1, 0.5, 0.5); // GHC09 Rxz group, conjugate to S.
    FieldSymmetry sztxz( 1, 1,-1, 0.5, 0.5);
    FieldSymmetry sxz  (-1,-1,-1, 0.0, 0.0); 
    
    FieldSymmetry sx(-1,-1, 1);
    FieldSymmetry sz( 1, 1,-1);
    FieldSymmetry tx( 1, 1, 1, 0.5, 0.0);
    FieldSymmetry tz( 1, 1, 1, 0.0, 0.5);
    FieldSymmetry txz(1, 1, 1, 0.5, 0.5);

    Real unorm2 = L2Norm2(u);
    cout.setf(ios::left);
    cout << "fractions of energy in symmetric/antisymmetric subspaces:" << endl;
    cout << "   s1  symm: " << setw(12) << project(u, s1, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, s1,-1)/unorm2 << endl;
    cout << "   s2  symm: " << setw(12) << project(u, s2, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, s2,-1)/unorm2 << endl;
    cout << "   s3  symm: " << setw(12) << project(u, s3, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, s3,-1)/unorm2 << endl << endl;

    cout << "sxtxz  symm: " << setw(12) << project(u, sxtxz, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, sxtxz,-1)/unorm2 << endl;
    cout << "sztxz  symm: " << setw(12) << project(u, sztxz, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, sztxz,-1)/unorm2 << endl;
    cout << "sxz    symm: " << setw(12) << project(u, sxz, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, sxz,-1)/unorm2 << endl << endl;

    cout << " sx    symm: " << setw(12) << project(u, sx, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, sx,-1)/unorm2 << endl;
    cout << " sz    symm: " << setw(12) << project(u, sz, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, sz,-1)/unorm2 << endl;
    cout << "   tx  symm: " << setw(12) << project(u, tx, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, tx,-1)/unorm2 << endl;
    cout << "   tz  symm: " << setw(12) << project(u, tz, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, tz,-1)/unorm2 << endl;
    cout << "   txz symm: " << setw(12) << project(u, txz, 1)/unorm2;
    cout << "       anti: " << setw(12) << project(u, txz,-1)/unorm2 << endl;
    cout << endl;
    cout.unsetf(ios::left);
  }

  if (all || norm) {
    cout << "--------------Norms---------------------" << endl;

    cout << " L2Norm(u)    == " <<  unorm << endl;
    for (int i=0; i<u.Nd(); ++i)
      cout << "  L2Norm(u[" << i << "]) == " << sqrt(energy(u,i)) << endl;
    for (int i=0; i<u.Nd(); ++i)
      cout << "LinfNorm(u[" << i << "]) == " << LinfNorm(u[i]) << endl;
    cout << "L2Norm3d(u)   == " <<  L2Norm3d(u) << endl;
    cout << " divNorm(u)   == " << divNorm(u) << endl;
    cout << "  bcNorm(u)   == " <<  bcNorm(u) << endl;
    cout << endl;
  }

  if (all || spec) {
    cout << "------------Spectrum--------------------" << endl;
    int kxmax=0;
    int kzmax=0;
    int kymax=0;
    cout << "u.kxmin() == " << u.kxmin() << endl;
    cout << "u.kxmax() == " << u.kxmax() << endl;
    cout << "u.kzmin() == " << u.kzmin() << endl;
    cout << "u.kzmax() == " << u.kzmax() << endl;
    cout << "u.kxminDealiased() == " << u.kxminDealiased() << endl;
    cout << "u.kxmaxDealiased() == " << u.kxmaxDealiased() << endl;
    cout << "u.kzminDealiased() == " << u.kzminDealiased() << endl;
    cout << "u.kzmaxDealiased() == " << u.kzmaxDealiased() << endl;
    
    Real profmax = 0.0;
    for (int mx=0; mx<u.Mx(); ++mx) {
      const int kx = u.kx(mx);
      for (int mz=0; mz<u.Mz(); ++mz) {
        const int kz = u.kz(mz);
        BasisFunc prof = u.profile(mx,mz);
        if (L2Norm(prof) > profmax) {
          kxmax = Greater(kxmax, abs(kx));
          kzmax = Greater(kzmax, abs(kz));
        }
        for (int ky=0; ky<u.Ny(); ++ky) {
          Real sum = 0.0;
          for (int i=0; i<u.Nd(); ++i)
            sum += abs(prof[i][ky]);
          if (sum > eps)
            kymax = Greater(kymax, ky);
        }
      }
    }
    cout << "u.padded() == " << (u.padded() ? "true" : "false") << endl;
    cout << "Energy over " << eps << " is confined to : \n";
    cout << " |kx| <= " << kxmax << endl;
    cout << " |ky| <= " << kymax << endl;
    cout << " |kz| <= " << kzmax << endl;

    cout << "Minimum   aliased grid : "
         << 2*(kxmax+1) << " x "
         << kymax+1 << " x "
         << 2*(kzmax+1) << endl;

    cout << "Minimum dealiased grid : "
         << 3*(kxmax+1) << " x "
         << kymax+1 << " x "
         << 3*(kzmax+1) << endl;

    int kxmin = u.padded() ? u.kxminDealiased() : u.kxmin();
    kxmax = u.padded() ? u.kxmaxDealiased() : u.kxmax();
    kzmax = u.padded() ? u.kzmaxDealiased() : u.kzmax();

    int mzmax = u.mz(kzmax);
    int mymax = u.My()-1;


    // x truncation: Search over kz for largest spectral coeff of (kxmin,kz) or (kzxmax,kz).
    // Need to check both kxmin and kxmax, the neg/pos extrema of x Fourier modes
    // since the Fourier representation of real-valued fields includes both of these

    int kxtrunc = 0;
    int kztrunc = 0;
    Real maxtrunc = 0.0;

    for (int kx=kxmin; kx<=kxmax; kx += (kxmax-kxmin)) {
      int mx = u.mx(kx);

      // Search over kz values for a large mode
      for (int mz=0; mz<=mzmax; ++mz) {
        BasisFunc prof = u.profile(mx,mz);
        Real trunc = L2Norm(prof);
        if (trunc > maxtrunc) {
          kxtrunc = kx;
          kztrunc = u.kz(mz);
          maxtrunc = trunc;
        }
      }
    }
    cout << "Max x truncation == " << setw(12) << maxtrunc << " at kx,kz == " << kxtrunc << ',' << kztrunc << endl;


    // y-truncation: search over all (kx,kz) for the largest coefficient of the last Chebyshev mode
    maxtrunc = 0.0;
    kxtrunc = 0;
    kztrunc = 0;

    for (int i=0; i<u.Nd(); ++i) {
      for (int kx=kxmin; kx<=kxmax; ++kx) {
        int mx = u.mx(kx);
        for (int kz=0; kz<=kzmax; ++kz) {
          int mz = u.mz(kz);
          Real trunc = abs(u.cmplx(mx,mymax,mz,i));
          if (trunc > maxtrunc) {
            kxtrunc = kx;
            kztrunc = kz;
            maxtrunc = trunc;
          }
        }
      }
    }     
    cout << "Max y truncation == " << setw(12) << maxtrunc << " at kx,kz == " << kxtrunc << ',' << kztrunc << endl;

    // z-truncation: Search over kx for largest spectral coeff of (kx,kzmax)
    // In this case (as opposed to the x-truncation search above), negative kz values are
    // not stored but rather implicitly defined by complex symmetry for real-valued field.
    // So there is no search over (kz,kzmin)

    maxtrunc = 0.0;
    kxtrunc = 0;
    int mz = u.mz(kzmax);
    cout << "mz == " << mz << endl;

    for (int kx=kxmin; kx<=kxmax; ++kx) {
      int mx = u.mx(kx);
      //int mz = u.mz(kz);
      BasisFunc prof = u.profile(mx,mz);
      Real trunc = L2Norm(prof);
      if (trunc > maxtrunc) {
        kxtrunc = kx;
        maxtrunc = trunc;
      }
    }
    cout << "Max z truncation == " << setw(12) << maxtrunc << " at kx,kz == " << kxtrunc << ',' << kzmax << endl;

  }

  if (all || nrgy) {
    FlowField uU(u);
    uU += U;

    cout << "dissip (u)     == " << dissipation(u) << endl;
    cout << "wallshear(u)   == " <<  0.5*(abs(wallshearUpper(u))+abs(wallshearLower(u))) << endl;
    cout << "energy(u)      == " << 0.5*L2Norm2(u) << endl;
    cout << "dissip (u+U)   == " << dissipation(uU) << endl;
    cout << "wallshear(u+U) == " << 0.5*(abs(wallshearUpper(uU))+abs(wallshearLower(uU))) << endl;
    cout << "energy (u+U)   == " << 0.5*L2Norm2(uU) << endl;
    cout << "energy3d(u)    == " << 0.5*L2Norm2_3d(uU) << endl;
  }

  if (all || dyn) {
    cout << "------------Dynamics--------------------" << endl;

    // Construct Navier-Stoke Integrator
    DNSFlags flags;
    flags.timestepping = SBDF3;
    flags.dealiasing   = DealiasXZ;
    flags.nonlinearity = Rotational;
    flags.constraint   = s2constraint(meanstr);
    flags.baseflow     = LaminarBase;
    flags.ulowerwall   = ua;
    flags.uupperwall   = ub;
    flags.nu           = nu;
    flags.dPdx         = dPdx;

    if (nlmeth == "skw")
      flags.nonlinearity = SkewSymmetric;
    else if (nlmeth == "div")
      flags.nonlinearity = Divergence;
    else if (nlmeth == "cnv")
      flags.nonlinearity = Convection;
    else if (nlmeth == "alt")
      flags.nonlinearity = Alternating;
    else
      flags.nonlinearity = Rotational;


    // Define flow parameters

    const int Nx = u.Nx();
    const int Ny = u.Ny();
    const int Nz = u.Nz();
    const Real Lx = u.Lx();
    const Real  a = u.a();
    const Real  b = u.b();
    const Real Lz = u.Lz();
    //const Real dt = (T < 0.01) ? T : 0.01;

    const int N = iround(T/dt);

    // Construct base flow for plane Couette: U(y) = y
    Vector y = chebypoints(Ny,a,b);

    FlowField q(Nx,Ny,Nz,1,Lx,Lz,a,b);
    DNS dns(u, flags);

    cout << "computing du/dt: " << flush;
    FlowField dudt(u);
    dns.advance(dudt, q, N);
    cout << endl;
    dudt -= u;
    dudt *= 1.0/(N*dt);
    cout << "L2Norm(u)     == " << L2Norm(u) << endl;
    cout << "L2Norm(du/dt) == " << L2Norm(dudt) << endl;

    FlowField dudx;
    FlowField dudz;
    xdiff(u,dudx);
    zdiff(u,dudz);

    Real dudxnorm = L2Norm(dudx);
    Real dudznorm = L2Norm(dudz);

    cout << "L2Norm(du/dx) == " << dudxnorm << endl;
    cout << "L2Norm(du/dz) == " << dudznorm << endl;

    Real ueps = 1e-12;
    dudx *= (dudxnorm < ueps) ? 0.0 : 1.0/dudxnorm;
    dudz *= (dudznorm < ueps) ? 0.0 : 1.0/dudznorm;

    Real ax = L2IP(dudt,dudx);
    Real az = L2IP(dudt,dudz);

    FlowField Pu(u);
    Pu.setToZero();

    FlowField utmp;
    utmp = dudx;
    utmp *= ax;
    Pu += utmp;

    utmp = dudz;
    utmp *= az;
    Pu += utmp;

    Real Punorm   = L2Norm(Pu);
    Real waviness = Punorm/L2Norm(dudt);
    Real theta = acos(L2IP(Pu,dudt)/(L2Norm(dudt)*Punorm));

    cout << endl;
    cout << " waviness == " << waviness << endl;
    cout << "    theta == " << theta << " == " << theta*180/pi << " degrees\n";
    cout << endl;
    cout << "(waviness == fraction of dudt within dudx,dudz tangent plane)\n";
    cout << "(   theta == angle between dudt and dudx,dudz tangent plane)\n";
    cout << endl;
  }

  if (all | wall) {
    TurbStats stats(U, nu);
    FlowField tmp(u.Nx(), u.Ny(), u.Nz(), 9, u.Lx(), u.Lz(), u.a(), u.b());
    stats.addData(u, tmp);
    
    cout << "ustar == " << stats.ustar() << endl;
    cout << "hplus == " << stats.hplus() << endl;
    cout << "saving Reynolds stresses etc to turbstats.asc" << endl;
    stats.msave("turbstats.asc");
  }
  if (all | local) {
    FlowField e = energy(u);
    e.makeSpectral();

    PeriodicFunc exyavg(e.Nz(), e.Lz(), Spectral);
    for (int mz=0; mz<e.Mz(); ++mz)
      exyavg.cmplx(mz) = e.cmplx(0,0,mz,0);
    exyavg.makePhysical();
    exyavg.save("exyavg");
    Real emin = exyavg(0);
    Real emax = emin;

    for (int nz=1; nz<e.Nz(); ++nz) {
      emin = lesser(emin, exyavg(nz));
      emax = Greater(emax, exyavg(nz));
    }
    cout << "sqrt(min(exyavg(z))/max(exyavg(z))) == " << sqrt(emin/emax) << endl; 
  }
}

Real energy(FlowField& u, int i, bool normalize) {
  assert(u.ystate() == Spectral);
  assert(u.xzstate()==Spectral);
  Real sum = 0.0;
  ComplexChebyCoeff prof(u.Ny(), u.a(), u.b(), Spectral);
  for (int mx=0; mx<u.Mx(); ++mx) {
    Real cz = 1.0; // cz = 2 for kz>0 to take account of kz<0 ghost modes
    for (int mz=0; mz<u.Mz(); ++mz) {
      for (int ny=0; ny<u.Ny(); ++ny)
        prof.set(ny, u.cmplx(mx,ny,mz,i));
      sum += cz*L2Norm2(prof, normalize);
      cz = 2.0;
    }
  }
  if (!normalize)
    sum *= u.Lx()*u.Lz();
  return sum;
}

char symmetric(const FieldSymmetry& s, const FlowField& u, Real& serr,Real& aerr, Real eps) {
  Real sqrteps= sqrt(eps);

  FlowField su = s(u);
  FlowField us, ua;
  ((us = u) += su) *= 0.5;
  ((ua = u) -= su) *= 0.5;
  serr = L2Dist(us, u);
  aerr = L2Dist(ua, u);

  // Check for symmetry
  if (serr < eps)
    return 'S';
  else if (serr < sqrteps)
    return 's';
  else if (aerr < eps)
    return 'A';
  else if (aerr < sqrteps)
    return 'a';
  else
    return ' ';
};

Real project(const FlowField& u, const FieldSymmetry& s, int sign) {
  FlowField Pu(u);
  if (sign<0)
    Pu -= s(u);
  else
    Pu += s(u);
  //Pu *= 0.5;
  return 0.25*L2Norm2(Pu);
}

Compare with Previous | Blame