j-c-cook / cpgfunction Goto Github PK

Keeping track of what changes happened when is good practice. See https://github.com/j-c-cook/pygfunction/blob/master/CHANGELOG.md for example.

https://people.sc.fsu.edu/~jburkardt/cpp_src/lagrange_interp_1d/lagrange_interp_1d.html

https://github.com/rncarpio/delaunay_linterp

Write you own:

Using BLAS/Linpack write your own linear interpolation of matrices. https://en.wikipedia.org/wiki/Linear_interpolation

https://scicomp.stackexchange.com/questions/20021/blas-lapack-or-atlas-for-matrix-multiplication-in-c

Matrix matrix multiplication: http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html

This does matrix multiplication on an upper or lower triangular matrix... Does that mean I could use my reduced matrix directly in matrix multiplication? : http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaf07edfbb2d2077687522652c9e283e1e.html#gaf07edfbb2d2077687522652c9e283e1e

This looks promising: https://github.com/parsiad/mlinterp. He claims to have a fast 1D interpolation in C++ of arbitrary dimensions.

The H matrix is "squashed" yes, but creating a H * R matrix multiplication where R is the reciprocity matrix would setup for performing matrix operations.

The matrix is psuedo symmetric by reciprosity.

Collapse any usage of hij_dt down into where it is used.

Transition all existing linear algebra in the library to make use of the jcc::la API.

There may be a bug when use_similarities is set to false.

More than likely, I did not change the h_ij matrix to the new indexing scheme for use_similarities = false.

Currently, each multi-threading call has a commented out non-multi-thread call below it. Rather than relying on commenting and de-commenting, provide a boolean option.

• Run a range of tests, store the g-function output in .json
• Create testing scripts that run those exact g-functions
• open .json and verify that the g-function at each place in time falls within 1.0E-06

Implement a struct that provides quick access to L, Open, U, rectangles or custom borefields.

The coordinate generation functions should be separate from a borehole; perhaps a coordinates.cpp could be created.

No reason to create a whole new variable derived from h_ij. The variable h_ij is a horrible use of memory as is, let there only be one cube.

Is this function being called a lot?

Because the Borehole argument is being passed by value, which is making a copy of the structure in each invocation of it. The argument most likely needs to be a const reference.

The segment response matrix has been transformed from a triangular matrix to a vector.

Currently, the object only gets the h value (get_h_value). A function for storing the h value needs to be made. The code currently exists at the following location:

To reduce memory consumption, the segment response matrix has been packed, and is of size n * (n + 1) / 2. Find a technique to recover the entire matrix using a low level linear algebra like function.

The four components of the hash will be:

• time value
• distance
• difference of burial depth
• difference in height

A g-function interface struct is going to be necessary for EnergyPlus to interact with this library.

• borefield : struct
  • A borefield object needs constructed that provides easy interface to custom, L, U, Open and rectangular fields
  • Handles uniform height (in meters) H, uniform burial depth (in meters) D, uniform borehole radius (in meters) r_b
• time : vector
  • There are different geometric time routines available. EnergyPlus could call the same routine based on an input duration (in months) and optional number of time steps nt
• alpha : float, optional
  • The ground thermal diffusivity (m2/s)

The determination of the number of segments nSegments will be handled by the adaptive discretization scheme. If desired, the ability to turn off adaptive discretization and use a custom number of segments could be made possible.

Allocation of multidimensional vectors is a hindrance on performance. For speed performance, it is vastly superior to define multi-dimensional matrices in a 1D format.

q_reconstructed and SegmentResponse.h_ij are used in _temporal_superposition. Now that #30 is taking place, to obtain a speed increase without allocating (significantly*) more memory, the vectors used in the function need to be 1D.

*Note on significantly more memory: With this technique, it appears that a vector of size nSources * nSources will need to be declared to hold intermediate results. There is a balance between memory and speed.

The _borehole_segments function discretizes all the boreholes into segments of equal length. This occurs inside of the g-function calculation, so it is not feasible that a user calling the function from EnergyPlus could create a field of unequal segments. However, this solution does vary based on a grid. At some point a developer of this library could decide to adaptability vary the segment lengths based on proximity to the center of the field, so that the solution can be both converged and require less segments. At that point in time, it will be beneficial if that developer does find that the g-function calculation is accurate for unequal segment lengths in the field.

Commit fcc6ab5 brought this idea about, and is mentioned in #30 (comment).

The library makes use of a reduced segment response matrix that is inside of a struct

There used to be an h_ij vector that stood alone and was passed around. It needs removed.

The g-function calculation returns garbage results if duplicate boreholes are described in a borefield vector. Implement the find and replace duplicates functions.

MassimoCimmino/pygfunction#81

It would be nice if cpgfunction was both wrapped in Python and available via a pip installation through pypi.

This project is dependent on BLAS, which is used inside of LAPACK. Currently, the only BLAS routine used in cpgfunction is dgesv_ (though there will be more routines used in the future as specific sections of the code are transitioned from multi-threading to vector linear algebra). pygfunction makes use of the LAPACK routine _gesv.

At first, the Python matrix inversion in np.linalg.solve() was much faster than what I was making use of in C++. The root cause turned out to the use of Netlib BLAS library. Based on comparison of timing tests, it is possible the numpy package queries the chip-set to determine which routine is best optimized based on the structure of the chip. On my Intel i7, I was able to get matrix inversion speeds as fast as np.linalg.solve() by making use of the OpenBlas library.

I was reminded of this point when our UGRA Timothy West began computing g-functions on the cluster. He is building on the cluster using CMake, which currently makes use of the Netlib BLAS library. The submission scripts Timothy submitted were failing due to time limits. The reason why the jobs on the cluster can fail based on run-time are discussed in the following section.

This describes the creation of submission scripts for use on the cluster. The path to a directory of configurations is supplied. A wide range of timings have been previously recorded, and a function has been created that takes in the number of sources, and returns estimated amount of time to run. Then, based on those timings, the ideal subset-sums are found based on a given input number of days N. Typically, I will request each submission "bucket" to be one less day than what I am submitting to on the cluster, so that I know the jobs will not run out of time.

A possible outcome of this issue is an automatic selection of the fastest BLAS library, though that may not end up being the case. The following libraries will be tested, and the fastest will be requested in CMakeLists.txt:

• Netlib BLAS
• OpenBLAS
• ATLAS

The g-function calculation code currently computes g-functions using uniform borehole wall temperature (UBHWT) with accuracy. A set of test cases need to be computed and stored. The output of these test cases can be compared to the output of any code that has been modified on commit.

With the merging of #30, the code now has memory bloat. The following lines create a 1D segment response matrix, H_ij from the 2D SegmentResponse.h_ij. H_ij was constructed in this way to complete #30, however, the declaration of H_ij introduces nSources * (nSources + 1) / 2 * nt more memory than is required.

The response matrix SegmentResponse.h_ij has many dependencies throughout the code.

Currently, the library makes use of double precision for every variable. It is feasible that single (float) precision will produce results with high enough accuracy; this could result in a 1/2 memory reduction.

A typedef declaration could be made which applies to the entire library. This would take slightly more effort than simply overloading functions. One thing to keep in mind is that sgesv would need to be called when functions for single precision were overloaded. The call to the external lapack dgesv (double gesv) "C" function occurs here.