CFD Utilities Files

Traj_fit Description:

   This specialized curve-fitting scheme is implemented as a particular
instance of this subroutine:

   subroutine fun (n, x, f, ncall, fail)

   Subroutine FUN contains all of the application-specifics required to
perform unconstrained optimization via QNMDIF2.  The rest of the QNMDDRIVER2
framework should be aplication-independent.

   This version of FUN performs one nonlinear least squares fit of the form

                f(rho(t),V(t)) = C rho(t)**m V(t)**n

to a data table of flow solver and radiation solver results at a moderate
number of "CFD points" (flight conditions) along a trajectory for some body
point on an atmospheric entry vehicle.

   Here, linear coefficient C and exponents m and n are the three variables
being optimized and rho and V vary with time along the trajectory, for which
a full time history is read as data.  The function f(rho,V) typically forms
a pulse suited to this formulation.  For instance, f may be the pressure or
the heat flux at the nose/stagnation point of the vehicle, where this kind
of correlation is best understood.  In practice, the method will be applied
to conditions at other body points, including on the aft body where results
may be questionable at best.  Fitting of shear stress and heat transfer
coefficient data are also likely to be required as part of preparing the
surface environment data for a material response solver such as FIAT.

   Unlike its predecessor, PQfit, and the alternative Gnuplot approach, this
method provides an evaluated time history of the fitted curve, thus enabling
plots of good quality via (say) Tecplot.  (Gnuplot can plot its fitted
results - crudely - but does not seem to have a way to save them as output.)

   QNMDIF2 is chosen rather than a constrained optimization package because
it is not proprietary, and because automatic scaling of the data being fit
is expected to enable the same starting guesses for C, m, n (1., 1., 2. from
past experience) to suffice for any likely application, although provision
is made for supplying other starting guesses.  The lack of variable bounds
might be an issue, because C, m, n must be positive.  Experience, however,
suggests that this is not a likely problem.  Extremely sharp pulses, though,
may require trial and error with the starting guesses to achieve the best
fit.  Converging too quickly is another sign that other starts should be
tried.  Essentially the same solution from more than one start should be
reassuring.  Some datasets may be so un-pulse-like that many iterations may
be seen without proper convergence, and the same would be true of any other
fitting method.

   The dataset formats shown below arose from related efforts to automate
(or partially automate) the numerous steps leading up to this curve-fitting
step.  Using the # comment character in header lines is recommended for
Tecplot compatibility.

Namelist Control File Example (entered on standard input):

   $fit
   traj_file = '/home/dasaund1/VASE/EFPA-14/freestream-14.dat'
   CFD_file  = '~dasaund1/VASE/EFPA-14/Gamma-14-stag-point-data-to-fit.dat'
   CFD_col_fit = 4
   $end

   This example fits pw (column 4) using default settings.
   The full namelist controls are listed below.

Full Trajectory Data Format (as from program Traj, any # top header lines):

   #t,s V,km/s  pinf,Pa rho,kg/m^3 Tinf,K pstag,Pa Mach  qconv  qrad,W/cm^2
   0.00 11.938 2.80e-07 1.97e-12  258.68  2.80e-04 46.6    0.1    0.0
   0.10 11.938 2.84e-07 2.00e-12  258.68  2.85e-04 46.6    0.1    0.0
   0.20 11.938 2.88e-07 2.03e-12  258.67  2.89e-04 46.6    0.1    0.0
    :     :     :        :           :     :         :      :      :
   [Any embedded headers should be removed.  E.g., FILTER_LINES utility.]

CFD Data Table Format (as from DPLR/BLAYER/NEQAIR) (squeezed/split to fit):

   #t,s  rho,kg/m^3 V,km/s pw,Pa      qw,W/m^2   Tw,K       CH,kg/m^2.s
   47.00 1.18e-04   11.944 8.4162E+03 4.2334E+06 3.0626E+03 6.8289E-02 ...
   50.00 4.94e-04   11.740 3.3895E+04 1.1987E+07 3.9728E+03 2.0642E-01 ...
   52.10 1.25e-03   11.309 8.0023E+04 2.1407E+07 4.5926E+03 4.1149E-01 ...
     :    :           :     :          :          :          :

   ... He,J/kg    Hw,J/kg    tauw,Pa    qrad,W/cm^2  Hfree,J/kg
   ... 6.2704E+07 7.2142E+05 4.2301E+02 1.102117E+02 6.2714703872E+07
   ... 6.0445E+07 2.2254E+06 1.2268E+03 9.207419E+01 6.0296955560E+07
   ... 5.5515E+07 3.3082E+06 2.2875E+03 1.257758E+02 5.5332140482E+07
        :          :          :          :            :
   At present, the default is to weight all data points equally, but some
   physics-based weighting scheme may prove desirable.

Output Results Format (1):  Optimized coefficients C, m, n

   2.67037E-03              These are written to file fn.coefs, where
   5.43278E-01              n = CFD_col_fit.  Ideally, they won't be
   4.12809E+00              needed, but they may be informative.

Output Results Format (2):  Fitted curve evaluated along full trajectory

   #t,s  f(rho(t),V(t))     These are written to file fn.eval.dat, where
   0.00  1.23456E-07        n = CFD_col_fit.
   0.10  1.45678E-07
   0.20  1.56789E-07
    :     :

Output Log File 'qnmdif2.out':

   The QNMDIF2 iteration log appears in qnmdif2.out (hard-coded name).  Using
   grep "AT ITER" qnmdif2.out gives a quick clue as to the effectiveness of
   the fit.

Further Notes:

   o  QNMDIF2 (Quasi-Newton Method, finite-difference first derivatives) is
      a general purpose unconstrained optimization package closely related to
      the original QNMDIF which dates from the '70s.  Quasi-Newton refers to
      its use of not the Hessian matrix of second derivatives but rather a
      Hessian-like matrix that starts as a diagonal matrix here.  QNMDIF's
      use of Cholesky factors of this evolving matrix (possibly modified to
      ensure positive definiteness and hence descent at each quasi-Newton
      step) distinguishes it from other such methods.  The line search is
      derivative-free, while gradients are estimated with central differences
      for which the differencing intervals are computed at run time via a
      preliminary CENDIF2 scheme.  (These "h"s can also be input > 0. if
      desired, but here, that is not recommended for this application.)

Full Namelist Control Description:

   Defaults for most of the QNMDIF2-related controls are set in the reusable
   subroutine qnm_solution.f90 (see /optlib).  They may be added to the name-
   list $fit used here, but should not need to be.

   It is recommended that native solver units appear in the data tables.  The
   default column scales suit the indicated unit assumptions for rho, V, etc.

                                                                   Default
   CFD_col_rho          CFD dataset column for density             2; kg/m^3
   CFD_col_V             "     "     "     "   V                   3; km/s
   CFD_col_fit           "     "     "     "   f(rho,V)            4; pw,Pa
   CFD_filename          "     "    file name; no default
   CFD_col_scale(:)      "     "    column scale factors           See below
   CFD_row_weight(:)     "     "    row weights in sum of squares  1.
   traj_col_t           Traj. dataset column for time              1; s
   traj_col_rho         Traj. dataset column for density           4; kg/m^3
   traj_col_V            "     "     "     "     V                 3; km/s
   traj_filename         "     "    file name; no default
   niter                Iteration limit for QNMDIF2                1000
   h(:)                 Finite difference intervals for C, m, n;   -1. (auto)
   x0(1:3)              Override 1.,1.,2. C,m,n starts if > 0.     -1. (auto)
   
History:

   08/23/94  D.A.Saunders  Initial implementation of reusable framework.
   02/18/95  DAS  NCALL argument added; EPSMCH & STEPMX now initialized.
   08/30/11  DAS  Use the intrinsic EPSILON now instead of the old EPSLN.
   07/12/13  DAS  *.f --> *.f90 form and common /QNMCTL/ replaced by
                  qnm_module.
   09/15/14- DAS  Trajectory data fitting version of FUN.  This prompted the
   09/19/14       table_io.f90 module first used here.
   09/24/14  DAS  Initial applications compare well with Gnuplot's Levenberg-
                  Marquardt method applied to the same problems.  The missing
                  full namelist description has been added now (above).

Author:  David Saunders, ERC, Inc./NASA Ames Research Center, CA.