Quasicontinuum Code Revision History

Version 1.4 (November 2011)

A number of bugs were fixed:
1. Fixed a bug that could cause some crystals to be built incorrectly when the orientations of the global axes are along high-index crystal directions.
2. Corrected a bug that could cause the connectivity of repatoms to elements to be incorrectly computed for some rather special mesh designs.
3. Corrected a small bug in the calculation of the load-displacement curve for the nano-indentation example.
4. Corrected parsing of input lines in qcmain.f. Blank lines are no longer placed in qc.cmd
5. Corrected parsing of input lines in mod_head.f. Old code could crash if blank lines appeared in qc.cmd
6. Allocation of grad and bsave corrected in mod_solve.f. Old code would crash if "solv,nr" was called after a "solv,cg".
7. A bug in mod_repatom was corrected that results in a segmentation fault whenever the size of arrays used to find neighboring elements in BuildSeenElementsList has to be extended.
8. A bug in mod_repatom was corrected that occasionally led to missed atomic neighbors when there are multiple grains.
9. Bug in mod_plotting fixed. The old code did not correctly generate strain patch plots (no strains appeared in the file).

Two new features were added:
1. Added a macro "mark" that sets or returns to a marked position where the displacement field is stored.
2. Added a convergence test on the number of nonlocal nodes ('nnln'). This item can be used together with the previous item to back-up in the event of a nucleation event and to restart from the displacement field prior to the event.
```
loop,load,40
   bcon
   mark,set
   loop,,10
      mark,go
      solve,cg,1.0,500000,1
      status
      conv,nnln
   next
   solve,cg,1.0,500000,1
   plot,disp,gbshear,0,1.,1.
   report
   pdel,,p-delta
   restart,write,shear
   time
next,load
```
Several optimizations were made:
1. The routine to determine in a point lies inside a given polygon (PointInPoly) was optimized.
2. The routine that builds the connectivity array storing all elements touching a given repatom (BuildIrel) was optimized.
3. The adaption routine were improved so that the refinement of elements on which boundary conditions are applied is more robust (improved routine "fixcon").

Version 1.3 (May 2007)

The QC Code was largely rewritten and restructured to conform to the Fortran 90 standard. As a result the code is very portable. It has been tested with ten different Fortran compilers. In all cases it compiles and runs without modifications. Some of the main changes introduced are:
1. All functions and subroutines were placed within modules. This ensures explicit interfaces for called subroutines and functions.. As a result, a variety of potentially non-portable Fortran 77 "tricks" no longer work. These have been removed.
2. "double precision" variables were replaced with "real(kind=dp)", where dp is defined in mod_global as "selected_real_kind(p=15,r=307)". This ensures that the same precision is used on any machine. The change to "dp" means generic routines must be used such as sin(x), sqrt(x) and not the specific double precision routines dsin(x), dsqrt(x), etc. Also, type change to "dp" precision is done with real(n,dp), where n is an integer, and not dble(n).
3. Obsolete features in Fortran 90 were removed. For example, "character*n" were changed to "character(len=n)". Arithmetic if statements removed, etc.
4. Passing of array sections to subroutines were replaced by passing a temporary array containing the relevant section. Technically, this is not necessary since Fortran 90 allows passing array sections. However, in practice some compilers failed to do do this correctly especially, when the array was part of a derived type.
5. Dynamic memory allocation of variables was done automatically whenever possible, otherwise the "allocatable" clause was used. The use of pointers was limited to cases where this was unavoidable.
Greatly improved streamlined make architecture designed to allow easy porting from one platform to another. The new approach is very easy to use and robust. It has been tested on a variety of Unix and Linux systems as well as OS X 10.4.9.
The code was optimized in a number of ways. A speed-up of about 30-40% was observed for typical applications when switching from version 1.2 to 1.3.
The memory requirement of the code were greatly reduced. This was done by doing the following:
1. A number of memory leaks were detected and fixed.
2. The irel() array was reduced from storing all neighbors of a repatom within the cutoff radius to only its nearest neighbors.
3. The arrays xl() and ul() were removed.
Restart file size was greatly reduced due to the changes in item 4. Note that version 1.3 will still read in restart files from older versions of the code. This is done automatically without user intervention.
Numerous bugs fixed. Most are minor and will not be mentioned. The more major ones are listed below. The subroutine containing the bug appears in parenteses.
1. (StatusCalc) The calculation of deformation gradients and eigenvalues could be corrupted after adaption. As a result repatom status (local or nonlocal) could be incorrectly computed. This would be corrected in subsequent steps but could explain some "jerky" behavior that was sometimes observed in simulations. This was fixed. Note: The corrected code gives results that can be different from earlier versions because the set of nonlocal atoms can be different.
2. (pmacr) When calling restart b() was stored in blstat() for 'read' but not for 'txrd'. This was fixed.
3. (delaunay) Bug related to convex hull triangulation for the special case where ncb=0 and nce/=0 fixed.
4. (GetCellData) For "bad" crystallographic orientations with very large repeat distance the algorithm to find the repeating cell could fail. This was fixed.
A number of changes were introduced to the structure of the required user routines that interface with QC. The changes make the routines easier to work with and safer (i.e. less likely to "break" QC due to any errors in them). The main changes are:
1. All user routines are now in a single file called "user_APP.f" where APP is a particular application (e.g. user_gb.f). This file contains:
  
  user_mesh -- mesh generator
  user_bcon -- boundary conditions
  user_pdel -- force versus displacement
  user_potential -- user-specified external potential
  user_plot -- user-specified plot option
  
  The last two are new in this release. See the reference manual for details.
2. Direct access to many module variables is no longer possible to prevent a user inadvertently "breaking" the code. Instead routines were added to allow a user to query a module for the value of a particular variable. These routines include:
  
  In mode_grain:
  
  NumGrains - return number of grains
  GetGrainCellSize - Get the repeating cell dimensions for a given grain
  GetGrainNumCell - Get the number of atoms in the repeating cell for a given grain
  GetGrainCellAtom - Get coordinates of a given atom in the repeating cell for a given grain
  GetGrainRefStiff - reference stiffness matrix for a given grain
  GetGrainBvec - Get Bravais vector matrix for a given grain
  GetGrainBinv - Get inverse Bravais vector matrix for a given grain
  GetGrainRefatom - Get reference atom for a given grain
  GetGrainSpecies - Get species for a given grain
  GetGrainNumVrts - Get number of vertices for polygon of a given grain
  GetGrainVertex - Get given vertex of the polygon of a given grain
  GetGrainXlatvect - Get crystallographic axes for a given grain
  
  In mod_pload:
  
  GetLoadTime - Get current loading time
  GetLoadTimeStep - Get current loading timestep
  GetLoadPropFact - Get current loading proportional factor
  GetLoadOldPropFact - Get previous loading proportional factor
3. "only" clause added to all use statements.
4. Calling header to subroutine qsortr() (called from user_mesh) was changed to:
  
  subroutine qsortr(n,list,xkey,sign,listdim,keydim)
  
  The variables 'ndim' and 'ind' that appeared between 'xkey' and 'sign' in the old header were dropped.
A file user_template.f is included in the QC directory as an aid to porting user routines from older versions to version 1.3. The file contains a skeleton of the file along with instructions on where and how to place the older user routines.
The version 1.3 release is focused more on a restructuring of the code as a basis for future development. However, there are a number of minor extensions in this version:
1. Two user routines user_potential and user_plot were added, allowing the user to define a user-specified external potential and plot option. See the reference manual for more information.
2. The zone macro (formerly used to defined no-adaption zones) was generalized to describe "special-attribute" zones. This currently includes nonlocal zones and no-adaption zones. See the reference manual for details.
3. More input verification was added to the code to catch user input errors.
Version 1.3 is downward compatible being able to read input files used with earlier versions of the code.
The QC Reference Manual and QC Tutorial Guide were significantly expanded. The Reference Manual, in particular, has an entirely new chapter on the structure of the QC code and on writing the user routines.
The utility program "dynpots.f" was updated from f77 to conform to the f90 standard.

Version 1.2 (February 2005)

Support for multiple unconnected meshed domains (This was included in version 1.1 - but the new version is much more robust). A new example "Friction-example" is included to demonstrate this capability. (RM)
mod_crystal: Small bug in CryCount that would not have affected performance was corrected (ET).
New flag DummiesAllowed added. When set to .false. no dummy atoms are allowed, meaning that all neighbors of nonlocal atoms must be repatoms themselves. If the flag is set to false and dummies are detected in a status calculation, adaption is triggered until there are no more dummies. If using the NR solver no changes are necessary to the input file besides setting the appropriate flag. For the CG solver the solver must be placed inside an additional loop with a status,update command:
```
loop,,5
solve,cg,1.0,10000,1
status,update
convergence,force
next
```
The reason is that neighbor lists are updated inside the CG solver. This can result in dummy atoms being generated. It this happens and dummies are not allowed, the CG solver will exit with a warning. The loop above handles this through the status,update command which will trigger a full status calculation and refine the mesh to remove the dummies. (ET)
Before writing out a restart file, a backup of the last written restart file is saved (with "_LAST" added to its name). For example,
```
rest,write,qc_simulation
```
will copy the file qc_simulation.res (if it exists) to qc_simulation_LAST.res, before writing out the new restart file. The reason for this is that if something interesting happens in a simulation, it is often helpful to have the restart file of the step before the event occurred, in order to restart from there with some changes to the input file, e.g. more adaption, etc. (ET)
A new file called nonstandard.f has been added to the QC code. This file contains any nonstandard Fortran extensions used by the QC code to facilitate portability. Compilers which support these extensions differently from most or that need to use specific modules in order to access these extensions have their own nonstandard.f file (nonstandard_compiler.f) along with a corresponding make file. If you are porting QC to a new platform, modify nonstandard.f to support the new compiler. (ET)

Currently nonstandard.f contains the following routines:
1. flush command.
  
  Flush is a non-standard Fortran command that forces the output buffers to be written to disk. Without it the output files are only written when the buffers become full (the size of the buffers is system dependent). Because the command is non-standard, some compilers implement it differently than others. So far we have encountered trouble only with the SGI MIPSPro Fortran 90 compiler (versions 7.4 and higher) where the number of arguments to flush was changed. As a result the previous version of the QC code crashed with a bus error. (Thanks to Rafael Estevez at GEMPPM INSA Lyon for helping us figure this out.)
2. system command
  
  System is a non-standard Fortran command that lets the program issue an operating system command while running. The system command is used to make the backup copy of the restart file described above.
Various code optimizations were performed resulting in about a 25% speed-up (depending on the platform) (ET). These optimizations include:
1. The function dot11 which was called by the active column solver was removed and replaced with faster inline code.
2. The function trarea which computes the area of a triangle was eliminated and replaced with faster inline code.
3. The function intri which determines whether a point is inside a triangle was optimized. The part of the routine that computed area coordinates and determined whether the point was on the edges of the triangle was removed because this information was never used.
4. The routine fe_locate which finds out which element a given point lies in was optimized.
5. Some service routines in qclib (Dist2, distance, etc.) were removed and replaced with faster inline code.
The initial output generated by QC has been extended to include the stress tensor and elastic constants matrix (in Voigt notation) in the reference configuration. The stress tensor should be zero. The degree it differs from zero is a measure of the accuracy of the lattice parameter given to the program. Both the stress and elastic constant components are given in relation to the global c.s. (ET)

To prevent the code from crashing in the event that moduli calculations are not implemented in the potential being used, an internal error code was added to the routines Calc_Local and NonLocal and to routines calling them. This is transparent to the user and only affects developers.
Improved reliability in neighbor finding for nonlocal repatoms. In circumstances of very large relative deformation between two neighboring grains or for grains with a highly non-convex shape (for example with a long thin notch) the code could sometimes miss neighbors to some repatoms, leading to an error in the energy and force calculations. This has been fixed and the new version is now extremely robust at the expense of some slowing down of the code. (RM)
If a model contained grains with different periodicities in the out of plane direction, the code would not compute the neighbor lists of nonlocal atoms correctly. The new version can now handle this situation, although it is still constrained by the 2.5-D nature of the code and results must be interpreted with caution. (RM)
Multiple plots of the same data, for example the displacements after each load step, can now be appended to multiple zones of a single file, rather than multiple individual files. This is achieved with the "append" flag in the "plot" command. (RM)
The z-component of all repatom positions are now mapped onto Bravais lattice sites within the first periodic copy of the crystal structure. This doesn't change the results, but makes visualization of out-of-plane deformations easier. (RM)
The "nonlocal fully-refined" feature of the adaption routine is now more robust in that it is less likely to find circumstances where it cannot achieve the full refinement of nonlocal repatoms. In the event that it fails, it terminates with an error message. (RM)
The input file is now much more flexible in terms of accepting spaces at the start of commands and comment lines (denoted with a "%") are now accepted everywhere in the file. (RM)
A new option, 'type', has been added to the 'plot' command that plots the positions and atomic number of each repatom. (RM)

Version 1.1 (May 2004)

qclib.f: subroutine fe_locate was rewritten to allow for multiple separate meshed regions. Previous versions were limited to a single connected domain. (RM)
mod_pload.f: The proportional load table was modified to handle a step function loading without crashing (two prop values that have the same time value). (RM)
adaptgr2.f: Function GradientNorm() was modified so that an element with a nonlocal node is automatically targetted for refinement even if the centroid of the element lies inside a no-adaption zone. This is necessary to prevent nonlocal nodes at the edge of no-adaption zones representing more than themselves and generating an internal QC error. (corrected by ET following a bug report by Markus Beuhler)
mod_grain.f: (1) Target attribute added to grains so that pointers to grain data can be defined later in subroutine nearestbravais in polygrain.f. (2) Components of reciprocal lattice in global c.s. computed and stored for each grain to save computation time later in subroutine nearestbravais in polygrain.f. (changes introduced by Jon Wilkening, Courant Institute of Mathematical Sciences)
mod_repatom.f: Data structures used by subroutine FindNeighbors to store Bravais indices visited in neighbor search and dummy atoms were changed from linked lists to pointers into vectors. The resulting code is significantly faster. An additional speed-up was obtained by improving the initial guess for the element breadth search used to locate the element containing a specific atom. These changes will speed up status calculations, neighbor list recalculations triggered in the solver and restart operations. (changes introduced by Jon Wilkening, Courant Institute of Mathematical Sciences)

NOTES

(1) The new initial guess used by the element search routine (fe_locate) can result in slightly different results from those obtained using previous versions of the code. The reason for this is that nodes lying on the interface between two elements will sometimes be assigned to a different element that in previous versions. As a result, due to numerical precision issues, the energy and forces can be slightly different. Because the solution procedure in the code is iterative, even very small changes in forces can result in more pronounced changes after many iterations. Qualitatively, the results are expected to be the same, but small quantitative differences may be expected.

(2) A separate mod_repatom.f file for the Portland Group Compiler (PGF) (mod_repatom.f.pgf) is no longer necessary. The new mod_repatom.f file should work for PGF as well.
"character*xx" statements were replaced with "character(len=xx)" throughout the code, to prevent warnings in some fussy compilers. (ET)
mod_poten: Bug corrected at end of shell loop in subroutine Calc_Local. In the old version an incorrect error message was given in the rare case where rinf=rcryst and the simulation was terminated. Also the if statement checking whether to continue the loop did not work correctly at the limit where all shells were visited. Both issues were corrected. (ET)
mod_solve: a "call output" was added in subroutine cgstep to inform the user when the neighbor lists were updated during a CG step. The updates were carried out in the previous version, but no message was printed in the output stream. (ET)
mod_poten: a bug which caused the code to crash when multiple potential files were read in was fixed. The problem had to do with an incorrect print statement. (ET)
mod_poten: error handling was added so that if an incorrect potential format was specified the code would exit gracefully with a meaningful message. (ET)
refine.f: this file containing a duplicate definition of function refine was removed. (ET)
plotting.f: A time and property stamp were addded to the title of Tecplot plot files generated by the plot macro. This allows the user to know which load step he or she is looking at from within Tecplot. (ET)
Makefiles added for the Intel Fortran Compiler version 8.0. (ET)

Version 1.01 (August 2003)

mod_repatomb.f: a minor bug was corrected (incorrect loop upper bound on line 379).
mod_adaptgr2.f: a minor bug was corrected (inconsistency in the number of arguments in subroutine constraints on lines 248,252).
qcmain.f: a minor bug was corrected (an extraneous comma in the QC banner print statement was removed on line 157).
A bug in the input file for the GB example (gb_shear.in) was corrected. The crystallite radii for the grains was changed from 3.0 to 20.0 (lines 42 and 58).
Make files for the Portland Group F90 compiler were added.

Note: The Portland Group F90 compiler generates error messages when compiling mod_repatomb.f. This is because it does not recognize the NULL() statement used to nullify pointers. An alternative routine mod_repatomb.f.pgf has been included in the distribution which solves this problem. If you are using the Portland Group compiler you should replace mod_repatomb.f with this file.

Version 1.0 (February 2003)

Original release.