#####################
    #                   #
    #  VIPER CHANGELOG  #
    #                   #
    #####################

25 JUNE 2008:

- Fixed the symmetry boundary condition assignment - these were not 
  being detected following the alteration on 14 May 2008 to the way 
  in which Viper processed Dirichlet boundary conditions.

- Fixed a minor bug in the ROTATE implementation - when base flows have
  a non-zero w-velocity, the code now correctly includes the imaginary
  u & v components of the perturbation field when evaluation the Coriolis
  correction.


18 JUNE 2008:

- Upgraded ARNOLDI stability analysis solver to compute the stability
  of perturbation fields both with and without a w-component of velocity
  in the base flow (previously only w=0 flows were implemented).


13 JUNE 2008:

- Found a bug in the ROTATE routine - the x-direction advection RHS 
  vectors were being copied into both the x- and y-direction vectors.
  Tests of the fix to this suggest that the code is working better now.


12 JUNE 2008:

- Revising the rotating reference frame code for the perturbation field.
  Removed centrifugal acceleration component as this appears only in the
  integration of the base flow.

- Tidied up some of the outputting of floating-point numbers in scientific
  notation - tidier notation used for more values reported during 
  initialisation.


6 JUNE 2008:

- A small performance boost was achieved in parallel by using separate
  copies of the spatial derivative matrices for each processing unit.


3 JUNE 2008:

- Investigating the use of KMP_AFFINITY and extra OMP loops during 
  initialisation to get improved parallel performance on a per-processor
  basis. Thus far have sped up PARDISO computations, but other parts have
  slowed - may need to create per-thread copies of derivative matrices, etc.


30 MAY 2008:

- Early parallelization tests for the stability analysis code indicate
  a promising speedup of 4.7 times over 8 CPUs, and 6.1 times over 16 CPUs.

- Moved time evolution of perturbation fields into a subroutine - this will
  aid parallelization of the 2D time integration code.


29 MAY 2008:

- Condensed some redundant code in pressure and diffusion solves - now
  using one routine for all 2D & 3D static condensation and matrix solve
  phases of the pressure and diffusion substeps, rather than the 7
  separate routines that were being called earlier.


28 MAY 2008:

- Identified that the cylindrical SE/Fourier 3D solver was suffering from
  an incorrect calculation of theta-derivatives during the advection step.
  This was contaminating the results and leading to a spurious deviation of
  the flow near the axis, as well as instability. This has now been fixed.
  Tests against the Kovasznay flow solution verify that the solver is now
  functioning correctly.


22 MAY 2008:

- Tested the Cartesian spectral-element/Fourier 3D algorithm against an
  analytical test case - a Kovasznay flow - imposed oblique to the domain.
  The test verified that the Cartesian formulation is functioning correctly.
  Tests indicate that the cylindrical formulation is not yet correct, though
  is stable - error manifests as a slight discrepancy across the axis.

- Identified a stack access violation under Windows that was causing Viper 
  to crash during calls to "ARNOLDI". A function "second" was being called
  while no longer being defined - this has now been removed, as we are not
  utilising any timing statistics from the ARPACK package.

- Migrated compilers to Intel Fortran 10.1, and Intel MKL 10.0. This caused
  minor hiccups due to minor alterations in compiler options, and a problem
  with the unsymmetric matrix solver implemented in the MKL (which affected
  simulations with no Dirichlet pressure boundaries). These problems have 
  been corrected (using structurally symmetric matrices instead), and the 
  code is now functioning normally.


20 MAY 2008:

- For SE/Fourier computations in cylindrical coordinates, an instability can
  develop near the axis due to out-of-plane spatial resolution going to 
  infinity as r --> 0. A filter has been constructed that filters all modes 
  higher than the first mode towards zero approaching the axis. The filter 
  begins acting within the nearest 10% of the domain to the axis.


19 MAY 2008:

- Fixed a couple of small bugs in Tecplot output routine (involving
  sequenced numbering and uninitialised variables), as well as the 
  time stepping routine for SE/Fourier 3D computations. In particular,
  the Dirichlet boundaries were not being carried forward in the 
  transformed sense for cylindrical computations. This has been fixed.


18 MAY 2008:

- Fixed some bugs in the Tecplot plotting routine for SE/Fourier computations.
  Flow fields now closely resemble expected flow fields for test cases.


16 MAY 2008:

- Further work achieved a speedup of approx 7.5 over 16 processors. Still
  attempting further optimisations.

- Parallelization of pressure & diffusion substeps resulted in a speedup
  over 16 cpus of 3.2 times. When part of the advection term was 
  parallelized, a speedup of 6.7 times was achieved. Advection is the
  bottleneck (pressure and diffusion each achieve a speedup of approximately
  14 on 16 cpus, whereas advection currently achieves about 3.3). Work
  is continuing to improve these figures further.


15 MAY 2008:

- Added OpenMP parallelisation of the time stepping loop for the 3D
  spectral-element/Fourier algorithm. Running tests to determine the
  effectiveness of this parallelisation.


14 MAY 2008:

- Updated the way in which Viper reads boundary conditions. All conditions
  are input as mathematical functions, and Viper determines whether they
  are time-dependent, spatially varying, or constant, and computes them
  accordingly.


12 MAY 2008:

- SE/Fourier implementation coming together steadily. Save/load, Tecplot,
  L2, flowrate, monitor capabilities implemented. Still need to implement
  full boundary condition assignment capabilities - currently must read
  in as part of a 2D restart file. Testing of cylindrical solver still 
  required. Forces still required. 

- A rewrite of boundary conditions is being considered, whereby all
  Dirichlet conditions would be read in as mathematical functions, and
  Viper will ascertain which are time-dependent (to be evaluated during
  time integration), and which are time-independent (to be evaluated
  during initialisation). This will be incorporated during the roll-out
  of the SE/Fourier algorithm.


8 MAY 2008:

- Rolled out a provisional version of a spectral-element/Fourier 3D flow
  solver (in both Cartesian and cylindrical coordinates). Still need to
  complete boundary condition assignment, and postprocessing capabilities.


3 MAY 2008:

- Added command ROTATE, which is designed to correct the evolution equations
  for frozen perturbation fields during 2D Cartesian stability analysis for
  the rotating reference frame in which they would effectively be computed in.  


29 APRIL 2008:

- Modified the Tecplot output routines so that the binary files are generated
  by Viper, rather than using the supplied "tecio" library. Viper now generates 
  Tecplot binary files readable by Tecplot versions 10.2 or later. Users 
  with earlier versions should use the ASCII output option (see HELP TECP).


28 APRIL 2008:

- Tecplot 360 (2008) has altered the functions provided in the tecio.dll 
  routine from the previous 2006 version. If users run Viper (which was 
  compiled using the 2006 "tecio" library) while having Tecplot 360 2008 
  installed will be confronted with an error when calling "tecp". I am now 
  contemplating writing my own routines to save Tecplot binary data, which 
  could be made readable from a number of Tecplot versions.


5 MARCH 2008:

- Fixed a Tecplot plotting bug - the velocity divergence (grad_u) was being 
  incorrectly calculated in cylindrical coordinates. This has been fixed, and
  this problem did not affect the solver, which handles divergence correctly.


28 FEBRUARY 2008:

- Added an option "-t" to the Tecplot output routines, which enables ASCII
  output in Tecplot format rather than binary. These files can be read by
  older versions of Tecplot, and may be used by other visualization programs
  more easily.


15 FEBRUARY 2008:

- Fixed a bug in MONITOR, where 2D flows with a w-velocity component active
  were not written with the w-velocity component included.


5 JANUARY 2008:

- Used the modular storage of global variables in the ARPACK routines to save
  the state of an Arnoldi iteration at each call to ARNOLDI. This allows for
  resumption of a simulation at a later point.


4 JANUARY 2008:

- Cleaned up the ARPACK module, and replaced variables stored using the FORTRAN 
  "SAVE" command with variables stored in separate modules.


14 DECEMBER 2007:

- Added a command SAMPLE, which takes the spatial coordinates of a point in the
  computational domain, and outputs to file the interpolated values of velocity,
  velocity gradients, strain rate magnitude and pressure.


29 NOVEMBER 2007:

- Added convergence output to the ARNOLDI command, so that users are given some
  idea as to the progress of the computation.


28 NOVEMBER 2007:

- Added ability to specify a filename prefix for use with the ARNOLDI command
  to allow a unique set of files to be stored, rather than simply the default.

- Fixed a bug with the ia32 versions of the ARNOLDI command, where an unused 
  external timing command was being called as a subroutine rather than a function.
  Removed all calls to this function. 


19 NOVEMBER 2007:

- Implemented an Implicitly Started Arnoldi Method to find the complex eigenvalues 
  (Floquet multipliers), and eigenvectors (Floquet mode velocity fields) 
  corresponding to the matrix operator representing time integration of the 
  linearized Navier-Stokes equations. This is implemented through the "arnoldi" 
  command, which can be used in place of the "stab" command (after a call to 
  "floq") to perform linear Floquet stability analysis.


17 OCTOBER 2007:

- Fixed an issue with the Tecplot calculation of perturbation field vorticity.
  When a fully complex expansion was used, some of the z-derivatives were not
  being correctly computed, which caused havoc for the calculation of perturbation
  field strain & vorticity.


16 OCTOBER 2007:

- Added further simplifications to the function simplifier - now removes redundant
  brackets around variables as well as numbers.


15 OCTOBER 2007:

- A bug in the function parsing code was found, whereby components of variable
  names were being substituted in a function expression by a matching variable
  name - instead, only matching FULL NAMES should be substituted!!! This has
  been fixed, and the displaying of the functions to screen has been streamlined.


11 OCTOBER 2007:

- Found and fixed a frustrating bug in "LOAD" which was causing the w-velocity 
  component of the current-time perturbation field to be incorrectly read.
  The physical quantities were being replaced by garbage. This was also
  rendering attempts to load perturbation fields prior to calling INIT futile,
  as only zeroes were being input, which were promptly normalized to infinity
  when the code attempted to set itself up for Floquet analysis time integration.


10 OCTOBER 2007:

- Fixed a sequential file numbering problem. Calls to "TECP" and "TEC_FLOQ"
  were accidentally both incrementing the same counter within the code.


9 OCTOBER 2007:

- Improved plotting of the "max ds" quantity, the maximum change in the scalar
  field from one time step to the next. It was being incorrectly calculated
  previously, and was therefore always O(1) in size, even when the flow and
  scalar field converged to a steady state.

- Altered output format of time "t" and "max d?" quantities. These now use the
  "G" edit descriptor in FORTRAN, which produces more readable output (i.e.
  only uses scientific notation if numbers are very large or very small).


8 OCTOBER 2007:

- Found and fixed a small issue with the math parser function simplification
  routine. "If" commands (due to the comma) were being incorrectly processed.


7 OCTOBER 2007:

- Shifted searching through elements for Neumann pressure BCs from within time
  stepping code to the initialization code. Thus should result in a small
  speedup in time integration.


5 OCTOBER 2007:

- Scope for a further 10% speedup has been identified for the evaluation of
  high-order Neumann boundaries for pressure. These will be implemented shortly.

- Successfully implemented a speedup in the evaluation of time-varying user-
  defined boundary conditions. If the BCs are not spatially varying, then 
  they are only evaluated once, rather than for every node. This has resulted
  in measured speed increases of 82% in some vortex-interaction simulations.


30 SEPTEMBER 2007:

- Modified the source code varibale naming of w-velocity in the perturbation 
  field, to acknowledge that the w-velocity field is imaginary when a 
  z-velocity is not active in the base flow. This has made the code read more 
  clearly, and now the reduced and fully complex expansions of the Floquet 
  modes have been coded up consistently. Now both are producing the same 
  Floquet multipliers for some simple test cases, in both cylindrical and
  Cartesian coordinates.


29 SEPTEMBER 2007:

- Have attempted a swap of the sign of the change-of-variables in the
  diffusion step for a Floquet mode in cylindrical coordinates with z-
  velocity active - will see if this produces results consistent with
  simulations without an active z-velocity in the base flow.

- Fixed a small bug in the 3D Tecplot Floquet mode reconstruction routine.
  The scaling factor is now consistently applied to the perturbation field
  and the corresponding derivatives correctly.

- Performed some tidying up of the code: The separate eigenvalue/vector
  routine used to calculate the strain fields for ia64 Linux platforms has
  been replaced by the standard code, as the floating-point exception errors
  that were occurring have been avoided by explicitly setting the "-fpe3"
  flag at compile time. 


28 SEPTEMBER 2007:

- Corrected the vorticity components for the "TEC_FLOQ" routine when using
  cylindrical coordinates. The vorticity components are now calculated as
  rotation about each of the cylindrical coordinate axes, not (x,y,z) axes.

- Developed a new routine for outputting the perturbation fields of Floquet 
  modes. The new command "TEC_FLOQ" can be called to generate a three-dimensional
  Tecplot file containing a reconstruction of the vorticity field acquired
  when the 3D Floquet mode is superimposed on the 2D base flow. This command
  works in both Cartesian and cylindrical coordinates, as well as for base flows 
  with or without an active z-velocity component.  


27 SEPTEMBER 2007:

- Stabilized finite-z-velocity Flqouet analysis in cylindrical coordinates:
  problem lay in the change-of-variables applied in the diffusion term. Code
  still giving incorrect multipliers in cylindrical coordinates: testing
  indicates source of error is not likely the advection step. Problem currently
  under investigation.

- Modified the SAVE and LOAD commands to output and read in the 
  imaginary components of the complex perturbation fields when the z-
  velocity is active in the base flow.

- Tests showed that When a non-zero velocity was present in a 3-component
  base flow, the perturbation fields were unstable in cylindrical coordinates.
  Further testing underway.


26 SEPTEMBER 2007:

- Tests of the Cartesian formulation of the complex spanwise expansion
  in the Floquet solver (with non-zero spanwise base flow velocity)
  indicate that Floquet multipliers may be being correctly calculated.
  Further testing underway.


25 SEPTEMBER 2007:

- For some reason, the HELP entry for the SCALAR set of commands was
  missing. This has been updated, and the Viper Reference Manual has
  also been corrected to reflect this change.


24 SEPTEMBER 2007:

- Modified the command interface, so that if macro or loop commands
  are being executed, they are displayed on the command line as they
  are called.

- Implemented complex Fourier mode expansion of Floquet modes if a
  non-zero w-velocity is present in the base flow. Testing is underway.


18 SEPTEMBER 2007:

- Added a command STOPCRIT, which can be used to specify a stopping
  criterion on time integration based on a critical value of "max_du".
  Once "max_du" reaches the specified criterion, time stepping is 
  abruptly cessated, and the control passes to the next input command.


16 SEPTEMBER 2007:

- Missing terms in the perturbation field advection step for Floquet
  analysis of flows with swirl or w-velocity were identified; relating 
  to products of the w-velocity and z-derivatives of the perturbation. 
  These have been added, and tests are underway to verify the solutions.


12 SEPTEMBER 2007:

- Found and rectified a memory allocation-related bug that emerged during
  initialisation for Floquet stability anaysis. This bug was caused when 
  un-initialised CSR matrices were being DEALLOCATED. For some reason, 
  they were being initialised with arbitrary sizes, etc., and this was
  causing run-time errors and segmentation faults across all platforms. 


29 AUGUST 2007:

- The internal use of the "Re" variable in the FORTRAN source files of
  Viper has been rolled over to the name "RKV", consistent with the 
  recent alteration to the name of the corresponding parameter.


27 AUGUST 2007:

- Stuart Cogan identified a bug in the Tecplot plotting of strain for
  axisymmetric flows with swirl (only a 2-by-2 strain matrix was being
  allocated, whereas the code requires a 3-by-3 matrix).

- The correction to velocity vectors along a symmetry boundary has been
  modified to avoid correcting velocities along an axis of symmetry in 
  cylindrical coordinates: this condition is imposed naturally as a
  result of the Galerkin formulation of the problem in cylindrical 
  coordinates, so no further modification is required.


11 AUGUST 2007:

- Isolated and rectified the ia64 crashing when evaluating the 
  strain field. The fix is messy - basically the crash was happening
  whenever the input matrix had zeros in it, so the code performs
  a crude check for this condition, and returns a zero strain in
  these instances, though this is sometimes not appropriate.


10 AUGUST 2007:

- Added an extra timer on the evaluation of Dirichlet boundary
  conditions for the diffusion substep - this contributed about 45%
  to the overall compute time for vortex interraction flows, and
  is therefore a good target for optimization/parallelization.

- Changeover from "Re" to "RKV" when referring to the reciprocal
  kinematic viscosity. This includes "help" files, and user-defined
  functions.


9 AUGUST 2007:

- A number of bizarre errors have been encountered when using very
  large meshes on ia32 linux platforms. These include integer divide-
  by-zero, and others. Some of these may be related to the external
  libraries...

- Neatened the alignment of output from various routines to improve 
  appearance at the command line.


7 AUGUST 2007:

- Migrated ALLOCATABLE TARGET matrices to POINTER matrices in the
  INIT routine when building Laplacian matrices & their inverses. This
  reduces the number of ALLOCATE/DEALLOCATE calls, and should then make
  for faster and more efficient re-initialisation.


4 AUGUST 2007:

- Improved GETMINMAX routine - now can find multiple minima or maxima in each 
  element.

- Rectified a bug in the LOAD routine - a variable "floq_mode" was being used
  before being defined. This was causing occasional crashes.


27 JULY 2007:

- Attempted to locate a memory leak in INIT, with no success. Leak is only
  a problem if many INIT calls are made in a single Viper session, which is
  unusual.


23 JULY 2007:

- Added a "build date" to the title screen that appears when the code is 
  invoked, facilitated by the use of the "fpp" pre-processor "__DATE__" macro.


20 JULY 2007:

- Parallelized elemental Laplacian construction loop.

- Have now parallelized all pressure, diffusion & advection elemental loops.


19 JULY 2007:

- Parallelized an elemental loop at the end of the diffusion substep.
 

11 JULY 2007:

- Added some additional OpenMP declarations to loops in the pressure substep,
  to complement those in the advection step. Presently achieving good speedup
  from 1 to 2 threads, but no further benefit beyond 2 threads.

- Added separate timers for each of the advection, pressure, and diffusion
  substeps.


6 JUNE 2007 (D-DAY):

- Fixed a bug in the "forces" command - when pressure gradients (using "pgrad"
  were imposed on a flow, the pressure component of the force was being
  incorrectly calculated. 


5 JUNE 2007:

- Changed the eigenvector/eigenvalue calculation routines to use LAPACK 
  routine "dsyev" instead of "dsyevr", as the previous routine was problematic
  on the Itanium II architecture on APAC, and the performance benefit in 
  evaluating only selected eigenvectors/eigenvalues was considered
  negligible here, as our matrices are only order 2 or 3 anyway.


2 JUNE 2007:

- Added extra feature to the ADVECT command. Users can now use commands
  "ADVECT ON" and "ADVECT OFF" to switch on/off the calculation of the
  advective terms of the base flow. This could be used for time-varying
  creeping flow computations, with time-dependent boundaries, for example.


29 MAY 2007:

- On Itanium-2 ia64 processors, crashing was observed when evaluating user-
  defined functions as part of the "getminmax" and "int" routines, which
  incorrectly implemented this facility. This has been rectified, and the
  erroneous & unpredictable crashing (with floating-point errors) has been 
  remedied.


28 MAY 2007:

- Ironed out wrinkles in "getminmax" command - the routine had some bugs which
  either resulted in an incorrect turning point being found as a result of 
  polynomial stiffness, or turning points to be neglected completely.  


25 MAY 2007:

- Matt Fitzgerald detected an error in TECP, which relates to the computation of 
  "-o 3" level fields (strain or lambda_2). Problem stems from a "divide-by-zero"
  in strain directional component evaluation.


23 MAY 2007:

- Added a routine that finds and outputs both the physical location, and the value, 
  of a user-specified scalar function. This is invoked by the command "GETMINMAX",
  and has as an intended application the accurate identification of vortices in
  a flow (i.e., by calling /> GETMINMAX 'dvdx-dudy')


27 APRIL 2007:

- The plotting of strain vectors has been improved to more accurately compute vectors
  at element boundaries.


24 APRIL 2007:

- "-o" option added to Tecplot output routine. This option allows different levels
  of output to be selected, from 1 to 3, where 1 is the smallest number of variables
  (and hence the smallest file size) and 3 is the largest.

- Rolled eigenvalue and eigenvector routines being computed as part of the Tecplot
  routine over to LAPACK.


23 APRIL 2007:

- Tecplot output routine modified to fix strain field output. Routine now outputs
  vector components that give the direction of the strain field (eigenvector 
  corresponding to the largest eigenvalue of the strain tensor), and whose
  magnitude equate to the largest eigenvalue - the magnitude of the principal
  strain throughout the flow.


24 MARCH 2007:

- Added extra floating-point overflow and divide by zero checking in "gll.f90"
  and "tec_out.f90" routines to eliminate some run-time errors in linux. Code
  now more robust.


23 MARCH 2007:

- Added a config. file command "gvar_movref", which allows users to
  specify functions for velocity components corresponding to a moving frame
  of reference. The ability for this routine to handle time-varying functions
  makes it far more powerful and user-friendly than the "CHREF" command.


8 MARCH 2007:

- Fixed a problem with "math_parser" module - some code had not been updated
  from the latest fix - some incorrect evaluation behaviour was being generated.

- An error has been cropping up on APAC - isolated its cause as being a floating
  point error in the "find_circle_from_three_points" routine. Error due to argument
  of an inverse cosine being out of bounds. For some reason the Itanium compiler
  did not pick this up.


5 MARCH 2007:

- Dr Stuart Midgley issued a further fix which removed the "Operator following
  operator" error message in "if(" functions.


2 MARCH 2007:

- Made further cosmetic changes to clean up startup and initialization info.

- Implemented a fix from Dr Stuart Midgley for the function parser. This has
  eliminated the improper "Operator following operator" error with "if" functions.


28 FEBRUARY 2007:

- Made some (mostly cosmetic) changes to the code in "ti_advect.f90", which
  in some cases might cause a miniscule speedup in the computations.


27 FEBRUARY 2007:

- Dr Stuart Midgley (author of the Mathematics Parser) sent a fix for the 
  "if(" function error messaging (the correct function was still being 
  evaluated), but another error ("Operator following operator") has been
  produced with expressions of the form "if(x > -1)...". Currently awaiting
  a fix for this.


17 FEBRUARY 2007:

- Fixed a problem with the "loop" command, where long lines of text containing
  functions were sometimes incorrectly parsed because the inverted commas were
  being removed from the text in the scratch file.


16 FEBRUARY 2007:

- A problem with the "math_parser" library has been found - the package is not
  accepting nested "if()" functions - raising correspondence with the developer
  in an attempt to find a solution - otherwise will need to reverse-engineer 
  their source code...

- Removed change of reference shift of perturbation field velocity in CHREF command:
  This was causing irregularities in Floquet analysis. It is incorrect to 
  change perturbation velocity field as the change in reference frame is 
  described by the change made to the base flow velocities. 


24 JANUARY 2007:

- Confirmed with computation of sphere wake at Re=212 that Floquet problem now 
  resolved. Floquet multiplier 1.000 achieved for m=1 mode.

- Problem with axisymmetric Floquet code persisted after base flow axi fix. 
  The cause seems to be the incorrect application of radial and azimuthal 
  Helmholtz equations to the transformed coordinates. A fix of this is producing
  perturbation fields consistent with Natarajan and Acrivos (1993) results.


23 JANUARY 2007:

- Fixed problem with 2D/axi skew-symmetric formulation of advection terms.
  This had crept back in to the solver from an earlier fix. Additional axisymmetric
  terms were missing a 1/2 factor. Checked with wake behind a sphere - drag forces 
  now correctly computed.


22 JANUARY 2007:

- For Tecplot plotting, the code now uses L'Hopital's rule to approximate azimuthal 
  derivatives of perturbation field velocity components on the axis.

- Found likely cause of Axi Floquet problem - an attempt to zero transformed 
  v-velocity values on the axis was using the incorrect global numbering.

- Fixed a bug with the MONITOR routine - users can now perform multiple INIT calls
  in a single session without MONITOR-assigned output files causing the programme
  to crash. Users can also close existing files during a session, and create new ones.

- Identified a problem with the Tecplot routine: Radial component of z-derivatives 
  was not included for cylindrical coordinates. This affected most spatial 
  derivative-based output, including vorticity.

- Floquet test 2: NO phantom pert field vort present in Cartesian (no swirl) 
                  computations with multiple mode numbers.

- Floquet test 1: Phantom pert field vort present in AXI (no swirl) computations
                  with multiple mode numbers.

- Found some Nglobal constants which should have been Nglobal_u: this could have an 
  effect on flows with periodic boundaries...


20 JANUARY 2007:

- Increased amount of information displayed regarding matrix factorisation during
  initialisation phase.

- Detected a problem with perturbation field in axisymmetric computations - appears to
  be associated with Laplacian or Helmholtzian solve.

- Fixed a bug whereby "radial_dist" array was not being constructed for FROZEN Floquet
  analysis - caused perturbation fields to diverge.


10 JANUARY 2007:

- Begun adding !$OMP declarations to parallelise some regions of the code using OpenMP.


5 JANUARY 2007:

- Used IFDEF directives to differentiate Viper for Windows & Viper for Linux in "main.f90".

9 DECEMBER 2006:

- Used "(*,a)" formatting for string output rather than "(*,*)", which was looking
  terrible under the Intel compiler, with all lines breaking after 80-odd characters.
  Now initialisation and HELP facility are again legible! Yay!

- Detected a problem whereby if illegal variables were parsed in function strings,
  the code would obliviously continue with incorrect function evaluations. A firm
  warning message has been added to inform the user when these problems occur. This
  will drastically improve the robustness of the "usr_var" feature of the code.

- Upgraded "math_parser" module from a 2004 version to a 2006 version of the code.


6 DECEMBER 2006:

- Identified and rectified a potential bug in TECP routine in which z-derivatives 
  might have been miscalculated if WVEL and FLOQ were both active.

- Fixed bug in SAVE routine where w-velocity component was not being saved correctly
  due to a use of the near-obsolete "solverID" tag.


3 DECEMBER 2006:

- Fixed an error in Floquet analysis routine where pressure substep was incorrect
  due to incorrect enforcement of a condition on the pressure field for matrix
  consistency which was not required as the perturbation field matrices are
  not singular due to the Helmholtzian-style inclusion of the z-derivative
  component along the diagonal.

- Detected an incorrect elemental vector product operation in L2 routine, which is now fixed.

- Generalised the axisymmetric swirling flow solver to 2D Cartesian flows with non-zero w-
  velocity fields. This will be used for Floquet analysis of vortex interaction 
  studies where the vortices have non-zero axial flow in the z-direction.


17 NOVEMBER 2006:

- Fixed predefined preprocessor symbol declarations for Itanium (APAC) processors.

- Cleaned up some miscellaneous routines & messy statements throughout code.


13 NOVEMBER 2006:

- Detected and fixed a bug which had corrupted the 2nd-dimension elemental derivative
  matrix in 3D only.


12 NOVEMBER 2006:

- In "main.f90", VIPER is now recognized as "Version 3", reflecting the recent
  implementation of the optional scalar transport by advection-diffusion.

- Implemented a low-order scalar advection-diffusion scheme. Scheme uses a 2nd-order
  predictor-corrector scheme for advection of the auxiliary equation used to provide
  the scalar field at the "departure point" of the Lagrangian discretization of the
  advection-diffusion equation. The Lagrangian form of the transport equation
  is a simple diffusion equation, which is solved using a first-order implicit
  scheme. A variable number of advection steps per diffusion update are permitted.


9 NOVEMBER 2006:

- Detected a bug (through IA-64 machine runtime errors on APAC) in the sparse matrix 
  factorization routine. Problem caused by indefinite Laplacian matrices. Have now 
  altered the routines to solve an unsymmetric matrix with a single non-zero entry
  in an extra row, rather than using the previous full extra row & column of ones. 
  Should provide higher accuracy and avoid runtime errors.


8 NOVEMBER 2006:

- Ported code over to Intel Fortran 9.1 on Windows IA-32 machines (previously Lahey
  Fortran Express). Currently having trouble using the TECIO DLL library. Other
  aspects of the solver appear to work correctly. Minor array bounds errors were 
  detected during the migration.


7 NOVEMBER 2006:

- Detected an error in which program exits with untraceable EXCEPTION_ACCESS_VIOLATION
  error on a Windows IA-32 testbed. This occurs regardless of Intel MKL version used, 
  indicating that the Lahey Fortran compiler might be at fault. Tests are underway to 
  identify the source of the error, but this is problematic as the error has no easily
  identifiable pattern. Attempting to install Intel Fortran 9.1 compiler to cross-check.


6 NOVEMBER 2006:

- Added the ability to impose uniform pressure gradients (in x, y & z-directions) via
  the PGRAD command. The pressure gradients are imposed during the advection step,
  and are also included in the estimate of the high-order homogeneous pressure boundary
  condition. This facilitates an easy method for computing pressure-driven periodic flows. 


4 NOVEMBER 2006:

- Removed obsolete routines which provided Runge-Kutta stepping for the initial 2 time 
  steps - now exclusively employing 3rd-order backwards multistep time stepping.

- Modified SAVE and LOAD routines to work with velocity fields over previous 3 time
  steps, instead of only at current time. This avoids the annoying perturbation that
  was added when simulations were restart from saved files due to the solver estimating
  the previous velocity fields for the backwards-multistep time integration.


30 OCTOBER 2006:

- Fixed bug in TEC_BIN routine, where gradient quantities were incorrect due to 
  using the wrong numbering map. 

- Tightened code in static condensation routines: some redundant code remained from
  recent boundary numbering restructuring.

- The unique global node numbering for each variable means that implementation of 
  future variables such as a scalar field is simpler, and no convoluted remapping
  is required in the computations.

- Completed implementation of reorganization/procedurization of boundary condition
  allocation and node renumbering. Still some code to tidy up - no need to shuffle
  between base mesh numbering & reduced numbering when periodic boundaries are present.
  This is now intrinsic to the number maps. Extensive testing required to ensure no
  "global_number_?" or "bmap_?" indexing matrices are out of place. 


25 OCTOBER 2006:

- Began implementation of scalar advection-diffusion field implementation.


24 OCTOBER 2006:

- Implemented user-defined static and transient Dirichlet pressure boundaries.

- Replaced backwards integration estimates of initial 3D velocity fields with a direct 
  copy of current time, as previous method was too unstable.

- Fixed another bug in inter-mesh reinterpolation for LOAD routine. 

- Fixed a bug in the polynomial order reinterpolation routine by reducing ALLOCATE calls.


23 OCTOBER 2006:

- A bug yet to be fixed: diverged solutions develop from loaded (interpolated?) 3D fields.

- Fixed another bug in inter-mesh reinterpolation for LOAD routine. 

- Fixed a bug in the polynomial order reinterpolation routine by reducing ALLOCATE calls.


20 OCTOBER 2006:

- Finished cross-mesh interpolation feature of LOAD routine in 2D - 3D still to be finalized.


18 OCTOBER 2006:

- Partial implementation of an interpolation feature of the LOAD command which allows loading of
  data from a file generated from a different mesh - the data is interpolated onto the new mesh
  providing the file was saved with the "-m" option. Still need to make mapping derivative
  calculations portable.

- Added a "-m" option to SAVE and LOAD routines - this option saves or loads the vector field by
  interpolating the old data and old coordinates onto the new mesh (as per the configuration file).


9 OCTOBER 2006:

- Tweaked screen output so that numbers display more nicely during config & initialization.

- Fixed file sequencing so that SAVE or TECP calls dumping different Floquet modes are 
  independently sequenced.


6 OCTOBER 2006:

- Added a reinterpolation feature to the LOAD command to permit the loading of solutions with
  a different number of quadrature points.


19 SEPTEMBER 2006:

- Fixed a bug in the TRACK SEED routine, where particles were placed at +ve coordinates only, not 
  over whole domain.


18 SEPTEMBER 2006:

- Found & fixed an annoying bug in the Itanium compilation of the io.f90 module. "ifort" did not
  like a WRITE(buffer,*) call, where "buffer" was a CHARACTER string.

- Obseleted "gvar_monitor" and replaced with runtime command MONITOR, which works the same way,
  but allows for multiple simultaneous point velocity monitoring into numerous filenames.

- Modified "TRACK SEED" routine so that particle distribution performed in physical space.


17 SEPTEMBER 2006:

- Fixed a bug in LOAD routine - perturbation fields were not being read from the correct position.


16 SEPTEMBER 2006:

- Fixed a bug in LOAD routine - actual vector fields were not being retrieved for perturbation fields.

- Fixed a bug where the SAVE routine was not defaulting to the "k=0" base flow mode on a SAVE call.

- Removed unnecessary OPTIONAL specifiers for some parameters passed to routines which execute SAVE 
  and LOAD commands - makes for neater code.


14 SEPTEMBER 2006:

- Added "-f" option for filename selection for LOAD routine, and added "-k" option to both SAVE and LOAD
  routines to allow Floquet analysis perturbation fields to be stored and retrieved as well as base flow
  fields.

- Added optional file name selection for time history monitoring. Useful for monitoring the time histories
  of several simulations in the one directory.

- Fixed an error where an ALLOCATABLE array in the particle injector locator routine was not being DEALLOCATED.


13 SEPTEMBER 2006:

- Fixed a bug in the particle tracking routine where injector points which lay near an element boundary
  were occasionally being missed as the search for the elemental coordinates of the injector location 
  only searched the first element found which contained the nearest global node to the injector location.
  If the node is on an edge, face, or corner of an element, several elements need to be tested to find
  the actual location. 

- Fixed a bug in the user-specified variables facility in which runtime errors were encountered when
  no user-specified variables were defined.


7 SEPTEMBER 2006:

- Added a new command "INT", which performs an integral over the computational domain of a user-
  specified function. The function can include time, spatial variables, velocities and their spatial 
  gradients, the Reynolds number, and any user-defined variables from the configuration file.


5 SEPTEMBER 2006:

- Added an option for SAVE and TECP routines which adds a sequential number string to the output
  filenames to facilitate compact macro loops when collecting large series of data files for animation, etc.


31 AUGUST 2006:

- Added a feature to the configuration file "gvar_usrvar", which defines a user-specified variable
  for which a supplied function is evaluated.


19 AUGUST 2006:

- Fixed Gaussian blur of particle densities in Tecplot routine - have now implemented a correct
  expression compensating for the circular/spherical distributions in 2D/3D. The integrals of
  these functions have been verified to be unity in each case.


16 JULY 2006:

- Combined Newtonian and non-Newtonian solvers for consistency. Non-Newtionian computations are
  invoked whenever a non-linear viscosity is entered in the "viper.cfg" file through the
  "gvar_nnvisc" function.

- Added an unmissable warning message to the non-Newtonian solver to clearly identify when zero 
  or negative viscosities are computed, which invariably cause simulation instability.


29 JUNE 2006:

- Added user-defined filename functionality to output of Floquet multipliers.

- Altered header of "flowrate.dat" files so that the boundaries are numbered "bndry1", "bndry2", 
  etc., rather than simply "1", "2", etc. This allows the files to be input directly into
  Tecplot, preserving correct columns.


24 JUNE 2006:

- Resolved an issue where re-starting from a saved axisymmetric flow field was impossible
  due to incorrect treatment of the diffusion terms during prediction of previous time step
  velocity fields.

- Fixed problem with advection of axisymmetric flows, tested on flow past a sphere, and now
  predicting correct drag forces and Floquet mode stability. 


23 JUNE 2006:

- Detected errors in rotational and skew-symmetric forms of advection term in cylindrical 
  coordinates. This has caused an effective reduction in Reynolds number for flow past a 
  sphere.

- Disabled rotational form of acceleration terms for axisymmetry.


18 JUNE 2006:

- Added a function "FREEZE", which cessates time integration of the base flow. This reduces the
  cost of steady-state Floquet stability analysis or particle tracking.

- 2D Cartesian Floquet linear stability analysis implemented, and tested on a circular cylinder
  wake - successfully predicted the emergence of Mode A instability.


16 JUNE 2006:

- Fixed bug in 3D code which resulted in scrambled meshes following recent "%n()" alteration.

- Improved error trapping to accommodate computations with 1st-order elements (no interior nodes).


13 JUNE 2006:

- Implemented advection step for Floquet analysis perturbation fields.

- Modified "max du" reporting during time integration to include perturbation field
  maximum velocity change.


6 JUNE 2006:

- Began implementing Floquet linear stability analysis for 2D and axisymmetric solvers.

- Added "gvar_?" commands to VIPER "HELP" facility.

- Implemented a facility for setting a user-defined initial velocity field in "viper.cfg".


30 MAY 2006:

- Replaced "%n1, %n2, ..." and "%b1, %b2, ..." node number & boundary type 
  identifiers with arrays, allowing some reduction of loop sizes.

- Fixed enforcement of periodic boundaries in axisymmetric computations with swirl.


28 MAY 2006:

- Increased precision of "FLOWRATE", "FORCES" and "L2" routine output to 18 significant
  figures, and added "t" variable to "FLOWRATE".


22 MAY 2006:

- Improved functionality of "FORCES" routine - now permits user to specify a file 
  name to save to, in a manner aligned with functionality of the "FLOWRATE" routine.


19 MAY 2006:

- Added spatial coordinate variables to function for non-Newtonian viscosity.

- Repaired non-Newtonian solver to align with current Newtonian solver.


3 MAY 2006:

- Implemented exact periodic velocity boundary condition (periodic nodes only counted
  once instead of inexact averaging technique previously used).

- Consolidated redundant code in FLOWRATE, FORCES modules, as well as compressed
  Cartesian and axisymmetric diffusion substep routines into one module. Compile time 
  reduced by 75% and binaries approx 10% smaller.

- Consolidated redundant code in symmetry boundary correction routines.


2 MAY 2006:

- Consolidated redundant code in "make_bmap" and "make_bmap_3d" routines so that both 2D
  and 3D are computed using "make_bmap" routine, ensuring consistency in boundary implementation.


26 APRIL 2006:

- Fixed problem with 3D "forces" routine - the viscous terms were being dotted with incorrect
  normal vectors (always normal to 1st face, not correct face).


25 APRIL 2006 (ANZAC day):

- Tests showed that pressure BC was being subtracted instead of added - swapped this around
  with improvement in mass conservation & presumably temporal accuracy.

- Important! Need to check the sign of the pressure BC sfc integral in the pressure
  substep - this might now be negative due to correction of sign of mass matrix elements.

- Fixed bug in Tecplot output routine where compressible vectors were being allocated for
  incompressible 2D flows, causing a crash upon subsequent calls to "tecp".


23 APRIL 2006:

- Coded up global static condensation solver for Runge-Kutta vectors for compressible solver.
  Solver not yet operational - boundary conditions not implemented yet.


22 APRIL 2006:

- Found some errors in the construction of compressible flow matrices, but still results 
  in diverged solution.

- Detected an error whereby mass matrix was -ve due to a recent switch of element coordinate 
  definition. Fixed everything so now elements defined as follows: 
               (p,q,r) --> (vertical, horizontal, height) --> (xi2, xi1, xi3)


18 APRIL 2006:

- Implementing compressible flow solver following recent papers by J. Iannelli. Solver uses a
  diagonally implicit 2nd-order Runge-Kutta algorithm to compute compressible Navier-Stokes equations
  modified to incorporate upwinding in the differential equations being solved by means of an
  "acoustic-convective upstream resolution algorithm". Code compiles but boundary conditions not
  yet implemented.


10 APRIL 2006:

- Modified the Tecplot data output routine to compute 3D quantities on 2D mesh for Axisymmetric
  computations with swirl.

- Implmeneted axisymmetric computations with swirl (non-zero w-velocity evolved with z-direction
  gradients set to zero.

- AXI toggle command changed to a cycling command, allowing selection of three 2D modes:
  2D Cartesian, Axisymmetric without swirl, and axisymmetric with swirl.


9 APRIL 2006:

- Axisymmetric version of 2D solver implemented using cylindrical coordinates. Formulation as per
  Blackburn & Sherwin (J Comp Phys, 2004). Axisymmetry invoked with "AXI" command, and employs the
  x-axis as the axis of symmetry, with y being the radial direction. 


8 APRIL 2006:

- Extended periodic boundaries to 3D - still x-direction only.


7 APRIL 2006:

- Rewrite of boundary definition code successful - new "btag" format now required in "viper.cfg" files.

- Planning to re-design boundary definitions in viper.cfg file: user specifies 2-digit type code for each
  boundary number, with 1st & 2nd digits defining velocity & pressure conditions, respectively. Bndry code
  numbers 1, 2, 3, 4 & 5 are fixed and user-defined static and transient Dirichlet boundaries, periodic and 
  symmetry boundaries, respectively.

- Added periodic velocity boundary condition (restricted to 2D, x-dir only at this stage).

- Added pressure "p" as a variable for transient user-defined boundary functions.


5 APRIL 2006:

- APAC simulations reveal that the current version of the code has greatly improved temporal characteristics on
  an initial-condition-independent computation: Testing the temporal convergence of the shedding period of a 
  cylinder wake.

- Replaced CONSERVE command with ADVECT, which now allows three options for the form of the advection operator:
  Convective, rotational and skew-symmetric. Skew-symmetric is now the default as it conserves linear momentum
  and kinetic energy in the inviscid limit, and it introduces minimal aliasing errors.


3 APRIL 2006:

- Further refined the displaying of infomation to the user as VIPER loads and initialises.

- Fixed bug in rotational form of the advection operator (employed after a CONSERVE command).


1 APRIL 2006:

- Cleaned up configuration and initialisation reporting to screen.

- Implemented a fancier header logo displayed upon program invocation.


31 MARCH 2006:

- Still a problem with temporal convergence tests based on point velocities, etc in the flow. New tests being
  performed based on frequency measurements taken from time history traces of a Karman shedding wake.

- Added a "L2 norm" functionality, where an integral of the velocity magnitude over the computational domain is
  computed. 

28 MARCH 2006:

- Careful inspection of variation of flow quantities over the first few time steps revealled that much of the 
  temporal order issues were due to an effectively 1st-order initial step due to the incorrect previous time step
  velocity fields. Two virtual velocity fields at times t0-dt and t0-2dt have been constructed using an RK4 
  backwards time integration of the advection and diffusion parts of the Navier-Stokes equations. New temporal
  accuracy results to follow.

26 MARCH 2006:

- Tests showed that pressure Neumann BC being incorrectly applied - the gradients were negative of what they were
  intended to be. With correct sign, solution still only first order.

24 MARCH 2006:

- Problems detected with normal boundary vector construction - had become incorrect as a result of Dxi1 & Dxi2
  matrices being defined in wrong directions - this also affected forces/flowrate calculations.

- Found error in pressure substep - was incorrectly imposing Dirichlet velocity boundary conditions during
  intermediate update of velocity field - this might be the cause of the mysterious persisting first-order
  temporal convergence of velocity despite 3rd-order time integration.

- Implemented and tested 3D particle tracking - results appear promising - multi-element simulations correctly
  tracking across elements.

- Cleaned up code in symmetry boundary condition module - compilation time extremely slow in here, so hopefully
  this will help.


21 MARCH 2006:

- Added particle tracking functionality - can now switch off injection of particles mid-computation.

- Modified distribution variance in Tecplot data file output particle density variable - now set depending
  on local density of sample points - the finer the mesh, the lower the variance. This should produce more 
  quality in regions of interest.


20 MARCH 2006:

- Added a Tecplot output variable "particles", which plots a sum of a normal Gaussian distribution of the 
  proximity of particles to each data point in the plot. Currently the standard distribution (diffusion) of 
  the particles is fixed at 0.05 length units.

- Fixed a bug in Tecplot dumping routine where extra data was being added due to old number of nodes being
  passed to TECIO routines after cleanup and reduction of duplicate nodes.

- Implemented 2D particle tracking, and laid the code base for 3D (just need to code up hexahedral element 
  face adjacency information).

16 MARCH 2006:

- Moved looping in "step" command to within the stepping subroutine.

- Implemented a 4th-order Runge-Kutta scheme to advance the advection term for the first three time steps, 
  to accurately fill the previous step velocity vectors required for the 3rd-order backwards multi-step scheme.

- Fixed a small indexing bug in pressure boundary condition evaluation (could this be the 1st-order issue?).


15 MARCH 2006:

- Reorganised allocation, initialisation, and setting or RHS vector routines for time stepping so that the 
  save/load routines function correctly. 


13 MARCH 2006:

- Fixed a bug which resulted in singular global Laplacian matrices when no Dirichlet pressure boundaries
  were present. 

- Added a function "chref", which adjusts the non-Dirichlet velocities to simulate a change in reference 
  frame motion (e.g. a body coming to rest, or accelerating).


12 MARCH 2006:

- Implemented circular arc-based curvature algorithm for 2D elements: still need to implement continuous 
  curvature around closed curves.  


9 MARCH 2006:

- Implemented some minor changes to pressure substep, which should slightly increase speed by reducing 
  number of conditional statements in loops when evaluating pressure Neumann boundary condition.

- Further tests show convergence still 1st-order in time... where is the frakking problem?

- Problem with pressure boundary condition: This Neumann condition was also being "enforced" on Dirichlet 
  pressure boundaries (actually slightly altering the Dirichlet pressure value on the RHS of the equations).
  This might be the reason for the sub-optimal temporal convergence often observed (convergence often little
  better than 1st-order rather than 3rd-order).


8 MARCH 2006:

- Successfully implemented a 3D curvature procedure based on blended circular arc segments. This produces a
  very natural looking curve, and importantly, reduces to perfect circles when one is being represented in 
  the mesh. Tests show that a shape with surface area equal to a circle is described with 8 N=10 elements 
  around the circumference, to the limit of numerical precision. 


6 MARCH 2006:

- Decided to change the element curvature algorithm so that it is based on computed circular arc segments rather
  than polynomial segments. Formulating method now, to be implemented shortly - should get better accuracy
  when computing flow through circular pipes or past cylinders. Was getting flowrates for pipes off by 1% even
  with 8 elements around circular cross-section!

- Changed code so Dirichlet pressure boundaries are also enforced along edges of 3D elements which border
  a Dirichlet velocity region. This way the boundary is more accurately represented.


23 FEBRUARY 2006:

- 10% improvement in time stepping speed through use of global pointer
  storage for time step subroutines, as well as pointer reassignment as
  opposed to vector copies for flow vectors.


20 FEBRUARY 2006:

- Replaced some FORALL statements containing vector notation, which the Intel Compiler is not optimising
  correctly, with DO loops.


13 FEBRUARY 2006:

- By predicting the final CSR vector lengths prior to building them, was able to elimiate many "reallocate_csr" 
  calls, resulting in a slightly reduced memory overhead, and a drastically improved initialisation time (50% 
  improvement measured).

- Added a "trim_csr" routine to "mat_ops.f90" which trims off the unneeded extra elements at the end of the 
  CSR format column and value vectors.


12 FEBRUARY 2006:

- Attempted a speedup of global Laplacian & Helmholtzian Schur complement matrices. Only 1% or 2% improvement.


8 FEBRUARY 2006:

- Altered 3D curvature routine to permit separate curvature of differently numbered boundaries (e.g. better
  for branches in vascular systems - less chance of developing invalid elements).

- Fixed a bug where 3D Dirichlet pressure boundaries were being replaced by interior boundary types along 
  an edge. This was leading to slight errors in computation of pressure-driven flows on some meshes.


25 JANUARY 2006:

- Isolated issue with 64-bit Itanium version of 3D code incorrectly evaluating derivatives. Was due to
  vector notation in a WHERE construct within a FORALL loop. Replaced with another variable and now
  works. 

- With high levels of optimisation the code crashes between "make_bmap" and "make_global_number" routines.
  Attempts to identify the cause of this are unsuccessful - likely either a compiler bug or a problem
  with the external libraries.


24 JANUARY 2006:

- Stripped out all unused routines in "mat_ops.f90", and called previously INTERFACEd routines
  directly. Code now more streamlined, and have removed unused user-defined types from "typeconst.f90"

- Fixed an issue with the 3D solver's load/save routines. It was saving 2D fields due to an 
  incorrect flag setting in "mesh.f90".


20 JANUARY 2006:

- Increased precision of "lambda2" calculation in "tec_out.f90" routine, as intermediate steps can involve
  extremely large numbers in magnitude.

- Fixed bug in "make_bmap_3d" which related to an incorrect indexing of the Dirichlet pressure 
  boundary condition.


19 JANUARY 2006:

- Reverted to solenoidal component of diffusion only in the high-order pressure boundary condition. 
  The explicit treatment of the irrotational part created instabilities which were especially 
  apparent for pressure-driven flows.


18 JANUARY 2006:

- For low-Reynolds number confined flows, an error of up to several percent was detected in mass
  conservation.  This error has been reduced by an order of magnitude by reintroducing the velocity
  divergence component of the diffusion term in the high-order pressure boundary condition.  This 
  is the complete form of the term, and does not require a divergence-free field to be valid.


14 JANUARY 2006:

- Added time derivative term to high-order pressure boundary condition. Improves accuracy of 
  time-varying boundaries.

- Fixed bug in viscous component of drag force routines in 2D & 3D. Terms were ommitted from the shear
  stress calculation. Tested on circular cylinder wake & works.


3 JANUARY 2006:

- Moved the non-Newtonian non-linear diffusion explicit step to before the pressure correction, 
  ensuring that the incompressibility condition is satisfied in the same manner as in the 
  Newtonian solver.


31 DECEMBER 2005:

- Removed the last "TYPE(vector)" statement from the "tec_out" module.


15-31 DECEMBER 2005:

- Replaced 3rd-order Adams-Bashforth advection / Crank-Nicolson with theta mod diffusion with 3rd-order
  stiffly-stable backwards multistep scheme. Testing promising, but issue detected at higher numbers of
  nodes per element - some edge boundary nodes develop a wierd pressure value which perturbs solution.


14 DECEMBER 2005:

- Fixed error in surface integral determination of pressure boundary condition. Now have 2nd-order 
  convergence of velocity fields, but pressure convergence from 1 to 0th order as time step decreases.


13 DECEMBER 2005:

- Fixed small bug in "forces.f90" 2D routine, which arose after Newtonian & non-Newtonian solvers were
  amalgamated.

- Fixed issue with pressure BC (there was a -ve not +ve vorticity in the 2D formulation. Computations now 
  slightly more stable. Consistency as time step reduced needs to be checked.


10 DECEMBER 2005:

- Altered global numbering so that fixed Dirichlet boundaries take precedence over other boundary types, 
  especially symmetry boundaries.

- Altered "save" and "load" functions to store fields elementally. This way the boundary types, etc can be 
  changed between loadings (which changes global numbering) without scrambling flow field data.


8 DECEMBER 2005:

- Fixed bug in 2D "forces" routine for non-Newtonian computations - shear rates now taken from nonlinear viscosity.


6 DECEMBER 2005:

- Implemented a 3D version of the non-Newtonian solver (not yet tested).

- Ironed some wrinkles out of non-Newtonian version of code - small errors in shear rate determination
  and employment of linear reference viscosity. Stability and physical appearance of solutions improved.

- Employed a run-time parsing of relationship for non-linear viscosity as a function of shear rate.

- Added 2D & 3D shear rate variables to Tecplot binary dumps.


3 DECEMBER 2005:

- Reorganised time integration into separate modules - the goal is to be able to switch integration
  schemes at compile time rather than recode the routines.

- Implemented a non-Newtonian time integration scheme, which splits the viscous term into a non-linear 
  part and a linear part. Tests are currently underway.


30 NOVEMBER 2005:

- After testing, derivatives, integrals, 2D & 3D sims with all boundary types working.

- Fixed bug in boundary type categorisation routine - was failing when a user-defined boundary was listed
  prior to a fixed Dirichlet boundary.

- Fixed bugs in the 3D bmap & Jacobian routines which were due either to typos
  in reference texts or from loop optimisations I had implemented.

- Included an alternative time integration module, "time_int_blood.f90", which
  will be used to facilitate non-Newtonian fluid modeling.


29 NOVEMBER 2005:

- Fixed a small bug in tec_bin routine to avoid interpolation errors when the
  interpolation point corresponds to a data point.

- Massive speedup to GLL quadrature point evaluation due to more clever recursive evaluation of Jacobi polynomial.
  Now almost instantaneous. Also removed unnecessary routines in gll.f90 file.


23 NOVEMBER 2005:

- 50% speedup of initialisation phase, 26% speedup of time integration thanks to reducing unnecessary code and
  adding new sparse BLAS commands provided in the Intel MKL 8.0.


2 NOVEMBER 2005:

- Replaced many intrinsic MATMUL calls with BLAS equivalents - a sizeable performance improvement.


26 OCTOBER 2005:

- Formalized "slip-wall" bndry condition to enforce both the tangential flow and zero tangent vorticity.
  This boundary condition has been observed to function well in both 2D and 3D.

- Included tangent vectors in "nml_vect" module, which are employed in "slip-wall" boundary calculations. 

- Assigned "INOUT" INTENT to "p_global" pressure vector in pressure substep routines to avoid crashes
  detected with Dirichlet pressure boundary conditions.


2 OCTOBER 2005:

- Implemented "symmetry" or "slip-wall" (zero normal velocity) boundary condition (number 4).

- Set up routine to pre-compute normal vectors out of each boundary to speed up time integration.


30 SEPTEMBER 2005:

- Fixed bugs with 3D "tec_bin" and time integration routines. 

- Implemented untested lambda_2 data field in "tec_out.f90" routine, in both 2D and 3D.


28 SEPTEMBER 2005:

- Using full storage for all Laplacian & Helmholtzian matrices (except Schur complement) --> much faster initialisation
  and time stepping.


27 SEPTEMBER 2005:

- Regressed to a REAL(DP) pointer array type for velocity/pressure vectors in "time_int.f90".


26 SEPTEMBER 2005:

- Attempted more speedups with no success.


25 SEPTEMBER 2005:

- Simplified vector generation in "tec_bin" routine.

- Approx. 50% speedup by using CSR format derivative matrices for matrix multiplies over CSC format in time integration.

- Slight speedup by swapping index order in "bmap", "global_number" arrays.


23 SEPTEMBER 2005:

- Speedup optimisations to code performed throughout "time_int.f90". Approx. 35%-40% increase in speed.

- PARDISO solve calls now in "time_int.f90" routines for efficiency (approx 5% better). 


22 SEPTEMBER 2005:

- Added variable filename functionality to "flowrate" routine.


20 SEPTEMBER 2005:

- Implemented 2D & 3D static Dirichlet pressure boundary condition (with zero normal velocity gradients).
  This will be useful for modeling pressure-driven flows.

- Rewrote "make_bmap" to conform with simpler "make_bmap_3d" routine: now significant code redundancy.


18 SEPTEMBER 2005:

- Implemented accurate dp/dn boundary condition to ensure mass conservation is accurately enforced.

- Added velocity magnitude data set "vel_mag" to Tecplot output files. 


14 SEPTEMBER 2005:

- Fixed "init" so that solutions can be re-initialised successively.

- Tried a "tecp" speedup to no success.

- Implemented toggle between convective and rotational forms for acceleration terms.

- Identified non-smooth boundaries (i.e. corners) as the culprit for inaccuracy in viscous computations.


13 SEPTEMBER 2005:

- Rewrote diffusion routine: still haven't found why mass not conserved in branched meshes.


8 SEPTEMBER 2005:

- Rewrote 3D "forces" routine - results appear more acceptable.

- "flowrate.f90" produces incorrect results if compiled using "-o2" option - compile using "-o1".

- 3D flowrate implemented - testing incomplete.

- Debugged 2D flowrate monitor routine. Now for 3D...


7 SEPTEMBER 2005:

- Implementing 2D boundary flowrate routines. 2D routine near completion...

- Altered DP double-precision definition to "KIND(1.0D0)" in "typeconst.f90", and tweaked size of "small"
  parameter in "tec_out.f90" to remove annoying artificial internal boundaries from Tecplot data files in 2D.


31 August 2005:

- Sped up "tecp" routine by a factor of approximately 2.5 times.


30 August 2005:

- Fixed problem with user-defined boundaries.

- Detected problem with user-defined Dirichlet boundaries (showing up on reverse face of elements as well)...

- Fixed bugs in face-biased blending component of interior node blending routine!

- Detected bug in interior node location blending - partially rectified, but still buggy.


29 August 2005:

- Removed problem with first global to local spatial coordinate vector transform which was causing 3D curvature
  to be discontinuous across element faces - 3D curvature now working!


27 August 2005:

- Multi-face 3D curvature implemented - still need to improve smoothing across element edges.

- Improved 3D curved face blending: Preserves edge curvature throughout face (good for 
  curvature in a single direction).


26 August 2005:

- Curvature (single-face) implemented. Bugs still present. 


25 August 2005:

- Created vectors to store element, face and boundary surface numbers for curvature 

- Simplified 3D element blending routines based on global indexing in readiness for 3D curvature.

- Built global vertex, edge and face numbering for 3D curvature in "mesh.f90", routine "make_jacobians_3d"

- Changed source code file extensions to standard ".f90" from ".f95" for conformance across platforms.


24 August 2005:
- Fixed tec_out.f95 so that each element is not surrounded by its own boundary.


Pre 24 August 2005:
- Wrote code and made lots of changes.