Change log

31 August 2022
 Added a vertical vvelocity forcing capability to the 2D scalar transport equation. This is required to model a background vertical temperature gradient, and is activated by the "GVAR_SCALAR_VVEL_FORCING" command in the "viper.cfg" file.
10 August 2022
 Removed a text output from the math function simplifier code advising the user of the function being simplified (this had accidentally been left active from an earlier code update and was unnecessarily verbose).
27 May 2022
 Fixed a bug in the function simplifier that was failing to capture duplicated negative signs when a negativevalued constant was being subtracted, if the instance was positioned at the end of the function string.
13 May 2022
 Fixed some logic in the command interface of the "LINE" and "SHEET" commands.
22 November 2021
 Added a new command, "KE_EVOL", that calculates and outputs to a text file the terms of the kinetic energy evolution equation. The KE evolution equation is constructed by taking the dot product of velocity with the momentum equation. If buoyancy or electric potential fields are active, the respective gravity and Joule dissipation contributions are included in the output. This command is particularly helpful as it includes the domain integrals of the viscous and Joule dissipation terms.
17 November 2021
 Added sensitivity field (as per equation 8.11 from Giannetti & Luchini, J. Fluid Mech. 581, 167197, 2007) to the output from the "ENDO" endogeneity routine for completeness.
14 November 2021
 Fixed a bug in the threecomponent (nonzero outofplane base flow velocity) linearised adjoint solver (cylindrical and Cartesian versions) associated with a scalar (temperature) transport term that migrates across to the adjoint linearised momentum equation. There was an indexing error located in those terms since a recent code maintenance exercise (where the two and threecomponent cylindrical and Cartesian linearised time integration routines were separated into four separate routines.
11 November 2021
 Completed upgrading of the "ENDO" endogeneity feature. Now fully supports both Cartesian and cylindrical coordinate systems, plus endogeneity contributions arising from linear forcing terms in the momentum equation, Boussinesq buoyancy, and the Lorentz force in the quasistatic MHD solver.
8 November 2021
 Upgraded the "ENDO" endogeneity feature to handle cylindrical coordinates in addition to Cartesian, and endogeneity contributions from (Boussinesq) buoyancy and Lorentz force.
6 November 2021
 Implemented the endogeneity calculations as per Marquet & Lesshafft (2015), activated using the driver command "ENDO". Saved 2D base flow as well as forward and adjoint eigenmodes are required to execute this routine. Currently this is only implemented in Cartesian coordinates.
21 October 2021
 Added viscous and Joule dissipation fields as an output variable in the Tecplot output command "TECP". This is activated using the "vars dissip" option.
21 October 2021
 Added the electric current density field as an output variable in the Tecplot output command "TECP". This is activated using the "vars j" option.
 Removed the ASCII format feature and the interpolation onto uniform grid feature from the Tecplot output routine.
29 June 2021
 Fixed a bug that had crept into the threecomponent linearised solvers due to a complex constant representing i × wavenumber accidentally being assigned a real type, hence reducing to zero, thus screwing up outofplane derivatives.
28 June 2021
 Fixed a minor numbering error in the 3component Cartesian linearised solver involving the diffusion correction for the real part of the wvelocity field.
18 June 2021
 Added a command "ADD", that adds userspecified mathematical expressions to velocity and/or scalar fields (base flows or a specified linearised perturbation field.
12 May 2021
 Implementing a variant of Barkley's eigenvalue solver whereby converged eigenvectors are subtracted from the Krylov subspace and perturbation fields with the goal of converging larger numbers of eigenmodes if required (avoiding the situation where the dominant mode swamps all other modes in the perturbation field).
4 May 2021
 Fixed a logic error that was incorrectly prohibiting transient growth analysis via the "TG" wrapper if an analytically expressed base flow was being used.
21 April 2021
 Modified behaviour of Barkley's eigenvalue solver  each iteration, eigenvalue estimates are now written to file (sorted by largest magnitude), along with their respective convergence residuals, and a statement on whether each had converged subject to the supplied tolerance.
 Modified the reporting to screen of the Barkley variant of the eigenvalue solver to be a little more streamlined and to report each of the residuals rather than the lessuseful differences in eigenvalues each iteration, as they can shuffle around depending on the convergence behaviour.
 Altered the code that prepended the Tecplot output filename with "ioerr_" in the event the originally requested filename was busy or unable to be written to. Previously Viper endlessly attempted to prepend the filename, which could result in an endless loop and voluminous display output if the system was constraining writing of any files. Viper now only attempts a single prepend. If that fails, no Tecplot file is written and control is returned to the Viper command line.
16 April 2021
 Fixed a typo affecting linearised solver evolution of scalar perturbation fields due to an accidental power operator (**) appearing instead of a product operator (*) in the scalar perturbation advection update.
15 April 2021
 Added an option to the "CORIOLIS" command that varies the strength and sign of the Coriolis term in 2D Cartesian simulations conducted in a rotational reference frame. The "beta" parameter takes values ranging from 0 to 1, where 0 corresponds to a flat geometry having a constant prefactor to the Coriolis term, through to 1 modeling a latitudinal variation of the Coriolis term where the south pole lies a radial distance of 2 units from the origin.
1 April 2021
 (Not an April Fool's joke) Switched the generic linearised perturbation field time integration routine to four separate routines dedicated to specific combinations of cylindrical or Cartesian coordinate system, 2component or 3component base flows. The intention is to reduce the amount of conditional branching in the code making it easier to maintain/modify, and perhaps to reduce the memory footprint or overhead of the routine.
31 March 2021
 Added a command %quot;FMLF" (Fundamental Mode Linear Friction), that invokes a linear friction term applied to the fundamental Fourier mode of Cartesian spectral elementFourier 3D simulations only. This could, for instance, force the fundamental mode to behave as though subjected to linear friction similarly to the SM82 model for quasi2D MHD flow, for instance, while still permitting the 3D modes to evolve with unimpeded mechanics.
21 March 2021
 Added a "viper.cfg" command "GVAR_CIRCULAR" that takes x and ycoordinates, and deforms into a circular arc about this coordinate any element edge in the mesh that has a constant radial distance from the supplied coordinate. This can be useful for achieving a mesh that fully conforms to a circular boundary, which was not the default functionality when using Viper's automated boundary curvature via "GVAR_CURVE".
15 March 2021
 Fixed a dividebyzero error in "ELECVTX" and "ELECVTXGRID" commands when vortex lay precisely on a mesh node.
13 March 2021
 Added command "CORIOLIS" to apply Coriolis acceleration to 2D Cartesian simulations to model rotating reference frame.
 Replaced the code used to compute the exponential integral function Ei(x) used by "ELECVTX" and "ELECVTXGRID" with a routine from the GNU Scientific Library 2.6, which is approximately four to five times faster.
12 March 2021
 Updated "ELECVTX" and "ELECVTXGRID" so that the electrically driven vortices are based on a quasi2D averaging of Sommeria's (1988) viscous pointelectrodedriven vortex velocity profile solution, rather than applying the series solution in the confined domain  this update has been made to overcome the excessively slow calculations of large summations at every point on the mesh that the series solutions required.
5 March 2021
 Added commands "ELECVTX" and "ELECVTXGRID", that invoke the forcing associated with an electrically generated vortex or a chequerboard grid of electrically generated vortices in a quasi2D SM82 MHD flow in a rectangular domain.
2 March 2021
 Picked up an "already allocated" error when a spatially varying forcing function was being used. This has been corrected.
23 January 2021
 A bug in the code that constructs the shared global elliptic operators for various active flow fields was being tripped by an unusal combinations of boundary conditions that happened to see the u and wvelocity fields sharing one operator while the vvelocity field used a unique operator. This has been fixed.
21 January 2021
 Fixed an allocation error in a vector containing reciprocal of radius recently introduced that obstructed the cylindrical solver.
20 January 2021
 Renamed the "floquet" module and corresponding .f90 source file to "linpert" and "linpert.f90".
 Refactored variable naming in the source code to better disambiguate usage of "field", "mode" and "wavenumber".
 Corrected a missing imaginary unit in one of the nonBoussinesq advection term zderivatives.
19 January 2021
 Replaced some divides by sum of quadrature weight vectors by multiplication by precomputed reciprocals of those integrals to streamline code.
 The time integration routines has a legacy quirk whereby advection RHS vectors were evaluated in one set of vectors before being copied to a second set. These have been streamlined into a single set of vectors to reduce memory footprint and eliminate unnecessary vector copy operation.
17 January 2021
 Removed "ITERATE" command, which complicated the 2D time stepping code, and didn't generate a benefit for accuracy or stability.
 The time integration routines has a legacy quirk whereby advection RHS vectors were evaluated in one set of vectors before being copied to a second set. These have been streamlined into a single set of vectors to reduce memory footprint and eliminate unnecessary vector copy operation.
16 January 2021
 Corrected a memory allocation error created by a new loweroperationcount global gather routine introduced in the 14 January 2021 build.
14 January 2021
 Implemented the optional centrifugal buoyancy contribution in the thermal transport equation.
8 January 2021
 Added a command "BOOSTCONV" implementing the stabilisation algorithm developed by Citro, Luchini, Giannetti & Auteri (J. Comput. Phys. 2017), that can be used to drive a naturally unstable solution towards a steady state, or to accelerate the convergence of a steadystate solution.
 Added a command "ELECVORT", that adds an oscillating electrically generated vortex into a quasi2D MHD duct flow (based on an equispaced placement of electrodes in a duct).
4 January 2021
 Modified the "L2" command so that the reported L2 norm is calculated as the square root of the domain integral of the square of the magnitude of velocity. Previously the result was returned without the square root.
 Added a command "RESID", that computes the square root of the domain integral of the time integration residual for the current time step (i.e. the absolute difference between the flow fields at the current and previous time steps).
23 December 2020
 Implemented nested OpenMP parallelisation in the construction of the elemental Laplacian operators, and rewrote the code used to generate the operators. This has significantly simplified that part of the code during initialisation, and should make near optimal use of parallel threads where available.
11 December 2020
 Disabled the "SVV" spectral vanishing viscosity feature, which was not being used and was making the construction of elliptic operators less efficient.
9 December 2020
 Implemented significant speedups in the routines used to find the element and parametric coordinates corresponding to a given physical coordinate, which underpins the "SAMPLE", "LINE" and "SHEET" commands. E.g. for a 10,000point line on a mesh with 800 elements with polynomial order 23, the code previously took approximately 41 seconds to execute the required iterative searching and interpolations. The improved code executes the same operation in 10 seconds on a single thread, with solid speedup over multiple OpenMP threads (e.g. taking just 1.5 seconds with 16 threads).
6 December 2020
 Fixed an issue with the "LINE" and related commands ("SAMPLE", "SHEET", etc.) where points on element interfaces could be missed if the polynomial order was high.
5 December 2020
 Rewrote the duplicate elemental elliptic operator identification code that has produced a significantly faster initialisation time.
4 December 2020
 Revised the "MAKE_BMAP()" routine to improve efficiency. Large 3D hexahedral meshes were taking several seconds to scan through the mesh for face matches.
2 December 2020
 Replaced separate evaluation of derivatives (e.g. du/dx, du/dy, du/dz, dv/dx, etc.) with a unified "gradient" evaluation routine, that facilitates speedups by invoking a single instance of OpenMP parallel threads, and by avoiding the repeated evaluation of derivatives with respect to the parametric coordinates within elements. Efforts to achieve further efficiency gains using a single 3*NperElem by NperElem sparse matrix to evaluate the three parametric derivatives in 3D (2*NperElem by NperElem in 2D) was thwarted by what appears to be a bug in the MKL routine that executes sparse CSR matrixfull matrix multiplication.
 Introduced a set of global mass matrix vectors based on the various global numberings. These avoid the more costly division by the inverse global mass matrix that was being used in various places to conveniently integrate a globally numbered flow field variable.
27 November 2020
 An alternative derivative evaluation strategy, whereby flow field variables on all elements are packed into a single array and their products with the sparse elemental parametric coordinate derivative matrices are formed, from which the spatial derivatives in physical coordinates are recovered via the mapping vectors. This was found to be consistently faster than the other approaches (one benefit is that it eliminates a loop through elements, therefore replacing many small matrixmultiply calls with few large calls  this may have benefits for improved parallelism). This approach has been implemented in the 3D hexahedral solver.
 Testing of speed of derivative evaluation methods permitted a memory and timesaving threshold to be enforced, whereby full and alignedfull matrixvector derivative evaluation methods are removed from the timing tests during initialisation for large jobs, as they are orders of magnitude slower than the tensorproduct approaches.
26 November 2020
 Reordered the derivative evaluation tests so the typically faster methods are tested first  this allows slower methods to be cessated early, which can expedite initialisation time in large jobs.
12 November 2020
 Extended the "LINE" command to the 3D hexahedral solver.
 Implemented a new command, "SHEET", which interpolates a surface between 4 supplied coordinates in a 2D or 3D domain, and interpolates the flow solution and gradients on that surface, similarly to the "SAMPLE" and "LINE" commands.
7 November 2020
 Rewrote the global matrix assemly routine. The new routine uses a compacted storage approach that is more streamlined and performs more than 100 times faster for large 3D hexahedral jobs than the previous routine.
 Reduced the time taken to perform the derivative evaluation tests during initialisation, which was becoming cumbersome for very large jobs.
3 November 2020
 Fixed the nonBoussinesq implementation in the spectral elementFourier 3D code (the nonlinear buoyancy update was being computed in Fourier space rather than physical space by accident).
 Implemented some OpenMP parallelisation to the static condensation code for Poisson and Helmholtz solves  speedup is reasonable to several cores for the 3D hexahedral code.
14 October 2020
 Refined an estimate of the size of the global boundaryboundary sparse matrices, which reduces the number of memory reallocation calls invoked during initialisation. This has introduced a significant speedup of initialisation, particularly in 3D hexahedral jobs.
25 September 2020
 Rewrote linearised perturbation field advection term code, for fields comprising the full complex form (activated when a wvelocity is active in the base flow  "WVEL"). The revised forward and adjoint term code uses complex flow field vectors rather than separately evaluating the real and imaginary parts, to reduce the prospect of errors in the various cylindrical/Cartesian, 2component/3component, forward/adjoint combinations.
 Working on a steadystate solver based on the Uzawa algorithm. Not yet debugged.
23 August 2020
 Corrected a bug producing rapid divergence in the Boussinesq adjoint solver, where a term appearing in the adjoint thermal transport equation (originating from the buoyancy term in the momentum equation) had the wrong sign.
 Fixed a bug in the implicitly restarted Arnoldi eigenvalue solver, called with the "arpack" option of the "LSA" command. This rarely encountered issue arose when multiple perturbation fields were being concurrently computed, and one happened to modify the number of eigenvalues being sought, which then corrupted the restart file loading for the next perturbation field. Further error checking has also been added to the "LSA" restart file loading routine.
19 August 2020
 Added a missing data structure to the "VARTG" restart file that was throwing an error if a job happened to resume during a backtracing iteration (where underrelaxation factor is reduced to maintain monotonic increase in growth estimate).
4 August 2020
 Renamed the vorticity components in Tecplot output to be in terms of the components rather than the plane on which the rotation acts.
26 July 2020
 Added calculation and outputting of the turbulent dissipation rate to the "TIMEAVGRS" command. This is implemented for the 2D quad, 3D hex and spectral elementFourier 3D solvers.
25 July 2020
 Implemented commands "FLOOR" and "CEILING", used to enforce minimum and maximum values on any of the velocity or scalar fields during time integration.
26 June 2020
 Updated the functionality of "INTF" so that the integrated function is evaluated in physical space, before being Fourier transformed, at which point the absolute value of the Fourier coefficients of each mode are integrated over the domain and written to file. This is helpful to get an idea of how much each Fourier mode contributes to the 3D field being evaluated.
12 June 2020
 Fixed an error with the assembly of the electric potential global Schur complement of the Laplace operator matrix.
 Fixed an unallocated vector issue with the 3D hexahedral solver.
12 May 2020
 Corrected the cylindrical version of the SEF3D implementation of the "TIMEAVGRS" command so that the avg(v) and avg(w) components are no longer converted from cylindrical to Cartesian.
11 May 2020
 Corrected implementation bugs (unallocated variables, etc.) with the "TIMEAVGRS" command.
9 May 2020
 Added options to the "TIMEAVGRS" command to enable the user to specify the number of zdirection interpolation planes for Tecplot output from a spectral elementFourier 3D simulation, and to decide whether they want Viper to explicitly calculate the perturbation product (Reynolds stress averages for Tecplot output). These options permit the size of the Tecplot files to be controlled. The Tecplot output code for this command has also been restructured to reduce peak memory consumption.
 Added an option ("d") to the "BUOYANCY" command that specifies a Froude number to act on the scalar convection term similarly to how the "c" option specifies a Froude number to act on the advection terms. This is currently implemented in the 2D quad, 3D hex and SEF3D solvers, but not the linearised perturbation field solver.
5 May 2020
 Added a command "BNDRYSAMP" that outputs flow field data along a specified boundary.
26 April 2020
 Added a command "TIMEAVGRS" that computes time averages of products of velocity components in addition to the velocity components themselves, allowing Reynolds stress terms to be computed.
22 April 2020
 Added the quasistatic MHD / electric potential field / Lorentz term to the hexahedral 3D solver.
20 April 2020
 Fixed an error in the job resumption functionality of the "VARTG" command, where previously the resumption would fail if the restarted job happened to correspond to a stage in the iterative scheme where a backtracing reduction in the underrelaxation factor was required.
14 April 2020
 Altered the behaviour of "TIMEAVG" so that it saves a restart file each time it is called allowing timeaveraging jobs to be run over long multirun time integrations.
12 April 2020
 Added ":ALTERNATING" option to the "ADVECT" command, which activates an alternating convective and divergence form of the advection operator. These are the two parts of the skewsymmetric form of the operator, which conserves both linear momentum and kinetic energy in a collocation discretisation, whereas the default (convective only) form conserves neither. This is implemented in the 2D Cartesian and cylindrical solvers (with or without the wvelocity field being active).
27 March 2020
 Extended the new "FREEZE VEL" and "FREEZE S" capabilities to the spectral elementFourier 3D solver.
26 March 2020
 Altered the behaviour of the "FREEZE" command. This now accepts parameters "VEL" or "S" to freeze only the velocity or scalar fields. Omitting the parameter freezes both fields.
17 March 2020
 Added "efirst" and "efloor" options to the "VARTG" command for variational transient growth analysis. These options are used to specify an initial value of the "epsilon" underrelaxation factor in the iterative scheme, and the minimum value it can be adaptively reduced to, respectively.
12 March 2020
 Added command "ZNF", that toggles on/off an imposition of zero net flowrate (ZNF) on a 2D (zerowavenumber) linearised perturbation field. In some circumstances, such as comparison between a linearised evolve of a 2D disturbance in a periodic domain and a nonlinear evolve with "FORCEFLOW" imposing a fixed net flow rate, this can be helpful to ensure the same constraints are imposed in both cases.
10 March 2020
 Updated the nonlinear transient growth iterative scheme to include a backtracking feature where the "epsilon" variable is reduced by a factor of 4 if the current value does not result in an increase in the growth estimate.
5 March 2020
 Fixed a bug in the spectralelementFourier 3D Tecplot output "TECP" routine, where some velocity gradients were being corrupted due to incorrect initialisation, in turn corrupting the vorticity fields.
3 March 2020
 Fixed a bug affecting one of the derivative evaluation routines for yderivatives having complexvalued input, in recent optimised versions of the code. This was due to a complex/real data type mismatch.
 Added a test of the accuracy of complexvalued derivative evaluation methods to trap any bugs that might affect derivative calculations in spectral elementFourier 3D simulations (previously only the realvalued routines were tested).
2 March 2020
 Fixed a bug with the "SAVE" command when called from spectral elementFourier 3D simulations and the "s" sequential file numbering option activated, where the numbering array was being incorrectly indexed, possibly giving an incorrect first number in the sequence.
26 February 2020
 Fixed a bug screwing up the fundamental mode of the outofplane velocity field in the spectral elementFourier 3D solver (it was becoming zero when it shouldn't have). Because the fundamental Fourier mode is real only, Viper was branching to the real (rather than complex) version of the elliptic solver called to solve the viscous Helmholtz equation for the viscous diffusion substep, passing temporary real vectors rather than the usual complex vectors. Reverting to complex solves for all modes has eliminated this bug. The issue may have been due to length mismatches between the general temporary vector and the vectors containing each of the different flow fields.
 Fixed a bug in the nonoptimised build of the code that manifested in the spectral elementFourier 3D code due to an unused allocatable array that was unallocated was passed into a multipurpose subroutine.
24 February 2020
 Corrected a fault with the "TIMEAVG" command where the Tecplot output file was not being written when called from a spectral elementFourier 3D computation.
6 February 2020
 Added "SAVE" and "TECP" output from "VARTG" iterative solver for the solution state at the beginning of the forward pass (u_0).
31 January 2020
 Fixed the input parser for a couple of Viper's commands that were not accepting mathematical expressions evaluated at calltime for consistency across all other commands. This included "DERIVTYPE" and the "E" option of the "VARTG" command.
 Added "SAVE" and "TECP" output from "VARTG" iterative solver for the solution states at the end of the forward pass (u_T) and the end of the adjoint pass (v_0) before the new guess of the optimal initial disturbance (u_0) is calculated.
29 January 2020
 Corrected "FORCEFLOW" forcing output from spectral elementFourier 3D solver.
28 January 2020
 Extended the "FORCEFLOW" capability to the spectral elementFourier 3D solver.
21 December 2019
 Corrected an error with the "VARTG" algorithm that was screwing up attempts to resume a partially converged run.
20 December 2019
 Added the elemental mapping Jacobians to the output of "TECP" if called prior to "INIT". Previously, only the mesh information was written. This is largely for curiosity's sake, but may be of interest in curvilinear and/or highly distorted meshes, as a constant Jacobian within each element is desired for accuracy of the mappings. If this reveals any issues, these may be able to be addressed via modification of the placement of quadrature points within offending elements.
 Added checking for invalid elements via a test for negative Jacobian values during startup. This is not foolproof  elements can be invalidated outside of the discrete quadrature points at which the Jacobian is explicitly evaluated, but it is a good trap for aggregious errors that can emerge from bugs in mesh generators, etc.
19 December 2019
 Fixed an array out of bounds bug in the implementation of Barkley's eigenvalue solver where the number of requested eigenvalues was greater than the current iteration (and size of the Krylov subspace).
 Testing the spectral elementFourier 3D "VARTG" solver is underway.
16 December 2019
 Extended the "VARTG" command to the spectral elementFourier 3D solver.
 Corrected a sign error in the adjoint equation implementation in the 2D quadrilateral "VARTG" solver.
10 December 2019
 Extended the "VARTG" command to the 3D hexahedral solver. The spectral elementFourier 3D solver implementation is pending.
 Added a command "VARTG", that performs a nonlinear transient growth analysis based on a variational method. This is currently implemented in the 2D quadrilateral solver, and will be extended to the 3D hexahedral and spectral elementFourier 3D solvers in a future update.
4 December 2019
 Added a command "BASE", that extracts current flow state in a 2D quadrilateral / 3D hexahedral / 3D spectral elementFourier simulation and converts the evolving flow to a disturbance to the extracted base flow. This also sets Dirichlet boundary conditions for the disturbance to zero. This can be used with "varTG" to conduct a nonlinear variational transient growth analysis.
29 November 2019
 Added forcing terms to the scalar field in the spectral elementFourier 3D solver.
20 November 2019
 A bug was identified that affects linearised perturbation fields with nonzero wavenumber and spectral elementFourier 3D simulations, only in Cartesian simulations with no Dirichlet pressure boundary conditions applied. This has now been corrected.
19 November 2019
 Replaced array range references "A(:)", "B(:,:)", etc. in the source code with wholearray references, e.g. "A", "B". This can offer the compiler greater opportunities to optimise the code as array ranges and whole arrays are treated differently.
15 November 2019
 Implemented a command "varTG" that performs a transient growth analysis using a variational method similar to that described by Pringle, Willis & Kerswell (J. Fluid Mech. 702, 415443, 2012). The analysis can be conducted either on two or threedimensional linearised perturbation fields, or on twodimensional nonlinear perturbation fields (invoked using "NONLINEAR").
9 November 2019
 Implemented an alternative Krylov subspacebased eigenvale solver as an alternative to the Implicitly Restarted Arnoldi Method implemented via the ARPACK package Viper has used previously. The alternative method follows the scheme proposed in Barkley, Blackburn & Sherwin (Int. J. Numer. Meth. Fluids 57, 1435–1458, 2008), and appears to offer favorable convergence performance. It is now set as the default eigenvalue solver in the "LSA" driver routine, etc. The ARPACK approach can be invoked by adding the flag "arpack" to the driver routine call in Viper.
8 November 2019
 Added a text output buffer option "b" to the "SAMPLE" command, to give users the option of buffering data for a number of calls before doing a less frequent bulk output of data to file. If this is successful, it may be worth extending this to other output routines.
1 November 2019
 A bug was detected when using "reconload" to reconstruct a flow field, due to the recently added "gvar_analytic_?" functionality, that was accidentally zeroing the velocity field components every time step.
30 October 2019
 Added a command "NONLINEAR", which when combined with "PERT" and a zerowavenumber perturbation field, adds the nonlinear advection term to the evolution of the perturbation field. This changes the linearised solver to a nonlinear solver, and is akin to decomposing a full 2D flow into a stable base component and a fluctuating part.
28 October 2019
 Changed the behaviour of the "VERBOSE" command in the verbose builds of the code  now it acts simply as a toggle, switching verbose output on/off. A "GVAR_VERBOSE" command has been added to the "viper.cfg" command list to extend verbose output to the parsing of the configuration file.
 Separated the 2D, 3D hexahedral and spectral elementFourier 3D time integration routines into separate modules for easier navigation while editing/debugging.
 Created a source code file "stability_drivers.f90" that now contains the driver code formerly in "cmd_interface.f90" for the linear stability analysis, linear operator SVD from leading eigenvalues, linear transient growth, etc, capabilities. This will be more convenient for the new minimal seed nonlinear capability currently under construction.
16 October 2019
 Corrected a missing premultiplication of userdefined Neumann boundary conditions by radius when in cylindrical coordinates.
14 October 2019
 Fixed an error with the linearised adjoint solver due to an uninitialised indexing vector during the pressure substep that crept in during work on Viper's backend in September 2019.
13 October 2019
 Fixed a divergence error in the linearised solver caused by uninitialised base flow vectors being passed to the perturbation field time integration subroutine when "FREEZE" was called.
10 October 2019
 Removed the "SYMMETRY" boundary condition type  an old, inexact attempt to enforce a stressfree condition on an arbitrarily curved boundary.
5 October 2019
 Corrected an error that had crept in from 30 August 2019 build in the advection term of the axisymmetric cylindrical solver with swirl ("WVEL" active).
 Implemented the adjoint linearised perturbation field evolution equation for scalar transport and Boussinesq convection (2 and 3component solvers). This involves the exchange of a term from the scalar transport equation to the momentum equation when a scalar field is active, and a byproduct of the buoyancy term alters the scalar transport equation when buoyancy is active.
4 October 2019
 Extended "OVERINT" capability from 2component/2D velocity advection to 3component/2D velocity advection and their scalar transport convection counterparts.
 Unrolled 2component/2D velocity advection, 3component/2D velocity advection and 2D scalar convection calculation routines back into the main 2D time integration routine. Also did the same for the linearised perturbation field advection/convection calculations, and the 3D hexahedral advection/convection calculation routines.
2 October 2019
 Fixed an error in the math function simplifier where the function fragment immediately following a boolean operator was accidentally being skipped over.
18 September 2019
 Source code housekeeping: Removed dimension2 localtoglobal numbering arrays in favour of the "indx_el_to_?()" indexing vectors introduced more recently in the code's life.
29 August 2019
 Added commands "gvar_analytic_u", "gvar_analytic_v", "gvar_analytic_w" and "gvar_analytic_s" used to specify analytic expressions for a 2D base flow (it is expected that the expressions are timedependent, otherwise the user could use the "gvar_init_field" capability along with the "freeze" command.
 Replaced a large number of large character string data structures used to parse content in the "viper.cfg" with a single temporary structure.
24 August 2019
 As part of ongoing elemental/gather changes, replacing "NperElem*(e1)+1:NperElem*e" index syntax in loops over elements with preevaluated range "N1:N2".
 Replaced "radius_el(1:Nelemental)" with "local_y(1:Nelemental)" which it was a duplicate of to streamline the source code.
 Completed implementation of the variable magnetic field direction, low_Rm MHD capability in the 2D Cartesian twocomponent solver, the 2D Cartesian threecomponent solver, their corresponding linearised perturbation field solvers, and the Cartesian SEFourier 3D solver.
23 August 2019
 Continuing to work on modification to code relating to changes to elemental vs globally gathered storage of flow field vectors.
 Found and fixed a bug in the v and wvelocity component implementation of Neumann boundary conditions.
16 August 2019
 Changed the time integration algorithm so that the advection RHS is not gathered globally before the pressure substep to avoid error propagation from element interfaces where the discontinuity from the advection derivatives across elements leads to wiggles in the intermediate veloicity field whose divergence forms the RHS of the pressure Poisson equation. The gather is delayed to after the Poisson RHS velocity divergence operation is computed. Results of this modification are pending.
 Removed the "OUTLET" boundary condition type as it was experimental and is not used.
14 August 2019
 Changed the "MHD" command inputs so that users can supply the magnetic field vector.
8 August 2019
 Implementing a feature that extracts highest mode from 1D basis within elements as a possible marker for regions of underresolution in a solution.
27 June 2019
 Switched sign of pressure gradient term in adjoint of linearised solver to align with published work contrasting the form given in Barkley, Blackburn & Sherwin's 2008 works.
23 May 2019
 Used more generic preprocessor symbols to more consistently output compiler version and build date information when Viper opens.
 Added an error message to alert users of the Windows build of the code that the spectral elementFourier 3D algorithm is not currently implemented (requires a serial replacement of the MPI_ALLTOALLV() routine).
 Removed the "xHost"/"/QxHost" option from the optimised compilations of the code to avoid possible issues if an executable is not run on the same CPU architecture.
21 May 2019
 Modified the endoffile loop exit when reading "viper.cfg" file, so that a line of text is parsed correctly when it is the last line of a file, which previously triggered a loop exit as endoffile reached when line was read, rather than after contents of the line were interpreted.
23 April 2019
 Removed a "STOP" call from the "fn_u_ut()" subroutine in the 2D Nusselt number calculation, that was previously activated when a point was identified as being outside the domain. In some cases (such as an obstacle in a channel), this was inappropriate.
19 April 2019
 Isolated all MPI routines to the msg_passing.f90 module, which is now the only module to call the MPI library. This module contains driver routines that are called by the code and which fall back safely for singleprocessor/nonMPI execution of the code. Note: I have not yet made the spectral elementFourier 3D Windows MPI routines nonMPIsafe.
5 April 2019
 Added a command "VERBOSE" to the verbose build of the code. This command restricts verbose output only to routines specified by the user. This is useful to reduce output file size when diagnosing issues.
13 March 2019
 Tidied up the 2D force calculations, by taking constant factors out of the vector summations.
 The 2D/SEFourier 3D "forces" calculation routine only calculated outofplane shear forces when "is_w_vel" was active (via the "wvel" command in 2D simulations). This meant that in spectral elementFourier 3D simulations with a net outofplane flow, the nonzero shear force would not be output. This has been corrected.
21 February 2019
 Removed duplicate output from each MPI process mesh is being checked for similar element shapes.
 Found some DZGEMM() statements in an unused spatial derivative subroutine that supplied DBLE not DCMPLX parameters. This has been fixed in case the routine is called in a future code update.
 Added a missing conditional statement during initialisation to avoid an unnecessary wvelocity matrix operator copy operation during initialisation.
20 February 2019
 Fixed a divergence error in the 2D solver activated when "advect off" is called and a scalar field is active, due to an uninitialised vector holding the calculated advection term.
23 January 2019
 Modified the logic of the streamfunction outputting facility in the "TECP" routine. Previously the streamfunction was only output for 2D/2component flows, but it is also valid for 2D/3component flows.
17 January 2019
 Improved checking and output during factorisation of global sparse matrices to detect rare cases where a mesh has no homogeneous global nodes (e.g. one element with all Dirichlet boundary conditions).
 Fixed a bug in SAMPLE where a MERGE() function was leading to an unallocated array "floq_modes()" being addressed when linearised fields were not set up. Fixed by reverting to a conditional statement.
15 January 2019
 Refined the degrees of freedom calculation displayed during startup (Dirichlet boundary nodes are removed so only nodes containing values forming part of the solution are included from each field).
14 January 2019
 Fixed some typos in Viper's "Help" facility.
 Added an option "minmax" to the "line" command, which can be used to output the mimimum and maximum values of the function along the interpolated line (using quadratic curve fitting for a highprecision estimate of the turning points.
3 December 2018
 Established that the test for matching elements of elemental diagonal mass matrices ("We(:)") in addition to elliptic operator entries was redundant, as the elliptic operator already scales with "We(:)" entries along the diagonal. Thus the test for duplicate elliptic operators has now been streamlined to remove the mass matrix test.
 Tzekih detected an error in the duplicate elliptic operator identification routine, where slight differences could see one element recorded as a duplicate of an earlier element, but a later element only shows up as a duplicate of the earlier duplicate, not the original element. This is now fixed by removing duplicates from the tests against later elements.
23 November 2018
 Updated Viper's title banner displayed when the code opens to say "Viper" rather than "Vmpir".
19 November 2018
 Completed compaction of elliptic operator matrices. Currently one compaction is used for all operators (pressure, viscous and scalar diffusion). However, it is possible that extra compaction might be possible for the pressure Poisson operator as it doesn't have the additional components on the diagonal; there is also scope for having elements of different construction (e.g. defined by vertices starting bottom left, top right, etc, if interior nodes are also counted in an ACW direction.
15 November 2018
 Beginning compaction of elliptic operator matrices to remove duplicates for more efficient elemental operations through gathering (matrixvector into matrixmatrix multiplications, solving linear problems for multiple simultaneous RHS rather than repeated singleRHS solves).
11 November 2018
 Remaining outstanding issues: 2) Chris also reports some reloads producing slightly different subsequent time integration, 4) Still need to get to the bottom of the time accuracy degradation in the linearised perturbation solver. Scope for further speed increases lies in duplication of PARDISO solve calls for complex fields  can replace these with an nRHS=2 call with real and imaginary parts in two double precision RHS vectors. May also be possible to achieve greater performance in tensorproduct derivative evaluation if we use DGEMM() over two elemental dimensions rather than Nqdpts DGEMV() calls. May be worth grouping elements with same Laplace/Helmholtz matrices, and conducting operations with DGEMM() and multiRHS DPOTRS() rather than DGEMV() and singleRHS DPOTRS(). Will also reduce memory requirement and data transfer to/from CPU each time step.
 Noticed that "alpha" and "beta" coefficients in new "dzgemv()" routine for computing a product of a doubleprecision matrix by double precision complex vector multiplication being used in the linearised solver and the spectral elementFourier 3D solver were being passed as double precision reals, not double precision complex numbers. This may be stuffing up the routines  checking now. Yes, this has fixed it  both the linearised solver and the spectral elementFourier 3D solver now give stable computations, and SEF3D solver tests correctly against a 3D Kovasznay flow test case.
10 November 2018
 Corrected a bug in the function simplifier where a nested bracket pair was incorrectly identified across separate bracketed terms. This was picked up during a test of Kovasznay flow, where the function parser reported an error from 8 June 2018 builds through to corrected 10 November build.
 Fixed periodic boundary issue. This was due to a logic fault in a loop that was used to merge pairs of boundary nodes  an overly eager incrementing index sometimes skipped past nodes, meaning that at worst no pairs were merged. The correction to this then found an error in the subsequent loop that deflated global node numbers so that the missing numbers from the periodic node merging weren't skipped  this too has been fixed.
9 November 2018
 Outstanding issues needing rectification: 1) Chris reports periodic boundaries not being properly enforced, 2) Chris also reports some reloads producing slightly different subsequent time integration, 3) New static condensation approach seems to have screwed up the complex static condensation routine, 4) Still need to get to the bottom of the time accuracy degradation in the linearised perturbation solver.
 Replaced the choice of three matrixvector product operations chosen based on speed with a single solution strategy for the static condensation working around the PARDISO() sparse matrix solve for the global Schur complement arising from the boundary nodes. The previous approach required two matrixvector multiplies; one using the inverse of the elemental innerinner parts of the Laplace/Helmholtz operators, and the other using a precomputed product of the inverse and an off diagonal part of the elemental operators. As inverse matrixvector products are less precise than the equivalent linear solve opreation, the code now replaces the two matvec products with one cleaner matvec product and elemental matrix solves.
26 October 2018
 Updated PARDISO "iparm(:)" vector to start from zero as per current documentation. Stuck with default values for all parameters, except calling for output of FLOPs and number of nonzeroes.
 Slightly streamlined output relating to sparse matrix construction/factorisation.
25 October 2018
 The function simplifier had an incorrect implementation of the "hat(x)" function when preevaluating functions. This has been corrected.
16 October 2018
 Saw negligible performance improvement with POINTERtoALLOCATABLE modification.
 Fixed the divergence issue with the linearised perturbation fields: this was due to an incorrect matrix being supplied in a DGEMV() call in a branch of the code only activated if the DGEMV() routine was timed to be faster than the intrinsic MATMUL() and BLAS DSYMV() routine.
14 October 2018
 For performance improvement, began converting POINTER arrays to ALLOCATABLE arrays wherever possible.
9 October 2018
 Revised the divergence checking routine to use ISNAN() intrinsic rather than the cruder IF x /= x check in the hope that this eliminates the compiler optimization bug throwing incorrect divergence detections.
7 October 2018
 Restructured the divergence checking routine which appears (hopefully!) to have corrected that wierd compileroptimizationrelated bug where divergence was being incorrectly activated.
 Fixed an output formatting issue with the "int" routine for perturbation fields.
 Revised divergence detection routine so that it (a) reports the field as well as the processor exhibiting divergence for perturbation fields, and (b) an MPI "OR" operation is used via an MPI_ALLREDUCE() command to inform all processors tha divergence has been deteced to facilitate a graceful shutdown.
1 October 2018
 Merged two loops into one for efficiency in each of the global node renumbering routines.
 Fixed an issue with the renumbering for cylindrical Dirichlet notes arising from the recent speed increase updates.
29 September 2018
 Added a new command, "pert2", used to create a sequence of wavenumbers for linearised perturbation fields (linear or logarithmic).
 Removed the acceptance of the obsolete command "floq" in place of "pert" for establishing linearised perturbation fields.
27 September 2018
 Profiled code on a large job with 16x16 elements, finding that time spend in BLAS "dsymv" routine was 10 times that spent in "dgemv" routine, despite the calls appearing together in the part of the code that does the innerinner part of the global Poisson and Helmholtz operations elementbyelement. "dsymv" was costing approximately 25% of the overall compute cost. Implemented code to test the speed of calculation of FORTRAN intrinsic MATMUL(), and the BLAS DGEMV() and DSYMV() routines, along with branching to select the faster method (and this is MPIsafe: different processes can choose what is best for them). Unfortunately the timing process gives markedly different results from one run to the next, and the difference between the BLAS routines is seldom anything like the 10fold difference profiled  it may be that the routine is also called by the sparse solver or another external routine. Certainly the PARDISO() sparse solver consumes a decent amount of the overall compute time.
26 September 2018
 Found a bug in "floquet.f90" in relation to MPI when multiple perturbation fields are simultaneously evolved over multiple processes. The indexing vectors were being incorrectly constructed. Fortunately, this capability has not typically been employed by users; users preferring to run many singleCPU jobs rather than a single large parallel job.
25 September 2018
 Implemented faster code for numbering the global boundary nodes, particularly the code responsible for merging the pairs of nodes shared across periodic boundaries. 100fold and better speedups have been achieved, streamlining and reducing the precompute time as the code opens before accepting command input.
 Added extra checking to the "LOAD" command so that the code no longer blindly attempts a physical interpolation between meshes just because LOAD is used with "m". It now also confirms that the mesh informtion is actually in the file.
24 September 2018
 Removed the MPI broadcast call when "freeze" has been called that spreads the base flow solution across MPI processors every time step. This removes an MPI bottleneck at the time step level.
 Fixed an issue where nonVerbose versions of the code could crash due to the "step" index not being defined when it was used for text output reporting on the state of a time integration at the end of a block of time steps.
 Misc. code cleanup: Removed unused internal "SolverID" flag; added output once the "Ndimension" variable has been set while reading a mesh file; replaced the user defined TYPE for mesh vertices and mesh elements/boundaries with simple integer arrays.
20 September 2018
 Corrected displayed text in "load" command when element number doesn't match while "m" option active.
 Corrected minor text outputs in command interface module that were out of date.
22 August 2018
 Added a command FIXSCALAR, which corrects a scalar field each time step to maintain a specified domain integral value.
 Corrected the "stopcrit" command, which was not working in the nonverbose version of the code in 2D simulations.
18 August 2018
 Shifted divergence checking routine so that divergence is checked every time step. Also streamlined the procedures so they should be quicker.
 Altered the check for infinities to be of the form "TRUE when number x = x1.0", which also traps numbers with very large magnitudes once the roundoff error exceeds unity. This is used both in the "neat_real_number()" function to avoid trying to format an infinity, and the divergence checking routines.
16 August 2018
 Corrected a bug in the "neat_real_number()" function that involved an overflow error in an internal write operation used to construct a formatting string for a subsequent write operation.
14 August 2018
 Fixed an uninitialised variable bug with logical "isSignLeft" variable in the function simplifier that was rarely triggered.
8 August 2018
 Simplified "neat_real_number()" function used to format printed floatingpoint output to avoid a rare record overflow error.
7 August 2018
 Corrected a bug arising from 6 August domain interpolation routine fix. Routine now tests all nodes in mesh from nearest outward until a hit is found. It does not rely on a restricted subset of nodes.
6 August 2018
 Fixed a bug in domain interpolation routine whereby the fixed search through nearest 5 elemental points could fail if there were more than 5 elements sharing the nearest global mesh node. Changed the code to test only the nearest global node, but all elemental nodes it shares.
25 July 2018
 Cleaned up a blank line in "init" output in the nonverbose version of the code.
6 July 2018
 Corrected a bug detected by Zhi in the Cartesian form of the perturbation kinetic energy evolution calculation involving part of the dissipation term.
15 June 2018
 The previous upgrade dated 13 June 2018 permitting functions/variables to be supplied as floatingpoint inputs to Viper commands has been extended to integer inputs (the evaluated floatingpoint value is rounded to an integer at the time that the command is called).
 Corrected some bugs involving +/ signs and addition/subtraction operations in the function simplifier code.
13 June 2018
 All floatingpoint values of options/parameters supplied with commands on the Viper command line are now parsed as userdefined functions (of time "t" and "RKV" and userdefined variables) rather than fixed floatingpoint values. Note: these functions are evaluated at the time the command is called.
8 June 2018
 Added a function "saw(x)" to the math parser. This produces a periodic "saw" wave, defined as "saw(x) = xfloor(x)", which features a line rising from 0 to 1 every unit increase from integers along the xaxis.
 Rewrote parts of the function simplifier routine with clearer syntax and logic to avoid illegal/incorrect simplifications.
7 June 2018
 Added a function "hatsmth(x, xs)" to the math parser. This produces a smoothed hat function, where the width of the smoothing applied to the discontinuities at x=0.5 is supplied by "xs". This is useful for applying regularised boundary conditions.
6 June 2018
 Fixed a writetoscreen bug that occurred when a perturbation field of zero wavenumber was specified (thrown by the "neat_real_number()" function when number was out of bounds, e.g. infinity or NaN).
5 June 2018
 Added "floor()" and "ceiling()" functions to the maths parser; also upgraded the function simplifier to operate on many more functions.
 Added a compiler flag "Verbose", and "ifdef Verbose" statements to only compile in text output on various debuggingrelated outputs when the compiler flag was set. The goal is to have a leaner compiled code if this is adopted more widely, rather than relying on conditional statements, etc.
21 May 2018
 The horizontal Nusselt number routine was being found to sometimes converge endlessly. This appeared to be caused by the halving tolerance at each level outpacing the convergence of the integrals, which floored out at around 1e19. To avoid this, a convergence trap was included which cessates iterations if convergence is below 1e15.
14 May 2018
 Tightened the output of the "load" routine to remove unneeded blank spaces now that percentage counter is no longer used.
10 May 2018
 Reduced volume of screen output in the available potential energy calculation during the Tecplot "tecp" routine. Percentage status update replaced by simple status bar.
 Similarly streamlined the verbose "% complete" output from the "load" routine.
9 May 2018
 Corrected a bug in the recursion floor implementation for adaptive GaussKronrod and adaptive trapezoidal routines where the specified floor value was being ignored.
 Modified the screen output of the "save" routine to reduce the amount of text output. Large 3D jobs can result in hundreds of data fields being saved to file. Rather than a line of text and a percentage complete value being written to screen for each one, the code now writes a single "." to screen for each field as a rudimentary status bar.
5 May 2018
 Added a floor to the recursion in adaptive integration routines for horizontal 2D Nusselt number calculations. This is specified as a positive integer with the "m" option in the "nu_horiz_2d" command.
4 May 2018
 Revising NewtonRaphson interative scheme for finding parametric coordinates corresponding to given physical coordinates due to the code seemingly excessively throwing a "Warning: Divergence detected in physical to parametric coordinate conversion routine." message. It was found that convergence routine was flooring out at distances in the order of 1e8 to 1e9, so the original inbuilt convergence criterion of 1e9 wasn't always being triggered before this noise floor was reached. Once at the noise floor, some iteration steps were slightly larger than their predecessor, leading to the divergence warning. Solution has been to modify the convergence check for consecutive larger steps over two iterations instead of one, and to relax the convergence threshold to 1e7 from 1e9.
 Found that the 2D horizontal Nusselt number routine can get stuck in an infinite loop if a dividebyzero happens in evaluation of either a bulktemperature or local Nusselt number calculation, which throws a NaN that never gets recognised as convergence. Added traps to set the result to zero if these circumstances arise, though it shouldn't happen provided there is a net nonzero flow, and a temperature difference between the wall and interior.
26 February 2018
 Fixed a bug in "ONLYW" command so that it is correctly applied to 2D/axi velocity fields.
15 February 2018
 Added a command "ONLYW" which switches off the u and vvelocity fields in 2D and axisymmetric computations so that flows with only a nonzero wvelocity field can evolve without growth of instabilities associated with nonzero u or vvelocities.
19 January 2018
 Peyman detected an error with the 18 October change to how floatingpoint numbers are written: small negative numbers (having negative signs before the significand and the exponent) were the same size as the space allocated to display them, which then removed any spaces between numbers output to screen or text file. This has been corrected in the present release of the code.
18 October 2017
 In rare instances where the magnitude of the exponent of a floatingpoint number exceeds 99 (e.g. extremely small nonzero values), the Fortran number formatting for screen or data file output was producing a number format that was not recognisable as a floatingpoint number in postprocessing software (e.g. Tecplot, Excel, etc.). This has been fixed by specifying the number of character spaces used to write the exponent part of the number as 4 up from the default 3.
24 August 2017
 Fixed an indexing error in the centrifugal buoyancy terms in the linearised solver leading to segmentation faults.
28 July 2017
 Added zdirection force components to "FORCES" command output, which can be nonzero when there is flow in the outofplane direction.
27 July 2017
 Corrected the implementation of the centrifugal buoyancy terms in the linearised solver.
9 July 2017
 Corrected an error in the "nu_horiz_2d" command, where the integrated Nusselt number was not being divided by the horizontal length.
8 July 2017
 Added an option "arnoldi" (or "a" for short) to the "LSA" command, which replaces the forward linearised equations with their adjoint, which are integrated backwards in time. This can then be used for sensitivity analysis of the flow.
27 June 2017
 Corrected an error with the implementation of forcing of the form "G*u" that occured only when evolving the adjoint linearised equations of a perturbation field  it was incorrectly being subtracted rather than added.
24 June 2017
 Corrected a bug in domain integration routine ("int", "l2") that resulted in infinity being output for linearised perturbation fields having zero wavenumber (as the integral over the 2D plane is multiplied by the outofplane wavelength of the perturbation. Instead, a unit (Cartesian or 2*Pi (cylindrical) multiplier is used for zerowavenumber fields.
10 May 2017
 Removed the debugging text output from the Arnoldi eigenvalue solver routine.
8 May 2017
 Fixed a typo in the linearised quasistatic MHD electric potential field calculation that was throwing a segmentation fault due to an attempt to access an unallocated indexing vector.
3 May 2017
 Modified "STEP" to accept negative integer input. A negative number of time steps performs backwards time integration (for evolution of the adjoint linearised governing equations). It should not be used otherwise.
28 April 2017
 Fixed an error in the implementation of the quasistatic MHD equations in cylindrical coordinates, where the magnetic field is assumed to act in the axial direction.
18 March 2017
 Added option "o" to FORCEFLOW command, facilitating output of the forcing applied at the current time forcing.
17 March 2017
 Added a new functionality to the evolution of 2D scalar fields; a forcing involving multiplication by the uvelocity component is facilitated by invoking the "GVAR_SCALAR_UVEL_FORCING" command in the "viper.cfg" file. This command is useful for computing heat transfer in 2D channels where the background horizontal thermal gradient has been removed from a computation (meaning that the scalar field represents a fluctuating temperature).
9 March 2017
 Added a new command, "FORCEFLOW", which fixes the flowrate through a specified boundary (2D only).
 Added a new command, "SPREADSCALAR", which rescales scalar fields to avoid them collapsing to a constant value when only one temperature is imposed on a boundary while other boundaries are insulating.
4 March 2017
 Modified behaviour of "SAMPLE" so that for linearised perturbation fields with a nonzero wavenumber, the zcoordinate is also used to permit extraction of the outofplane variation in the perturbation field.
 Wrapped the I/O "reread()" routine in an MPIsafe code to ensure that it is only ever called by the master process.
28 February 2017
 Updated the quasisteady MHD solver (2D; LSAl; SEF 3D) to handle computations in cylindrical coordinates  in these simulations the magnetic field is taken to act in the axial direction.
27 February 2017
 Bug detected in SEFourier 3D simulations where a scalar initial field was not being set up if there was no corresponding initial velocity field specified (i.e. "gvar_init_scalar_field" being called with no corresponding "gvar_init_field" call).
31 January 2017
 Added Neumann boundary condition evaluation for the scalar field in spectral elementFourier 3D simulations. NOTE: Neumann boundary conditions are NOT CURRENTLY IMPLEMENTED for the velocity field, only the scalar field.
18 January 2017
 Upgraded the "line" command so that it works in SEFourier 3D simulations when the "a" option is called (which generates average values along the line rather than values at all points along the line, and thus is still meaningful in a 3D domain).
1 December 2016
 Modified the output from Tecplot when the buoyancy option is activated. Previously the code would output a text update for every point on the mesh when computing the APE density. This made for very large and slow Tecplot file computation when the mesh is large. This output has been modified to output a maximum of 100 updates during this calculation.
30 November 2016
 Added two new options to the "line" command: "u": used to specify a userdefined function to be evaluated along the line (if this is called, only the values of this function are output rather than the full set of flow variables and their gradients. "avg": when called the routine outputs only the average of the flow/function values interpolated along the line.
8 November 2016
 Fixed a file output bug in the underlying procedure behind the "int" and "l2" commands where file and screen output was being skipped when the commands were called from spectral elementFourier 3D simulations.
7 November 2016
 WARNING: The cylindrical coordinates solver has been broken since 4 May 2016. There was an error in the vvelocity and wvelocity Helmholtz matrices that are used for the viscous diffusion substeps. This error was inadvertently introduced as part of the Spectral Vanishing Viscosity implementation in June 2016. This has now been fixed.
 Added more detailed output when divergence is detected; the code now reports which specific field(s) have diverged.
4 November 2016
 Corrected the behaviour of the "load" routine when the "r" is specified. Previously if the loaded flow was a base flow, the time and time step were being changed to whatever their stored values were from the loaded file, which is undesirable behaviour. The "r" option is used to add a stored flow field to an existing solution  the time and time step should not be changes.
 Added a centrifugal buoyancy capability to the 2D quadrilateral, 3D hexahedral, linearised and spectral elementFourier 3D solvers. This is activated by specifying a value for the new "c" option in the "buoyancy" command.
25 October 2016
 Added extra text output during invokation of the "load" command. This is to better inform the user of how the code has interpreted the options given when calling the command.
29 September 2016
 Fixed the "int" routine so that integrals are evaluated correctly for perturbation fields. This requires transforming to physical space, evaluating the integral in physical space, and taking the mean of that evaluation to recover the integrated value. Note that this means that the integrated value will be zero for linear functions of perturbation fields with nonzero wavenumber, because these fields oscillate around zero in the spanwise direction. Nonlinear functions integrate to nonzero results, such as the energy of a perturbation field (u^2+v^2+w^2), etc.
5 July 2016
 Fixed a bug in the enforcement of timedependent Dirichlet boundary conditions on scalar fields in the 2D solver.
4 May 2016
 Fixed some typos in the output file header at om the "line" command, where cylindrical coordinate system variable names were appearing in place of their Cartesian counterparts in Cartesian simulations.
6 April 2016
 An allocation error was found with the quasistatic MHD part of the linearised solver. This has been fixed (an unallocated vector of spanwise waveunmbers was being referred to).
 Fixed a bug where problems would arise if PERT was called before WVEL.
17 March 2016
 Changed the code that interpolates the flow fields from previous times steps to new time intervals when the time step is changed. Under some circumstances (such as when "wvel" is called but before "init" allocates the wvelocity vectors) not all fields have yet been allocated. The routine now explicitly checks that all fields are associated before attempting to access them.  Added the quasistatic MHD functionality to the linearised perturbation field solver: this enables linear stability analysis and transient growth analysis of quasistatic MHD flows.
29 February 2016
 Found the linear stability analysis error: the error was due to an incorrect additional term that had appeared in the advection calculation for the imaginary wcomponent of the perturbation field. The additional term that appears in this calculation in cylindrical coordinates was not being correctly switched off in the Cartesian solver.
16 February 2016
 Modified the way that the zderivative of pressure was being calculated for the pressure substep of linearized perturbation field time integration (was previously being calculated in a loop through each element, with dp/dz numbered using global wvel numbering taking the value I_IMAG * wavenumber * local pressure coefficients). Now uses the same elementbyelement loop, summing the weighted contributions before dividing by the global weight matrix.
10 February 2016
 Zhi has detected a possible fault with the linear stability analysis algorithm. The fault appears to lie with the solution of the linearized NS equations. The code gives the same results if using STAB or LSA, so it is not due to the Arnoldi eigenvalue solver implementation, but likely resides in the time stepping routine. He has found for a circular cylinder test case that quasiperiodic eigenvalues are significantly different to those computed by Blackburn et al. (2005), Barkely & Henderson (1996), and he has tested on the same mesh size as their studies. He is running further tests on a square cylinder test case, which can be compared to results from Robichaux, Balachandar & Vanka; Blackburn et al.; Sheard, Fitzgerald & Ryan (2009). ##### Fixed on 29 Feb 2016 #####
 Removed the zerod/dz linearised NavierStokes solver capability from the code (it was never being used anyway)  trying to clean things up in order to isolate the aforementioned error.
8 January 2016
 Corrected a Tecplot output error from spectralelementFourier 3D simulations where the pressure was not output correctly if a userdefined field was included in the Tecplot output.
 Fixed the "int" output from spectralelementFourier 3D simulations, which had not been implemented since the MPI upgrade.
19 December 2015
 Altered behaviour of ENERGIES command to avoid large text output when command is called frequently and/or the number of mesh points is large.
27 November 2015
 Added a command ENERGIES, that calculates kinetic energy, potential energy, background potential energy and available potential energy integrated over the 2D/axisymmetric plane (note for SEFourier 3D simulations, these are calculated on the zerowavenumber fundamental component of the solution only).
26 October 2015
 Added APE conversion production terms from base flow thermal gradients to the linearised perturbation field energetics output.
11 June 2015
 Added a "buoyancy" (accepted abbreviation "b") to TECP command. This option outputs either the available potential energy (APE) field if TECP is called to output a Boussinesq base flow simulation, and if outputting a linearised Boussinesq perturbation field, this option includes 7 terms from the out ofplane fluctuation kinetic energy equation in the Tecplot output: the energy conversion from potential to kinetic energy in the fluctuation, and the 6 nonzero outofplane averaged fluctuation energy production terms from the base flow.
31 May 2015
 Added "tic" and "toc" commands to display timing information to the user emulating the behaviour of the corresponding commands in MATLAB.
27 May 2015
 Added a command "nu_xsect_2d" which is used to compute the Nusselt number for a 2D outofplane channel flow with a heated righthand wall.
 Added a command "nu_horiz_2d" which computes the Nusselt number for a 2D horizontal channel flow with a heated bottom wall.
26 May 2015
 Added a command "nusselt" which is used to compute the Nusselt number for a 2D channel flow (heated bottom wall, horizontal channel walls).
19 May 2015
 Intiated implementation of quasistatic MHD time integration algorithm.
4 February 2015
 Fixed an instance where if multiple "init" calls are made in an SEFourier 3D simulation, the code attempted to allocate an already allocated array.
10 December 2014
 Added a command "timeavg" which can be used via successive calls to construct a time average of the flow field.
31 July 2014
 Tzekih detected an "output statement overflows record" error during time stepping for the linearized Boussinesq NavierStokes solver used for linear stability analysis. This has been corrected by making the string used to store the text prior to displaying to screen.
16 April 2014
 Added the spatial gradients of the scalar (temperature) field as allowable variables for the integrand expressions for the "int" and "intf" commands.
 Added the spatial gradients of the scalar (temperature) field as allowable variables for the userdefined field variable (specified using the "u" option in the "tecp" command.
 Added scalar gradients to the Tecplot output from SEFourier 3D simulations.
12 November 2013
 Fixed output from the "LINE" command, where "z" and "r" were being written to the output file header for Cartesian simulations, instead of "x" and "y".
*** Note: Tony has encountered crashing during SEFourier 3D MPI simulations (infrequently, and on the VLSCI system). This appears to be due to multiple processes attempting to simultaneously write to the output file. This has not yet been fixed ***
4 November 2013
 Upgraded the "RAND" command to work for all flow fields (2D, 3D, linearised perturbation fields), in addition to the spectralelement Fourier 3D solver.
19 August 2013
 Found the error with Boussinesq computations. This was due to an incorrect scalar field (intermediate not endoftimestep) being used at a point in the time integration algorithm. This was especially damaging for simulations using cylindrical coordinates because there was a 1/r factor difference between the two fields.
16 August 2013
 In an attempt to track down the buoyancyrelated error that has crept into the code in the last few months, checking compilation of the code on different systems and with different compilers: 1) Using different versions of the Intel Fortran compiler on NCI Vayu and Raijin produces the same results for Kovaznay 2D flow test case. 2) Yet to check if buoyancy effects same on both compilations 3) Compiling the code on Vayu with GNU Fortran compiler found the following: a) Errors in "msg_passing" module if compiled as a nonMPI code b) Nonstandard use of .XOR. rather than .NEQV. for exclusiveor logical operations in the recent bandwidthminimisation experimental code in the messagepassing section of the code.
15 August 2013
 Detected a potentially serious error with the use of index mapping vectors of the form "indx_x_to_?", where statements of the form "f_x(:) = f_?(indx_x_to_u(:))" were not seeing the LHS being correctly filled with the content from the RHS. The same also applied to the "indx_el_to_?" arrays. This affected some 2D base flow derivatives, and several outputs, including the L2 norm. The L2 norm is now fixed, and all other affected areas have been replaced with more robust loops through all entries in the vectors. **** However, note that the Boussinesq buoyancy computations are currently incorrect *****
6 August 2013
 Fixed an issue where the buoyancy coefficient may not be set to a default zero value if the user did not explicitly set a value when calling "buoyancy".
24 June 2013
 Fixed a bug in the LOAD routine where interpolating a saved solution onto meshes with a different polynomial order was failing due to a vector already being allocated  this had crept in during the MPI upgrade.
 Reverted to the behaviour where executing the LOAD routine changes the time step to that stored in the file. This is required as the file contains fields at the previous time steps, that would be incorrect if the time step is not set to be that in the saved file.
4 June 2013
 Fixed a bug in the "line" command, which will now correctly output outofplane component of velocity and its derivatives if the wvelocity field is active in a 2D or axisymmetric simulation.
 Fixed a bug with the "rotate" option recently added to the "tecp" command, where it was not correctly defaulting to a zero rotation angle if no rotation angle was specified.
3 June 2013
 Changed "LOAD" routine behaviour  "dt" and "RKV" are no longer changed to values stored in the restart file, only "t" is taken from the file  the other parameters need to be set in the "viper.cfg" file.
 Fixed an output bug in the "flux" commands where multiple MPI processes might each try to write to the same output file, despite only one thread doing the calculation.
30 May 2013
 Added a command "LINE", which extracts flow field data along a line from a 2D simulation.
24 May 2013
 Added an option "rotate" to the "tecp" command that rotates 2D or 3D hexahedral meshes clockwise by a specified angle in degrees about the origin in the Tecplot output files. This is useful for aligning visualized results with a specified gravity direction vector.
20 May 2013
 Fixed an error reporting bug that was incorrectly reporting errors about linearized perturbation field pointers not being associated when linearized perturbation fields were not being evolved in a 2D simulation.
7 May 2013
 Added spatial derivative routines capable of acting on the entire elementbyelement flow field vectors, "dd?_all_elem()". *** Note these did not give a detectable speedup and so have been replaced with the previous code. ***
6 May 2013
 Added a command "AVG_ONE_DIR", which takes a flow field and averages it along a specified direction.
5 May 2013
 Modified the "MASK" command so that it accepts separate mask functions for each velocity field component, as well as the scalar field variable.
29 April 2013
 Fixed a bug in the new MPI_BCAST arrangement where an unallocated pointer storage was being accessed. The MPI_BCAST code segment should only be used when linearized perturbation fields are being evolved in addition to a base flow, but this was being executed when no linearized fields were being evolved accidentally.
8 April 2013
 Changed the MPI_BCAST operation during evolution of linearized perturbation fields  now all velocity components are gathered together into a single broadcast of the base flow field to all processors, rather than calling multiple MPI_BCAST calls each time step.
3 April 2013
 Completed the implementation of constant and linear forcing terms. This is implemented in 2D, 3D hexahedral and 2D linearized perturbation field solvers, including scalar fields. It is not yet implemented in SEFourier 3D simulations.
1 April 2013
 Began the implementation of "gvar_forcing_fX" and "gvar_forcing_gX" commands used to implement constant and linear forcing terms in each velocity component momentum equation or the scalar advectiondiffusion equation. Here "X" may be "u", "v", "w" or scalar field "s".
 *** Have disabled commands that overlap with this forcing facility, including "pgrad", "quasi2d", "gvar_lorenz" *** Still need to implement in time integration > constant forcing: SEFourier 3D, > linear forcing: SEF 3D
11 December 2012
 Ported Viper to MPI for parallelization. The MPI version, Vmpir, is not publicly available, but may be made available on request. Parallelization is employed both for SEFourier 3D simulations (tested beyond 1000 CPUs) and linear stability analysis evolving multiple linearized perturbation fields in parallel.
24 October 2012
 Fixed a bug in the advection for scalar fields in Cartesian coordinates relating to an incorrect spatial derivative term.
17 October 2012
 Fixed an incorrect vector size allocation issue for the Arnoldi solver for linear stability analysis and transient growth analysis that had emerged since the recent alteration to the numbering of velocity components. This bug only presented itself when different boundary conditions were imposed on each of the velocity components.
20 September 2012
 Fixed a divergence error in the evolution of linearized scalar fields for linear stability analysis in cylindrical coordinates  I had forgotten to premultiply the intermediate velocity field after the advection substep by the radial coordinate.
 A segmentation fault in linearized evolution of perturbation fields in cylindrical coordinates with nonzero swirl (wvelocity) was traced to a requirement that WVEL be placed before PERT.
17 July 2012
 Updated the code used to advect scalar perturbation fields. Previously the advection routines were incorrect/unimplemented for linear stability analysis with scalar fields active.
 Moved routines around in source code package. Previously most linear perturbation field routines and scalar field routines were in "floquet.f90" and "scalar.f90". Now advection, pressure and diffusion routines for base, perturbation and scalar fields are grouped in "ti_advect.f90", "ti_pressure.f90" and "ti_difusion.f90". This is intended to improve the logic and readability of the source code package.
17 June 2012
 Initiating a MAJOR change behind the scenes in Viper's code base. This will facilitate the ability to SEPARATELY specify Dirichlet OR Neumann boundary conditions for the individual u, v, & w components of velocity. This requires new mapping arrays to be built for each component, and separate Helmholtz matrices to be constructed for the diffusion solves for each component. New "btag" variables "u", "v" and "w" can now be specified instead of the combined "vel" specifier, though the legacy "vel" specifier will be retained for backwards compatibility. Users should compare results using new versions of Viper with their previous versions to ensure that nothing has been broken by this upgrade.
18 May 2012
 Fixed a bug in the 'int' command where an unassociated scalar field array was being addressed in spectralelement/Fourier 3D computations which did not feature a scalar field.
11 May 2012
 Found another bug in the 2D Boussinesq solver that crept in when the "iterate" command was implemented  this one was seeing the buoyancy correction use potentially the incorrect temperature field vector for the buoyancy term.
 Found a bug whereby a mask used to identify nodes lying on the axis (with global velocity numbering) could be corrupted in some scalar field simulations conducted in cylindrical coordinates, depending on the boundary definitions. This has now been fixed.
10 May 2012
 Fixed a possible error that emerged when "bouss" was called in an axisymmetric simulation, where a stressfree boundary intersected the symmetry axis. In this situation, it was possible that the inexact stressfree boundary correction could lead to a nonzero radial velocity on the corner node of an element sharing the symmetry axis and stressfree boundaries. This was leading to contamination of the velocity field, first across the symmetry axis, and then throughout the flow. Now, at the end of the symmetry boundary correction, the code explicitly sets the vvelocity to zero along the symmetry axis.
2 May 2012
 Changed the behaviour of "gvar_dt", "gvar_RKV" and "gvar_scalar_diff" in the "viper.cfg" file  they now accept functions of userdefined variables, so "gvar_usrvar" functions/variables may be created prior to defining those parameters, and the parameters can then be expressed in terms of those user defined quantities.
 Changed the functionality of the "bouss" command to permit separate specification of the sign of the required buoyancy (the first coefficient) and the magnitude of the buoyancy coefficient. The command now also accepts the buoyancy coefficient as a function (including userspecified variables).
17 April 2012
 Added a DEBUG command, which can be called to write more info to screen during program execution for bughunting purposes. Note that additional output will not be large, and will change over time as new features are added into the code.
 Added the capacity for iteration of 2D base flow fields to see if rates of temporal convergence could be improved for Boussinesq simulations. This is invoked using the new command ITERATE.
14 December 2011
 Added a new command "moments" which calculates the moment on a boundary in 2D simulations.
27 October 2011
 Further minor code tweaks to several routines to correct warning messages from "gfortran" compiler.
26 October 2011
 In a behindthe scenes alteration, some source code has been restructured to replace Intel Fortranspecific syntax with standard Fortran code, making it easeir to compile Viper on nonFortran systems.
13 July 2011
 Fixed the bug in the "TRACK TECP" binary output routine (I hope!).
 A bug in the "TRACK TECP" binary output routine was detected. Currently this feature doesn't work  use "ascii" option.
9 July 2011
 Added a new function, "track load" to load binary restart files for particle transport simulations that are created using the new "track save" function.
 Changed the operation of "track save". It now outputs a binary restart file for particle transport.
8 July 2011
 Added a new function "track tecp" to output particle data in binary or ASCII Tecplot point data file formats.
 Changed the "track save" functionality (ASCII output of particle locations) to "track saveold"  this is to facilitate the upcoming particle transport restart (save/load) functionality.
6 June 2011
 Added the gradients of the scalar field to the "vars ddx" output option of the "tecp" command. This is implemented for all 2D solutions and 3D hexahedral solutions.
11 April 2011
 Fixed a bug in the math parser routine, where the function simplifier code was returning functions with two negative signs next to each other when it pre evaluated an addition/subtraction of two constants appearing after a negative sign, that returned a negative number as the result of the evaluation. e.g. 'x  3  5' was becoming 'x2' rather than realizing that a "3" is being added to a "5".
4 April 2011
 Fixed a bug in the math parser routine, where the function "rand" was being incorrectly preevaluated by the function simplification code.
13 January 2011
 Found a bug that slightly corrupts the SEFourier 3D solver when a scalar field and Boussinesq computation are being carried out. This related to incorrectly mapping the differently numbered velocity and scalar fields during advection calculations.
5 January 2011
 Modified LOAD behaviour: If called before INIT, the "viper.cfg" Dirichlet boundary conditions are applied, overwriting whatever was on those boundaries in the restart file. If, however, the user calls LOAD after INIT, then static Dirichlet boundaries remain as whatever they were in the file, rather than what is in the "viper.cfg" file.
20 December 2010
 Added scalar field to SEFourier 3D versions of "int" and "flux" routines.
16 December 2010
 Added an option "hugh" to the "save" command to output 2D or perturbation field data to an ASCII file readable by SEMTEX, Hugh Blackburn's spectral element code.
13 December 2010
 Added Boussinesq modeling capability to SEFourier 3D solver.
11 December 2010
 Problem solved: turned out to be a boundary condition mismatch.
 Tested the "AXIROTATE" command in both the base flow and linear stability analysis solver  something not working in perturbation field solver  investigating.
10 December 2010
 Added the "AXIROTATE" functionality to the cylindricalcoordinate formulation of the spectral elementFourier 3D solver.
 Added a command "AXIROTATE", which allows 2D computations conducted in the cylindrical coordinate system (with "AXI" activated) to be computed in a rotating reference frame.
30 October 2010
 Added a command "INTF", which works like "INT", except it outputs integrals evaluated on each mode of a spectral elementFourier 3D computation, and the results are output similar to the "ENERGYF" command. Note that the integrand is calculated using the magnitude of the complex Fourier coefficients at each mode. e.g. This command could be used to generate output akin to ENERGYF, but with a conditional mask so that the energy is only integrated over part of the domain (say, radial values less than 2 units in cylindrical coordinates). The usage in this case would be: \> intf u 'if(y<2,u^2+v^2+w^2,0)'
 Added an option "nozero" to the Tecplot output routine "tecp", which removes the fundamental mode (zerowavenumber component) from spectral elementFourier 3D data when creating a Tecplot output file. This is especially useful for isolating 3D/nonaxisymmetric flow structures from swirling flows in cylindrical coordinates, for example.
16 September 2010
 Modified behaviour of LOAD routine so that static Dirichlet BCs are overwritten by what is in restart file rather than what is specified in "viper.cfg", provided LOAD is called AFTER "init".
2 September 2010
 Windows builds are now available for both IA32 and Intel 64bit architectures.
28 July 2010
 Added the capability for the SVD driver to optimize energy norms containing only a selection of the available velocity component / scalar field, rather than all of them (i.e. norm E = u^2+v^2+w^2, plus s^2 if active), which is the default behaviour.
15 July 2010
 Added the capability for GETMINMAX and SAMPLE to be employed on perturbation fields.
7 May 2010
 Altered the Tecplot output of SEFourier 3D data in cylindrical coordinates so that if the data set represented a complete revolution in the azimuthal direction, the data set would be continuous through the interface between 2*pi and 0 degrees (this eliminates the appearance of a boundary plane visible in the data).
5 May 2010
 Added the ability for the default SAVE and LOAD routines to interpolate onto different meshes (using the m option, just as with the old SAVE/LOAD functionality). To use, first the "m" option must be added when saving a solution, and then to interpolate onto a different mesh, include the "m" option when calling LOAD in a new simulation.
1 May 2010
 Added diffusion to 2D and 3D hexahedral particle transport algorithms, by way of a Gaussiandistributed random walk. The user specifies the Schmidt number using the new "TRACK DIFF" command.
1 April 2010
 Parallelized calculation of "track sample" output from SEFourier 3D particle tracking computations, and now carrying out file output of particle data only once per time step to hopefully improve the speed of this process.
29 March 2010
 Fixed an error in a formatting statement for output from "TRACK SAMPLE" in SEFourier 3D simulations.
 Added the scalar field to MONITOR output from SEFourier 3D computations.
21 March 2010
 An unallocated scalar numbering index was being incorrectly referenced in simulations in which the scalar field was inactive, leading to a segmentation fault. This has been fixed.
20 March 2010
 Added scalar field to Tecplot output from SEFourier 3D computations.
 Added treatment of scalar fields from SEFourier 3D jobs to save & load routines.
19 March 2010
 The option of adding a scalar field to the SEFourier 3D algorithm has been added to the code's functionality.
 The particle tracking "track save" routine has been rewritten for SEFourier 3D computations  now instead of calculating data, opening, writing and closing the output file for every particle (which can be many), the code now only opens and closes the file once  all calculation and output is conducted while the file is open. Have also removed the redundant transformations from parametric to physical space and back again. These alterations should make the "track save" command faster.
5 March 2010
 Added scalar field its spatial derivatives to output from 2D/3D calls to "SAMPLE".
6 February 2010
 Changed the norm being optimized using the SVD transient growth solver to "u^2+v^2+w^2+s^2" when a scalar field is active. Previously, a standard kinetic energy norm of "u^2+v^2+w^2" was implemented.
12 January 2010
 Altered behaviour of RAND for SEFourier 3D computations. By default, less randomization is applied to higher wavenumbers. However, this behaviour was also being applied even when the user was discretely randomizing a single mode (using the "k" option in RAND). This behaviour has been corrected  now when the "k" option is used, the specified level of randomization will be applied to whichever mode is specified.
6 January 2010
 Altered output from "energyf" command so that the fundamental mode was listed as mode 0, and corrected the behaviour of "rand", so that spanwise/azimuthal Fourier modes were randomized from number 0, not 1, when using the "k" option.
21 December 2009
 Improved the streamfunction output in Tecplot: now a Stokes streamfunction field is plotted from axisymmetric computations (thus the streamfunction contours align with streamlines for computations in both Cartesian and cylindrical coordinates.
18 December 2009
 Found an incorrect derivative term in the advection step of the linearized NavierStokes solver in Cartesian coordinates with an outofplane (w) component of velocity active in the base flow.
15 December 2009
 Fixed an error with particle tracking in cylindrical coordinates, in which an interpolation error in the azimuthal component of velocity caused particle trajectories to oscillate within elements touching the axis of symmetry.
1 December 2009
 Added the "psi" (streamfunction) variable to the Tecplot output routine. This allows a streamfunction field to be calculated and output in a 2D simulation without an outofplane velocity component.
27 November 2009
 Removed diagnostic output from SVD computation as testing has shown it to be working. This speeds up the algorithm.
 Corrected a pointer initialization affecting initialization of simulations of a scalar perturbation field on some platforms.
25 November 2009
 Modified the Arnoldi solver so that if a user initiates an Arnoldi solve from a restart file from a previously converged run, then the solver will recompute and ouput eigenvalues/eigenvectors from that run (this capability allows Tecplot or Viper restart files to be regenerated from an Arnoldi restart file, or when using the SVD solver, it will allow multiple transient growth calculations from a single set of computed eigenmodes.
 Fixed an eigenvalue normalization issue affecting the SVD driver routine  this now appears to work  for cases with and without a scalar field, this routine seeks to find the transient growth based on maximisation of a kinetic energy norm.
15 November 2009
 Two significant capabilities have been incorporated into the LOAD routine: (1) A new option "r" has been added which allows a user to add the stored flow field in a restart file to the present computation, rather than the default behaviour of overwriting the existing fields. (2) A new option "a" can be used to specify a scaling factor for the field being loaded.
8 November 2009
 Replaced two MAT_MUL functions in the static condensation evaluation routines with BLAS "dgemv" routines  this has resulted in a small speedup (approx. 5%) in some cases.
6 November 2009
 Now normalising ouput flow field from SVD transient growth routine to unit energy.
 Streamlined the time step timing statistics output further.
 Improved output during Arnoldi iterations  the code now provides more information about how many cycles remain in each Arnoldi iteration.
4 November 2009
 Fixed an error in the output of timing stats for sparse solves in parallel computations.
 Modified the QUASI2D command to allow the friction coefficient to be supplied as a function.
3 November 2009
 Altered appearance of "max du" monitors during time stepping.
2 November 2009
 Added file format checking to Arnoldi restart files. Viper will now report an error and abort an Arnoldi call if an existing restart file is incompatible (this follows the 25 October 2009 update to the Arnoldi solver).
31 October 2009
 Improved spatial derivative timing test so that more accurate timings are calculated when parallel jobs are executed.
30 October 2009
 Fixed the perturbation field time step routine so that zerowavenumber fields can be computed in axisymmetric computations with swirl.
27 October 2009
 Tweaked text output in LSA and TG driver routines.
25 October 2009
 Changed the format of vectors passed to the Arnoldi solver. *** Note that this modification makes new Arnoldi restart files incompatible with the previous format. *** Only the homogeneous node values are passed to the solver, rather than the full velocity vectors  the zero Dirichlet nodes have been omitted to reduce the size of the Arnoldi matrix problem. Also, the real and imaginary components of the perturbation field vectors are all sent to Arnoldi as entries of a real vector (as a complex vector sent to Arnoldi implies a temporal representation, whereas complex perturbation field vectors imply a spanwise spatial representation).
20 October 2009
 Corrected a flaw in the LOAD routine. It was incorrectly giving an error if a user tried to map a field onto a Fourier mode in an SE/Fourier 3D simulation.
15 October 2009
 Removed a superfluous blank line written to screen at the conclusion of the SAVE command. SAVE output should now more closely resemble TECP output.
12 October 2009
 Chris Butler detected an aliasing artefact in vorticity which arose in postprocessing for Tecplot output from SE/Fourier 3D computations. This was traced to incorrect indexing of conjugate (ve frequency) modes when vorticity fields were reconstructed from Fourier modes.
11 October 2009
 Added driver commands for both linear stability (Floquet) analysis and transient growth analysis. These are named "LSA" and "TG", respectively. These single commands can be used in place of a loop containing ARNOLDI and STEP commands.
8 October 2009
 An error with one advection term in the adjoint solver for transient growth analysis was uncovered after testing code against Hugh Blackburn's SemTeX code. This has now been fixed.
1 October 2009
 Added additional error messages to the PARDISO sparse matrix solver error reporting routine to conform to the Intel MKL 10.xxx implementations.
30 September 2009
 Added error messages to the LOAD command: If a user specifies a nonzero perturbation field number ("k" option) before FLOQ has been called, an error message is reported and the LOAD routine terminates. If a larger perturbation field number is supplied than what is active, a warning is printed to screen, and the fields are loaded onto the largest active field number.
 Joine identified a problem with stability analysis of frozen base flows with the ROTATE command active. It was identified that the rotating component of the base flow was not being subtracted in spatial derivatives of the base flow (this arose from a change to the code on 8 September).
27 September 2009
 Slightly modified the names given to output from the Arnoldi command to more closely resemble default output from SAVE and TECP, and to remove association with Floquet analysis, as ARNOLDI is also used for transient growth.
26 September 2009
 Added timing output during time stepping in transient growth calculations.
 Added code to the time integration routine for perturbation fields which is designed to minimise the error in the first couple of time steps due to the unavailability of the full number of velocity fields from previous time steps. This is achieved by reducing the order of extrapolations from previous flow fields so that only available fields are used. i.e., if the time integration order is 3rdorder, then for the first step a 1storder scheme is used, and for the second step a secondorder scheme is used. Thereafter, the 3rdorder scheme is used as normal.
23 September 2009
 The "max du" monitor in the SE/Fourier code was being incorrectly calculated thanks to misplaced brackets. This has been rectified.
22 September 2009
 Fixed a problem with 3D hexahedral particle tracking  some particles were occasionally misplaced during element interface traverses.
 In the particle tracking routien, added a CLOSE statement to release the injector coordinate file "track_pts" immediately after the points have been read.
 Changed the misleading reference to "Element order N = ..." during startup to "Element polynomial degree N = ...". The element order is technically N1.
 Repaired a recently arisen pointer access violation causing the 3D hexahedral solver to crash during the INIT call. No time integration algorithms were affected.
18 September 2009
 Added a runtime check of the speed of three alternatives for computing spatial derivatives: sparse matrixvector product, full matrixvector product, and tensorproduct construction (summed from derivatives in each elemental parametric coordinate direction). The fastest method is determined at runtime, and is used for that simulation. The sparse matrixvector product approach was used previously, but the other methods can be significantly faster, depending on the element polynomial degree used, etc. Users may notice decent speedups, particularly in the advection term, and to a lesser extent in the pressure term, which are spatialderivativeevaluations intensive.
9 September 2009
 Followed up on an issue reported by Chris Butler  the Intel Fortran compilers (ver 11) issue an internal compiler error when they encounter an "OMP ATOMIC" statement in a subroutine which may or may not be in a parallel region of code. Temporarily removing these calls (they are not solutionsensitive as they are only used for timing statistics).
 Attempting an explicit initialization of the elemental Laplacian matrices to zero at the beginning of its construction to see if this avoids the NaNs appearing in the matrix when it is build using the Intel MKL 10.1 libraries on the EM64 xe.nci.org.au machine.
8 September 2009
 Stability analysis of multiple wavenumbers now skip fields which have converged, permitting runs to speed up as modes converge. Note that to exploit this speedup, multiple modes will have to be computed on each parallel thread.
 Shifted computation of base flow spatial gradients from within base flow and perturbation field advection calculation routines to a single evaluation at the beginning of each time step. This will permit faster time integration in stability analysis computations. Speedup tests have shown improvements of approximately 10%.
7 September 2009
 Added an extra optional parameter to the ARNOLDI command which allows the user to specify the exponent of the convergence tolerance (the default is 7, i.e. 1e7). Numbers with smaller magnitude (e.g. 5, 4, 2, etc.) will converge faster, but with less precision.
6 September 2009
 Fixed an issue with the solver for zerowavenumber perturbation fields in cylindrical coordinates, which should allow the zeroth wavenumber (axisymmetric) mode to be evolved stably.
5 September 2009
 Added checking for infinities as well as NaN values when monitoring for divergence of the solution due to instability in time integration.
2 September 2009
 The Womersley profile feature has been implemented in both the 2D and SE/Fourier algorithms (AXI should be used in both cases).
1 September 2009
 Added a feature to WOMERSLEY command which allows the profile to be constructed given areaaveraged velocity Fourier coefficients instead of pressure gradient coefficients.
 Added a command WOMERSLEY which can be used to impose a Womersley velocity profile based on supplied Fourier coefficients of an imposed timedependent pressure gradient.
 Added an error report which picks up when a usersupplied function has imbalanced brackets.
24 August 2009
 When loading from stored flow fields and changing the time step, this causes the fields at the two stored previous times to be incorrect. A routine has been added to the "SET DT" command to interpolate these previous fields to the correct times. Theoretically, this should provide for a smoother restart for runs where the time step is abruptly changed.
22 August 2009
 Internal code change: On IA64 systems (such as APAC), the 2D time integration routine was not being optimized during compilation due to its large size. This has now been restructured, with parts of the routine being split off into separate subroutines, and compilation is again optimized. Users may notice a speedup for 2D simulations on APAC.
 Internal code change: Previously was using a condition "Nfloq_modes>0" to check if perturbation fields were active. This has been replaced by a LOGICAL condition "isFloq", which is quicker to evaluate and cleaner in the source code.
 Time integration now includes trapping of divergence due to time integration instability. Viper will immediately terminate if divergence is detected. Divergence is considered to have occurred if "NaN" is present in any of the flow field vectors.
16 August 2009
 Added the capability for RECONSTORE to store spanwiseaveraged snapshots from SE/Fourier 3D computations (previously only storage of 2D quadrilateral or 3D hexahedral solution data was facilitated). This will enable stability analysis of the 2D field corresponding to a spanwise average of a spanwiseperiodic 3D solution to be carried out (note that the spanwiseaverage is simply the fundamental mode of the spanwise Fourier modes of the flow field).
14 August 2009
 Added evaluation of the present time base field when starting transient growth calculations from a reconstructed base flow. This may have been leading to a slightly inaccurate startup of the perturbation field.
11 August 2009
 Added evaluation of the base flow velocity field at two previous times at the beginning of transient growth computations in which the base flow is being reconstructed (to ensure that the guess of the futuretime base flow used in advection of the perturbation fields is correct). This is also done to acquire the two futuretime fields for the adjoint (backwardstime) solve in transient growth analysis.
 Further speedups to RECONLOAD: The pressure field is now not computed by default, as it is not needed for stability analysis. Users can choose to compute this using the "p" option in RECONLOAD. Speedups for Fourier, polynomial and Akimabased reconstructions were 21%, 28% and 27%, respectively, for the test case used for the 10 August 2009 tests.
10 August 2009
 Added quasi2D friction term to perturbation fields for linear stability analysis: as with the base flow, this term is activated if the command QUASI2D is called. It should only be used in 2D Cartesian simulations, and the stability analysis will only be meaningful if the zeroth (2D) spanwise wavenumber is analysed (i.e. calling FLOQ 0).
 Added BLAS "daxpy" calls for vectorscalar product & sum operations in polynomial interpolation routine using RECONLOAD. In a test case, a speedup from 4.8 to 4.0 seconds was achieved (17% improvement).
 Added BLAS operations to substantially speed up the Fourier reconstruction code using RECONLOAD  a test case improved execution time from 31.9 seconds to 2.15 seconds (93% improvement).
 An error was found in the advection term for perturbation fields in cylindrical coordinates with a nonzero wvelocity in the base flow (e.g. Stuart Cogan's swirling flow in a cylindrical vessel driven by a rotating base). This was introduced with the June/July transient growth upgrade, and has now been fixed. Tests with a 4 May 2009 build of the code verified (a) that the Cartesian and regular axisymmetric solver was unaffected (both for base flows and perturbation fields), and (b) that there was no difference between predicted stability modes using either STAB or ARNOLDI.
9 August 2009
 Found a LAPACK work array which was not being properly deallocated before a new allocation call  this affected ARNOLDI calls with WVEL active only.
 Achieved substantial speedups to polynomial and Akima reconstruction of velocity fields using RECONLOAD. A test case had 25 elements of order 7, and 20k steps took 50 seconds to time integrate. Polynomial interpolation sped up from 77.8 to 4.70 seconds, and Akima interpolation sped up from 36.7 to 4.61 seconds.
8 August 2009
 Added Akima spline interpolation to the flow field reconstruction schemes available through RECONLOAD  this should be used in place of polynomial interpolation to avoid wiggle artifacts.
6 August 2009
 Fixed an array reference bug in RECONLOAD  the pressure field was being reconstructed incorrectly, which sometimes led to unexpected crashes due to an incorrect array being referenced.
5 August 2009
 Fixed a bug in Cartesian scalar transport  the 30 July fix to the cylindrical scalar transport algorithm broke the Cartesian solver. Now both are fixed.
 Adjusted output of TRACK SAMPLE commands so that the report to the user that data is being written is only performed once per call, not once for each particle.
3 August 2009
 Added to the output of the FLUX command the integrated absolute value of the flux in addition to the integrated flux along each boundary.
2 August 2009
 Improved ARNOLDI command  now works for multiple perturbation fields.
31 July 2009
 Repaired scalar field timer for 3D hexahedral simulations (it was not being calculated), and restructured timers so that the time to evaluate the scalar field would be counted for both reconstructed as well as integrated flows.
30 July 2009
 Added the potential to use polynomial interpolation as an alternative to Fourier interpolation when using RECONSTORE and RECONLOAD (formerly FOURIERSTORE and FOURIERLOAD).
 The previously unused cylindrical formulation of scalar advectiondiffusion contained an error which has now been rectified.
 Identified that the code used to eliminate duplicate interpolated nodes in the Tecplot output routine was particularly slow in larger 3D hexahedral jobs  attempts are underway to speed this up.
 Corrected output of "max ds" for scalar field computations in 3D hexahedral simulations  it now adheres to the minimum 1second between outputs.
22 July 2009
 Altered how output is displayed during time stepping. Now the step/time/max_du information is written at most once per second (thus not every time step is displayed). This will produce smaller output files, and will reduce the amount of I/O performed by the code  this should be helpful on queued systems such as SunGrid, APAC, etc.
 Changed the computation of the change in velocity and scalar field monitors during time stepping. Now, instead of outputting the maximum change in the magnitude of velocity in the flow, the monitor now displays the maximum change in any one of the velocity components (u, v, w). This should very slightly speed this calculation up, and will also avoid over or underflow of the calculated quantity.
21 July 2009
 An error in the Tecplot routine for SEFourier simulations was incorrectly rendering the pressure field  an alternative approach has been taken for reconstruction from the Fourier modes.
 An incorrect divisor was used in a MOD() function in the cylindrical formulation of the SEFourier code. This was causing incorrect radius values to be employed when computing the advection terms.
 To aid in bug diagnosis, all versions of the code are currently compiled with the "traceback" option set, which will report routine/line number information where an error might have occurred in the source code. Let Greg know if this noticably slows things down.
 Joine detected an error occurring during LOAD which was tripped by interpolation to a different element polynomial order under the optimized "em64" build of the code. This has been fixed by preallocating the pointer array causing the problem.
20 July 2009
 Fixed bug where a logical vector "radial_mask" was causing the code to crash in some 2D runs on Windows due to it not being allocated.
 ARNOLDI now uses the new restart file format (consistent with the standard SAVE and LOAD calls). Also, the extracted eigenvector fields are copied onto the fields at previous times to avoid an unexpected kick being added to a simulation upon restarting.
16 July 2009
 Changed the functionality of L2 and INT commands. They now take their parameters as options, and have the added capability of being used to evaluate integral norms on perturbation fields was well as base flow fields. Note also that the L_2 norm is now computed as the square of the velocity magnitude, i.e. u^2 = u^2 + v^2 + w^2.
11 July 2009
 Fixed a problem with the routine used to return the parametric coordinates within an element corresponding to a physical coordinate. This issue meant that on some occasions during particle tracking or returning values at a point (using SAMPLE), the coordinates were returned before the solution had fully converged. Note that the worstcase scenario would be that the returned point could be between the desired point and the nearest mesh node.
8 July 2009
 Tweaked the allocation of storage for the Fourier reconstruction of velocity fields. The facility now will use less memory, and may be a tiny bit quicker.
 Fixed an error in the routine used to compute strain rate at a point in the flow using either "SAMPLE" or "TRACK SAMPLE" in SEFourier computations in cylindrical coordinates  the returned strain rate may have been miscalculated at the axis.
1 July 2009
 Applied an indexing refinement to the FOURIERLOAD routine, which achieved a slight speedup.
30 June 2009
 Implemented new commands (FOURIERSTORE and FOURIERLOAD) which can be used to save a timevarying solution of a two or threedimensional base flow as a Fourier time series, which can be used in a subsequent simulation to supply a periodic base flow for stability analysis computations, for instance.
 Deleted some redundant code used to supply the radial distance in the global velocity numbering sense for computations in cylindrical coordinates.
23 June 2009
 The 2nd term in the advection operator for the adjoint operation in transient growth computations was incorrect. This has been rectified.
19 June 2009
 Added the capability of a scalar perturbation field to be used with ARNOLDI for stability analysis.
 Added the ability for the SAVE and LOAD routines to recognise scalar perturbation fields.
 Scalar perturbation fields should now be plotted in Tecplot files using TECP.
 Cleaned out some unneccesary matrix storage for base and perturbation scalar fields.
18 June 2009
 Added a scalar perturbation field for stability analysis, plus a Boussinesq term in the perturbation field time stepping routine. This will facilitate stability analysis of Boussinesq flows, and is activated by activating a scalar field and calling BOUSS as well as FLOQ.
 Modified the function simplifier to do a better job of trimming zeroes off numbers.
15 June 2009
 Added the capability to compute linear stability analysis of perturbation fields with zero wavenumber.
12 June 2009
 Implemented a trial algorithm for transient growth calculations of frozen base flows.
 Fixed an allocationrelated crash affecting Cartesian stability analysis jobs.
10 June 2009
 Fixed a minor bug in which scalar transport fields were receiving an undesired perturbation upon restarting from a saved file.
5 June 2009
 Reintroduced a little more output in the LOAD routine  in particular, the routine now reports if any parameter values are changed while reading from the restart file, and also reports on mismatches in resolution requiring reinterpolation.
3 June 2009
 Fixed a bug in the function simplification routine  functions which when simplified occupied a longer string were being accidentally truncated.
 Fixed a bug in LOAD where complex perturbation fields from stability analysis were not being mapped onto nonzero Fourier modes of an SEFourier 3D simulation correctly.
2 June 2009
 Detected a 2nd error with the LOAD routine in terms of mapping 2D fields onto SEFourier modes  all fields were being zeroed out each call to LOAD.
1 June 2009
 Detected an error in the LOAD routine caused by today's earlier fix, which could cause parameters to be loaded incorrectly. This has now been corrected.
 Added some further minor streamlining to the function simplification routine. Routine should now be a little less memory hungry.
 Fixed a bug in the LOAD routine which will now allow multiple 2D field to be mapped onto SEFourier modes using the "k" option.
31 May 2009
 Implemented some memorysaving and timesaving changes to the function simplification routine.
29 May 2009
 Made the output of the SAVE and LOAD routines more compact.
28 May 2009
 Restructured the LOAD routine for faster performance. Storage size is determined by making a pass of the file, avoiding excessive reallocation and copy operations.
 An error was identified in the SEFourier algorithm  the symmetry boundary correction was accidentally placed outside the loop through the number of planes, causing the "k" index to the plane number to be undefined, resulting in an outofbounds error.
 Replaced some conditional evaluations (e.g. "IF (Ninner > 0) THEN...") with a preevaluated LOGICAL parameter (e.g. "IF (isNinner) THEN...").
27 May 2009
 Stripped out nonNewtonian code and compressible flow code, which were unused, untested, or not implemented.
25 May 2009
 Added a timer for the scalar field in time stepping.
22 May 2009
 Implemented a new command, ORDER, which can be used to change the order of the time integration scheme used for velocity or scalar fields. Allowable integration orders from 1 to 3 are facilitated.
21 May 2009
 Identified that a use of the FORTRAN function MINLOC to determine the element and node numbers of a global number (used in boundary condition, initial field, Boussinesq updates, etc.) was extremely inefficient. These have been replaced by preevaluated indexing vectors.
 Implemented a command "FLUX", which outputs the flux of a scalar field through each boundary.
20 May 2009
 Fixed a bug in which scalar fields were not being properly loaded from restart files.
 Slightly restructured the parallelization applied to the advection substep of the SEFourier 3D code. Small speedups (1% to 5%) were observed in testing.
 Replaced some "divide by radial distance" operations with multiplication by a preevaluated 1/R quantity for efficiency.
19 May 2009
 Fixed a rare allocation error affecting 2D simulations in which there were more curved boundary filaments than there were numbered boundaries. Storage is now fully dynamic in the 2D curvature construction algorithm.
14 May 2009
 Improved the code used to print numbers to screen in a neat readable format. Now trailing zeroes are removed for numbers in scientific notation as well as decimal notation.
 Improved the function simplification algorithm. Functions of constants (e.g. "cos(45.23)" are now also preevaluated (trigonometric, exponential, log, absolute value and square root functions are facilitated. This can make userdefined functions smaller, and may result in a slightly shorter evaluation (hence compute) time.
 Fixed a bug in "gvar_movref"  if 3 functions were not explicitly supplied, the function simplification routine broke down attempting to process a null string.
12 May 2009
 Modified the functionality of STOPCRIT  now the current loop iteration that is being executed when the STOPCRIT threshold is breached is completed, and no further loop iterations are entered.
 Added a safeguard to avoid the endless output encountered when a "viper.cfg" file contains an "mesh_file" statement pointing to a file which doesn't exist. Now Viper will only prompt for a filename 3 times before terminating.
 Streamlined the SEFourier 3D algorithm: Only the positive frequency spectrum is stored in the velocity field vectors, reducing the memory footprint of the code. Also, all forward & backward transforms now computed at the higher resolution employed for advection (either 50% extra modes or 1 extra mode, if the "alias" option is employed in FOURIER or not, respectively). Thus the Tecplot output will be interpolated onto a slightly higher number of planes by default. These improvements to the SEF algorithm mean it is no longer necessary to perform tweaks to maintain stability  the pressure correction is now again enforced on all nonzero Fourier modes, and it is not necessary to explicitly enforce the fundamental mode (or the highest mode when the number of planes is even) to be real only.
8 May 2009
 As part of the improvements to the SEFourier treatment of Fourier modes and aliasing (i.e. all BCs and advection evaluated at higher Fourier resolution, and only lowest Mfourier_planes/2 nonzero modes being retained in the computations), the imposition of symmetry boundary conditions has been reverted to Fourier space (this was switched to real space in a recent attempt to hunt down the problem affecting Stuart Cogan's simulations, which in fact was caused by the incorrect order of the diffusion part of the highorder pressure BC).
4 May 2009
 Chris Butler reported some curious phenomena when complicated boundary conditions were imposed on 3D spectralelement/Fourier simulations  namely, if the boundary condition energized all Fourier modes in the discretization, then the solution was displaying degraded convergence, and the highest pressure mode was sometimes being contaminated with spurious artifacts. These effects were dependent on the relationship between the symmetry of the boundary conditions and the number of planes being selected. These problems were partially alleviated by switching off the pressure update of the highest Fourier mode in the computation (effectively relaxing the discretization of the pressure term). The contamination in the fundamental mode near the axis when an odd number of planes was used was found to be due to an insufficient Fourier resolution  care must be taken to ensure that sufficient Fourier planes are employed to fully resolve the boundary conditions as well as the interior flow.
1 May 2009
 Joine So reported a bug in GETMINMAX which was causing the code to abruptly crash. This was traced to a "log10" function that was sometimes being supplied a negative argument. This has been fixed.
22 April 2009
 Added an option "k" (which can replace the "n" option) to the FOURIER command. This allows the user to directly specify the number of Fourier modes they require in an SE/Fourier 3D computation  this will make parallel computations easier to set up, as the number of modes must be an integer multiple of the number of available threads for maximum efficiency.
 Added an option "alias" to the FOURIER command, which can be used to specify if the TwoThirds Rule is employed to avoid aliasing of the advection term. Without this option, advection is computed at the specified resolution, which can lead to aliasing on the highfrequency end of the Fourier spectrum if insufficient modes are included in the computation.
21 April 2009
 Repaired an error in the advection term calculations for SE/Fourier 3D simulations. The velocity field (u,v,w) components were being incorrectly backwardstransformed into real space (derivative terms were unaffected).
7 April 2009
 Reduced to 2nd order the diffusion term of the highorder pressure boundary condition (this preserves 3rdorder accuracy in time for the velocity field, Karniadakis & Sherwin, 2005).
 Attempted implementation of a "twothirds rule" to eliminate aliasing in Fourier space due to quadratic terms in the advection substep (SE/Fourier computations).
6 April 2009
 Modified RAND so that higher modes have a smaller level of random noise applied. This is intended to reduce the effect of instability caused by a cascade of whitenoise energy into the highest Fourier modes.
5 April 2009
 Altered the way symmetry boundaries are enforced in SE/Fourier 3D computations. Now the velocity field is transformed into real space, the symmetry conditions are enforced on each plane, and then the fields are transformed back into Fourier space. It is hoped that this is more stable than the previous implementation, which was causing problems when the highorder pressure BC was employed.
2 April 2009
 Performing the L2 norm calculations with double precision real arithmetic instead of complex arithmetic  this is intended to attempt to overcome the precision issues detected in l2 norm convergence studies from SE/Fourier simulations.
 Removed the highorder Neumann pressure boundary condition on symmetry boundaries, which should be exactly zero anyway. This was done to eliminate a spurious fuzziness in solutions detected by Stuart Cogan.
1 April 2009
 Fixed a rare bug in TECP for SE/Fourier simulations. If the velocity components were being plotted in addition to a userspecified field, the ucomponent of velocity was being incorrectly transformed from Fourier to real space. This has been fixed.
 Slightly tweaked the "RAND" command for SE/Fourier computations. Complexconjugate copies of positive modes to their negativefrequency counterparts is now done for the whole velocity field, rather than just the randomized parts. This ensures that the conjugate symmetry of the Fourier transforms are preserved.
31 March 2009
 Added a command "PBC" which toggles the highorder pressure boundary condition on or off (default on).
 Added a userspecified function plotting capability to the Tecplot routine (this is activated using the "u" option).
30 March 2009
 Chris Butler reported that the L2 norm is still displaying singleprecision convergence in SE/Fourier computations, despite other integral quantities demonstrating doubleprecision convergence.
 Fixed a problem reported by Nick Boustead with the new LOAD/SAVE commands: the code was incorrectly determining the number of data fields written to restart files for perturbation fields when no wvelocity was active in the base flow.
27 March 2009
 Added the particle ID number to "TRACK SAVE" output from particle tracking in hexahedral computations.
 A compiler optimization error in the Windows version of Viper was overcome by adding a "/Qtrapuv" option to force initialization of saved variables.
 In order to improve the precision of the L2 norm integration in SE/Fourier computations, all instances of SQRT(), ABS(), SIGN(), MIN() and MAX() in the code which operate on doubleprecision REAL or doubleprecision COMPLEX data types have been replaced by their specific doubleprecision routine (e.g. DSQRT in place of SQRT, etc.).
20 March 2009
 The SAVE and LOAD routines have been revised. A new format is now in use. Users wishing to SAVE/LOAD using the previous format can do so by calling SAVE or LOAD with the new "old" option specified. In addition, if a user attempts to call LOAD for an old restart file without this option, the code returns an error and attempts to load the data using the old routine anyway. The new routines are more flexible, allowing new facilities such as saving and loading solutions with scalar fields active, safer loading of 2D solutions onto SE/Fourier solutions and vice versa (e.g., loading an SE/Fourier solution onto a 2D simulation loads in the zeroth (or fundamental) mode). In general it is cleverer, as each field is clearly identified in the file, so Viper has a better idea of what it can possibly do with each field.
15 March 2009
 The SE/Fourier axis filter for cylindrical computations has been disabled by default. Users can still activate a filter by setting the "f" option in the FOURIER command.
14 March 2009
 The axis filter in SE/Fourier computations was found to slightly shift the Kovasznay flow solution from its exact value. Suggest users do not use the filter (i.e., explicitly set the "f 0" option in a FOURIER call in AXI computations.
 Fixed the SE/Fourier singleprecision convergence issue: problem was caused by generic FORTRAN functions for extracting real & imaginary parts of complex numbers, and for assemlbing a complex number from two real parts. Some of these were not preserving the KIND=8 double precision (64bit) floating point precision. Have replaced all REAL() with DBLE(), CMPLX() with DCMPLX(), AIMAG() with DIMAG() and CONJG() with DCONJG().
 While searching for the convergence issue, replaced the EXTERNAL PARDISO calls with a "USE MKL_PARDISO" call to the MKL FORTRAN 90/95 interface for the routine. This picked up a couple of INTENT statements that could have caused problems.
 Removed the "skip matrix solve if RHS=0" conditional check in sparse solver routines. This should never occur anyway.
13 March 2009
 Encountered an optimisation or compiler bug under Windows  incorrect (or at least imprecise) results were computed when degree 2 or 3 elements were used.
 Refined the code for computing integrals from SE/Fourier results using the INT command.
 Added further warning messages to the code used to parse the Viper configuration file  more errors are now detected and reported to help users confirm that their boundary conditions, etc. are being processed correctly.
12 March 2009
 Added a warning message to detect if a "btag" command is assigned to a boundary which is not included in mesh file.
 Added a warning message when parsing "viper.cfg" files  the code now complains if an unrecognised variable name is supplied in a "btag" statement.
 Streamlined the code calculating maximum changes in velocity. for display each time step. New approach should be faster (a square root is now performed only on the largest "du", "ds", etc., not the whole list of velocities.
10 March 2009
 Altered the implementation of the Boussinesq approximation  it is now based on a principle of calculating the effect of changes from reference temperature & density.
9 March 2009
 Chris Butler found a rare bug in the "simplify_function_string" code, which appeared for functions which had a number of components precisely matching the size of the storage arrays for this information. This has been fixed.
7 March 2009
 A couple of variables in the parallel part of the SE/Fourier time step loop were not designated as "PRIVATE" (or separate copies per thread). This may partially explain the "max du" convergence issue.
 Corrected a problem with scalar diffusion  a rouge "Nscalar_steps" variable was not deleted since the transition to a 3rdorder backwards multistep integration scheme for the scalar field. This was causing the scalar field to diffuse incorrectly.
3 March 2009
 Increased the size of the output of a function insertion routine in "usr_vars.f90" to permit the handling of very large compounded mathematical expressions.
 Replaced statically allocated storage vectors in "simplify_function" routine in "math_parser.f90" module for more efficient handling of large mathematical expressions supplied by the user.
2 March 2009
 Increased the allowable function string size to approximately 10,000 characters. This permits users to compound many userdefined functions into larger expressions. Storage for these arrays is now dynamic: the hardcoded allowance for 500 boundary is no longer in place.
 Fixed a bug in the boundary assignment code: when processing "viper.cfg" files, if multiple "btag" statements duplicated a boundary definition, the code incorrectly defined the boundary  the code now robustly assigns boundaries according to the last of any duplicated boundary assignments.
28 February 2009
 Corrected an issue where the pressure field was being incorrectly plotted in Tecplot dumps from SE/Fourier computations (the negative Fourier modes were being corrupted  this was particularly noticeable if a LOAD operation was performed.
27 February 2009
 Fixed a problem with the SE/Fourier Tecplot output algorithm: the "cartesian" option was not functioning properly when computations were carried out in cylindrical coordinates.
25 February 2009
 Success! Tests confirm that the new version properly indexes the highorder pressure boundary conditions, and also permits time integration to proceed without the erroneous restriction in time steps.
 Stu identified a problem in which time integration stability had decreased since the upgrade to the Neumann BCs implementation. This was traced to an indexing problem which was scrambling the boundary nodes on which highorder pressure BCs were applied. This has been corrected, and tests are underway to see if this rectifies the time step restriction problem.
24 February 2009
 Detected an allocation error in the scalartransport routine  this has been fixed.
 Improved the "gvar_movref" facility  the change is now calculated as a 3rdorder time derivative of velocity at the appropriate time, and this acceleration term is imposed on the intermediate velocity field during the advection solve.
23 February 2009
 Changed the scalar field advection/diffusion transport algorithm to a 3rdorder backwards multistep method (consistent with velocity field).
 Fixed a bug in the 2D Tecplot binary data file generation routine which only arose if the "SR" variable was requested.
 Fixed a pointer association error with a 3D pressure boundary condition vector.
21 February 2009
 Fixed an indexing error in the cylindrical SE/Fourier Tecplot routine, which was causing the strain rate and lambda_2 fields to be corrupted (particularly near the axis.
20 February 2009
 Strain rate magnitude and lambda_2 fields are still strangely speckledy, despite the velocity, vorticity and velocity gradient fields being smooth. This is only present in SE/Fourier plots generated from data in cylindrical coordinates... currently searching for the problem.
 Identified that the crashing in SE/Fourier computations in cylindrical coordinates was caused by a dividebyzero in a component of the highorder pressure gradient condition imposed on Neumann boundaries  initial tests indicate that this correction has fixed the SE/Fourier cylindricalcoordinates crash recently plaguing the code.
18 February 2009
 Identified that the thetaderivativeterms of the rateofstrain tensor used for strain rate calculations was being incorrectly computed. This is being rectified currently.
 Still searching for possible problem with SE/Fourier cylindricalcoordinates shear rate issue  lots of scrambled noise superimposed onto data for some reason.
 Tightened file creation in Tecplot routine: existing files are tested to see if they can be overwritten immediately before data written to file  thus there is less chance that this status will change between the file check and the new file creation.
16 February 2009
 Investigating a possible problem with the "SR" field in SE/Fourier computations.
 Added a command "gvar_init_scalar_field", which can be used to set an initial condition for a scalar field (similar to "gvar_init_field" for velocity and pressure).
 Fixed a bug with STOPCRIT in SE/Fourier computations  it was not stopping the time integration loop.
 Fixed some typos in the HELP facility (gvar_init_field had nothing to do with setting the nonNewtonian viscosity function).
15 February 2009
 Neumann BCs now activated for velocity components, pressure, & scalar fields in 2D and 3D computations. This is not yet active in SE/Fourier computations.
13 February 2009
 Began the rollout of Neumann boundary conditions for velocity, pressure and the scalar field. Neumann boundary conditions impose the outward normal gradient of a quantity (e.g. dp/dn) at a boundary, rather than the value itself (as with Dirchlet boundary conditions). These are selected using boundary type #2 (formerly used for timeindependent Dirichlet conditions). Neumann boundaries can be specified as a mathematical function, just as for Dirichlet (type #1) boundaries. Currently this is implemented for velocity only (requiring 2 or 3 components for normal gradients of u,v,w fields). Pressure and scalar field Neumann boundary capabilities are coming. This feature will be useful for imposing flux conditions in scalar field computations  particularly Boussinesq approximations of densitygradientdriven flows.
12 February 2009
 Identified and corrected a couple of errors relating to the calculation of velocity gradients and shear rate fields in the SE/Fourier version of TECP with the "vars" option active.
9 February 2009
 TECP "vars" option is now operational for 2D, 3D, and SE/Fourier simulations. Velocity gradient variable names are now correct for both cylindrical and Cartesian coordinates.
7 February 2009
 Velocity gradient variables added to 2D & 3D Tecplot files output from TECP using "vars" option.
 Removed the "o" option in 2D & 3D Tecplot output with TECP. This functionality has been replaced with the "vars" option. Still need to add velocity gradients in 2d/3d TECP routines.
6 February 2009
 Rolled out an option "vars" in the TECP command (currently SE/Fourier only), which allows users to select which of a number of different variables they wish to include in a Tecplot data file. The "sr" option has now been removed. This new option will eventually replace the seldomused "o" option.
5 February 2009
 Fixed a bug in the boundary curvature algorithm which was failing to correctly curve boundaries which comprised two 2D element edges.
15 January 2009
 A bug was detected with the Floquet solver when base flows were frozen, that had emerged during a recent alteration to the way that base flow velocity fields were distributed between threads. This has now been corrected. A vector copy was being skipped accidentally when "FREEZE" was active.
9 January 2009
 Added option "sr" to TECP command, which can be used to output the shear rate field in an SE/Fourier computation.
8 January 2009
 Added capability of outputting only a single Fourier mode in SE/Fourier computations with the new "m" option in the TECP command.
7 January 2009
 Added "cylindrical" and "cartesian" options to the TECP and TEC_FLOQ commands. For SE/Fourier 3D computations, and 3D reconstructions of Floquet modes, these options can be used to determine whether the velocity components in the Tecplot data files are the cylindrical or Cartesian components. The Cartesian form can be useful for plotting velocity vector fields in Tecplot (which natively uses a Cartesian coordinate system for plotting).
 Removed the particle concentration field from Tecplot plotting of SE/Fourier solutions as it was interminably slow to compute.
2 January 2009
 Parallelised the particle tracking time integration routine. Testing indicates this has been successful, and will produce a solid speedup for particle tracking simulations.
 Observed that the TECP output of SE/Fourier computations when particle tracking is active is interminably slow. Attempting some speedups and parallelisation to help, but may have to rely on TRACK SAVE for plotting.
 Fixed a couple of typos in the HELP documentation.
 Noticed that some OpenMP statements in the source code were incorrectly decorated as !OMP rather than !$OMP, and were therefore not activating. These were only intended to activate when singlefield computations were being performed, so this does not effect the Floquet analysis and SE/Fourier parallelizations.
26 December 2008
 Added time "t" to the variables that can be adjusted using the "set" command.
 Added the capability to include an "=" sign in "set" expressions, and have also now added the possibility of changing "RKV", "dt" and "t" variables directly on the Viper command line, e.g. "t = 50", "dt = 0.005", "rkv 100", etc.
23 December 2008
 Fixed a bug with the "gvar_lorenz" command, which was not correctly evaluating the function during time stepping.
15 December 2008
 Added underrelaxation and additional tolerancing to the routine used to determine the parametric coordinates within an element of a physical coordinate on the mesh. This is designed to resolve problem of some points very near to a boundary being erroneously identified as lying outside the domain.
14 December 2008
 Tests confirm that the parallelisation bug in the Floquet solver has been fixed.
13 December 2008
 Altered the way base flow velocity field data is passed across threads during parallel stability analysis  testing underway to determine if this fixes the aforementioned bug.
 Identified a problem with the Floquet solver when running in parallel. If multiple fields are being computed per thread (e.g. 8 fields over 4 threads), then the multipliers returned for the perturbation fields computed on the same thread as the base flow are nonsensical. This is being looked into currently.
11 December 2008
 Added the particle ID number to the output of "TRACK SAMPLE".
10 December 2008
 Recoded the "AUTOCORRF" routine  the returned autocorrelation is now normalised (that is, mean is subtracted from the input signals, and the resulting correlation is divided through by the variance). This means now that the returned correlation will always vary between 1 (correlated), through 0 (uncorrelated) to 1 (inversely correlated).
8 December 2008
 Added two new commands to be used to facilitate magnetohydrodynamical modeling: GVAR_LORENZ can be used in the "viper.cfg" file to express functions for the components of the Lorenz force, which is then calculated and added to the computation of a 2D fluid flow. QUASI2D can be used within Viper to establish a "quasi2D" treatment of a 3D flow by modeling it using the 2D solver with the addition of a friction term, subtracted to mimic the effect of friction on the walls of the domain in the outofplane direction.
5 December 2008
 The command "sample" has been improved  it should now work for SE/Fourier computations.
 Tightened scratch file creation/removal code in "loop" command  a file that was no longer written to since yesterday's update is no longer created.
4 December 2008
 Added a new function to the mathematics parser: "rand(x)". The function returns a random number between 0 and x, and is always treated as a time and spacevarying function. That is, a different random number is calculated at each node where the function is evaluated, and at each time that the function is evaluated.
 Further pointer association problems with particle tracking have been fixed. Now Viper is more robust to "track" calls appearing either before or after "init".
 Modified the "loop" code to create smaller scratch files storing looped commands. These are invisible on Windows systems, but are visible on Linux systems. Nevertheless, if a very large loop is constructed, a significant pause could result as the scratch file containing the unrolled loop commands was created. Now these files are much more compact, as the loop counting is handled by Viper, and the commands are not repeated in the scratch files.
 Added screen output during loops informing the user how many loop iterations remain.
3 December 2008
 Addressed a bug with "track seed" for SE/Fourier computations  the zdirection pointers were not being associated properly.
 Added particle concentration field to Tecplot output of SE/Fourier 3D computations.
2 December 2008
 Slightly modified "max du" output during time stepping: an extra blank has been added to clear any extra characters from the line when the next time step is performed.
28 November 2008
 In ARNOLDI the limitations on number of eigenvalues being sought and the number of Arnoldi vectors in a computation (formerly 12 & 30, respectively) have been removed. Storage is dynamic, so these limitations are redundant.
 Slightly tweaked the format of output of floatingpoint numbers to get rid of annoying blank spaces in output.
25 November 2008
 Have modified the Tecplot plotting routines and the particle tracking routines to rectify some issues with SEFourer particle tracking not following the flow properly.
21 November 2008
 Observed that particles were not following the flow precisely when flowing across an (rtheta) plane in cylindrical coordinates. This is being corrected.
20 November 2008
 Fixed some bugs relating to the plotting of SE/Fourier computations, as well as particle tracking in SE/Fourier computations  particularly in a cylindrical framework.
19 November 2008
 Implementation of the particle tracking algorithm in the SE/Fourier 3D flow solver has been completed. The particle concentration field in Tecplot output has not yet been implemented for these computations. The track_pts file should contain points in the coordinate system being computed in (i.e. xyz for Cartesian, zrtheta for cylindrical).
18 November 2008
 Beginning to roll out an extension of the particle tracking algorithm to the SE/Fourier code.
17 November 2008
 Implemented the Neumann highorder pressure boundary condition in computations of the perturbation fields in the stability analysis solver, and in spectralelement/Fourier 3D computations. These solutions should now be formally 3rdorder accurate in time.
13 November 2008
 The issue with ARNOLDI producing quasiperiodic modes may be a roundoff error issue. Test cases suggest the imaginary component is very small. Further investigation underway.
 For ARNOLDI, when base flows have 3 components (using WVEL), the complex matrix capabilities of the ARPACK package are now being used. It has been shown that this produces the same results as the previous treatment. This revision to the code has modified the Arnoldi restart file created by Viper (it is now substantially smaller)  users will not be able to restart from Arnoldi runs from earlier versions of the code.
12 November 2008
 Added extra output to ARNOLDI routine: the error/information messages reported by the ARPACK subroutines are now output to screen. This information is provided to the user as the ARNOLDI iterations can occasionally fail for some combinations of parameters, etc.
 Stuart Cogan reported that ARNOLDI is returning complex Floquet multipliers when the power method predicts real multipliers during cylindrical computations with swirl (nonzero azimuthal velocities). This is being investigated further.
13 October 2008
 Added a new functionality to RAND: now using the k option, randomisation can be applied to a single mode.
10 October 2008
 Fixed a small bug where some outofrange integers were incorrectly displaying as negative numbers during startup and initialisation for very large jobs. This did not effect the flow calculations.
7 October 2008
 A bug was detected in the SE/Fourier algorithm: When timedependent boundary conditions were used during a parallel simulation, the Discrete Fourier Transform package used to calculate the Fouriertransformed boundary conditions returned incorrect results. The Fourier transforms are now performed outside the parallel loop.
30 September 2008
 Fixed a bug in the Tecplot Floquet reconstruction routine  in cylindrical coordinates, the radial and theta (v & w) components of velocity are now output in cylindrical coordinates.
26 September 2008
 A bug in the Tecplot plotting routine has been addressed, where vorticity components in cylindrical coordinates were being incorrectly computed.
 The TEC_FLOQ routine has now been updated to include velocity components.
 Reintroduced a filter to the SE/Fourier modes  the highestfrequency modes are forced to be zero  the complete removal of the 2/3rds rule yesterday resulted in unstable computations.
 Added a command "AUTOCORRF", which outputs the spanwise autocorrelation of the velocity components at a sampled point on the spectralelement mesh.
 Found and rectified a bug in the "SAMPLEF" routine  an incorrect reference to v & w velocities was corrupting the results.
25 September 2008
 Removed the "2/3rd rule" filter applied to the outofplane Fourier modes in SE/Fourier computations. This was introduced during implementation of the scheme in an attempt to stabilize computations in cylindrical coordinates, though a separate radial filter is also applied in those cases, so this filter was unnecessarily limiting. The 2/3rd rule requires that the highestfrequency 1/3rd of the Fourier spectrum be rendered equal to zero to guarantee that aliasing is not occurring (aliasing is the erroneous appeaerance of energy in highfrequency modes due to insufficient resolution in a Fourier transform. *** Users should contact Dr Greg Sheard if this causes problems in SE/Fourier calcs ***
 Added a new command "ENERGYF", which outputs to a text file the integral norm of the energy in each Fourier mode of a spectralelement/Fourier computation.
23 September 2008
 Added a command SAMPLEF, which can be used to output the Fourier coefficients of the velocity field at a point in the computational domain during SE/Fourier computations.
 Altered the behaviour of the RAND command in SE/Fourier computations  now the imaginary components of the nonzero modes are also randomised.
16 September 2008
 An immersedboundary feature has been added to the code  this experimental feature currently synchronises with an external algorithm that tracks a boundary immersed in the flow.
9 September 2008
 A memory deallocation bug was preventing repeated calls to "init" since the OpenMP parallelization of recent months. This has been fixed.
8 September 2008
 A bug was reported with the Arnoldi routine  on Linux systems the code crashes while writing Tecplot files. The bug was isolated to an OPTIONAL argument of the Tecplot output subroutine, which was not being treated properly. The Arnoldi routine did not supply this argument, causing its value to be undefined. This has now been fixed.
 A bug was reported in the cylindricalcoordinates stability analysis solver, where an incorrect function argument was causing the code to fail without explanation. This has been fixed.
14 August 2008
 Fixed a bug in the SE/Fourier implementation of the "forces" command  previously zero forces were output  now forces are computed from the fundamental mode of the Fourier expansion  this is not currently correct for the cylindrical formulation.
13 August 2008
 Added a Boussinesq approximation for buoyancydriven convection to the 2D and 3D hexahedral solvers. This feature can be activated using the BOUSS command, and requires that a scalar field is being computed, which acts as the temperature field.
6 August 2008
 Added a counter to the output of the step command. This helps the user identify how much longer the currently requested steps have to complete.
 Fixed an error relating to flow at the axis which was affecting the axisymmetric solver. This error had emerged in the restructuring of the code during the SE/Fourier implementation. The code now strongly enforces zero Dirichlet boundary conditions on the axis for the v & wvelocity components.
4 August 2008
 Restructured calculation of advection substep throughout the code. Now the velocity fields are extrapolated to time t+dt, then the RHS of the advection term is computed. Previously, extrapolation was conducted from stored RHS vectors from previous timesteps. This approach reduces memory overhead slightly, and also simplifies the code when simulations are initialized and restarted.
1 August 2008
 Added a command, MASK, which multiplies a velocity field by a specified function. This command can be used to filter or scale velocity fields.
30 July 2008
 Two new commands have been added. MESHPTS outputs the global mesh coordinates to a text file, and VISMAT (under development) will generate image files showing the structure of the matrices being solved by the code.
16 July 2008
 Added OpenMP constructs to advection, pressure & diffusion time integration routines, with conditional "IF" statements to be called only if a single field is being computed. This should improve scalability of simple two or threedimensional base flow computations.
12 July 2008
 Tests have been showing that ROTATE was not functioning correctly. A test case was established whereby a vortex flow with and without a bulk rotation was analysed for threedimensional stability both with and without the ROTATE command being called. This enabled the functionality of the command to be corrected. Tests on unequal strength vortex pairs will continue.
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 nonzero wvelocity, 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 wcomponent of velocity in the base flow (previously only w=0 flows were implemented).
13 June 2008
 Found a bug in the ROTATE routine  the xdirection advection RHS vectors were being copied into both the x and ydirection 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 floatingpoint 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 perprocessor basis. Thus far have sped up PARDISO computations, but other parts have slowed  may need to create perthread 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 thetaderivatives 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 spectralelement/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 outofplane 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 spectralelement/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 timedependent, 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 timedependent (to be evaluated during time integration), and which are timeindependent (to be evaluated during initialisation). This will be incorporated during the rollout of the SE/Fourier algorithm.
8 May 2008
 Rolled out a provisional version of a spectralelement/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 wvelocity component active were not written with the wvelocity 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 NavierStokes 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 zderivatives 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 wvelocity component of the currenttime 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 highorder Neumann boundaries for pressure. These will be implemented shortly.
 Successfully implemented a speedup in the evaluation of timevarying 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 vortexinteraction simulations.
30 September 2007
 Modified the source code varibale naming of wvelocity in the perturbation field, to acknowledge that the wvelocity field is imaginary when a zvelocity 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 changeofvariables 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 zvelocity 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 floatingpoint 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 threedimensional 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 zvelocity component.
27 September 2007
 Stabilized finitezvelocity Flqouet analysis in cylindrical coordinates: problem lay in the changeofvariables 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 nonzero velocity was present in a 3component 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 nonzero 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 nonzero wvelocity 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 wvelocity were identified; relating to products of the wvelocity and zderivatives of the perturbation. These have been added, and tests are underway to verify the solutions.
12 September 2007
 Found and rectified a memory allocationrelated bug that emerged during initialisation for Floquet stability anaysis. This bug was caused when uninitialised CSR matrices were being DEALLOCATED. For some reason, they were being initialised with arbitrary sizes, etc., and this was causing runtime 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 2by2 strain matrix was being allocated, whereas the code requires a 3by3 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 userdefined 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 byzero, 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 reinitialisation.
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" preprocessor "__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 (DDAY)
 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 timevarying creeping flow computations, with timedependent boundaries, for example.
29 May 2007
 On Itanium2 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 floatingpoint 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 "dividebyzero" in strain directional component evaluation.
23 May 2007
 Added a routine that finds and outputs both the physical location, and the value, of a userspecified 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 'dvdxdudy')
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 floatingpoint overflow and divide by zero checking in "gll.f90" and "tec_out.f90" routines to eliminate some runtime 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 timevarying functions makes it far more powerful and userfriendly 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 reverseengineer 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 skewsymmetric 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 vvelocity 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 MONITORassigned 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 zderivatives was not included for cylindrical coordinates. This affected most spatial derivativebased 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 80odd 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 zderivatives might have been miscalculated if WVEL and FLOQ were both active.
 Fixed bug in SAVE routine where wvelocity component was not being saved correctly due to a use of the nearobsolete "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 Helmholtzianstyle inclusion of the zderivative 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 nonzero w velocity fields. This will be used for Floquet analysis of vortex interaction studies where the vortices have nonzero axial flow in the zdirection.
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 2nddimension 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 advectiondiffusion.
 Implemented a loworder scalar advectiondiffusion scheme. Scheme uses a 2ndorder predictorcorrector scheme for advection of the auxiliary equation used to provide the scalar field at the "departure point" of the Lagrangian discretization of the advectiondiffusion equation. The Lagrangian form of the transport equation is a simple diffusion equation, which is solved using a firstorder implicit scheme. A variable number of advection steps per diffusion update are permitted.
9 November 2006
 Detected a bug (through IA64 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 nonzero 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 IA32 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 IA32 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 crosscheck.
6 November 2006
 Added the ability to impose uniform pressure gradients (in x, y & zdirections) via the PGRAD command. The pressure gradients are imposed during the advection step, and are also included in the estimate of the highorder homogeneous pressure boundary condition. This facilitates an easy method for computing pressuredriven periodic flows.
4 November 2006
 Removed obsolete routines which provided RungeKutta stepping for the initial 2 time steps  now exclusively employing 3rdorder 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 backwardsmultistep 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 advectiondiffusion field implementation.
24 October 2006
 Implemented userdefined 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 intermesh 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 intermesh reinterpolation for LOAD routine.
 Fixed a bug in the polynomial order reinterpolation routine by reducing ALLOCATE calls.
20 October 2006
 Finished crossmesh 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 userspecified variables facility in which runtime errors were encountered when no userspecified 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 userdefined 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 userspecified 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 nonNewtonian solvers for consistency. NonNewtionian computations are invoked whenever a nonlinear viscosity is entered in the "viper.cfg" file through the "gvar_nnvisc" function.
 Added an unmissable warning message to the nonNewtonian solver to clearly identify when zero or negative viscosities are computed, which invariably cause simulation instability.
29 June 2006
 Added userdefined 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 restarting 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 skewsymmetric 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 steadystate 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 1storder 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 userdefined 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 nonNewtonian viscosity.
 Repaired nonNewtonian 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 RungeKutta 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 di 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 2ndorder RungeKutta algorithm to compute compressible NavierStokes equations modified to incorporate upwinding in the differential equations being solved by means of an "acousticconvective 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 (nonzero wvelocity evolved with zdirection 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 xaxis as the axis of symmetry, with y being the radial direction.
8 April 2006
 Extended periodic boundaries to 3D  still xdirection only.
7 April 2006
 Rewrite of boundary definition code successful  new "btag" format now required in "viper.cfg" files.
 Planning to redesign boundary definitions in viper.cfg file: user specifies 2digit 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 userdefined static and transient Dirichlet boundaries, periodic and symmetry boundaries, respectively.
 Added periodic velocity boundary condition (restricted to 2D, xdir only at this stage).
 Added pressure "p" as a variable for transient userdefined boundary functions.
5 April 2006
 APAC simulations reveal that the current version of the code has greatly improved temporal characteristics on an initialconditionindependent 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 skewsymmetric. Skewsymmetric 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 1storder initial step due to the incorrect previous time step velocity fields. Two virtual velocity fields at times t0dt and t02dt have been constructed using an RK4 backwards time integration of the advection and diffusion parts of the NavierStokes 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 firstorder temporal convergence of velocity despite 3rdorder time integration.
 Implemented and tested 3D particle tracking  results appear promising  multielement 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 midcomputation.
 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 4thorder RungeKutta scheme to advance the advection term for the first three time steps, to accurately fill the previous step velocity vectors required for the 3rdorder backwards multistep scheme.
 Fixed a small indexing bug in pressure boundary condition evaluation (could this be the 1storder 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 nonDirichlet velocities to simulate a change in reference frame motion (e.g. a body coming to rest, or accelerating).
12 March 2006
 Implemented circular arcbased 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 1storder 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 suboptimal temporal convergence often observed (convergence often little better than 1storder rather than 3rdorder).
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 crosssection!
 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 pressuredriven flows on some meshes.
25 January 2006
 Isolated issue with 64bit 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 userdefined 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 highorder pressure boundary condition. The explicit treatment of the irrotational part created instabilities which were especially apparent for pressuredriven flows.
18 January 2006
 For lowReynolds 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 highorder pressure boundary condition. This is the complete form of the term, and does not require a divergencefree field to be valid.
14 January 2006
 Added time derivative term to highorder pressure boundary condition. Improves accuracy of timevarying 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 nonNewtonian nonlinear 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.
1531 December 2005
 Replaced 3rdorder AdamsBashforth advection / CrankNicolson with theta mod diffusion with 3rdorder stifflystable 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 2ndorder 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 & nonNewtonian 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 nonNewtonian computations  shear rates now taken from nonlinear viscosity.
6 December 2005
 Implemented a 3D version of the nonNewtonian solver (not yet tested).
 Ironed some wrinkles out of nonNewtonian version of code  small errors in shear rate determination and employment of linear reference viscosity. Stability and physical appearance of solutions improved.
 Employed a runtime parsing of relationship for nonlinear 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 nonNewtonian time integration scheme, which splits the viscous term into a nonlinear 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 userdefined 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 nonNewtonian 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 "slipwall" 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 "slipwall" 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 "slipwall" (zero normal velocity) boundary condition (number 4).
 Set up routine to precompute 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 pressuredriven 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 reinitialised successively.
 Tried a "tecp" speedup to no success.
 Implemented toggle between convective and rotational forms for acceleration terms.
 Identified nonsmooth 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 doubleprecision 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 userdefined boundaries.
 Detected problem with userdefined Dirichlet boundaries (showing up on reverse face of elements as well)...
 Fixed bugs in facebiased 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
 Multiface 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 (singleface) 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.