[WIP] New particle interpolation schemes

[hlokavarapu]

Added bilinear and biquadratic least squares method as particle property interpolation schemes with a simple overshoot and undershoot correction. Mainly, a minimum and maximum property value of all local particles within a cell is stored. If the interpolated value is larger than the local maximum or is smaller than the local minimum, then the interpolated value is set to the local maximum or local minimum property value, respectively.

Contributed by @yinghe616, @egpuckett, and myself.

Part 1 of #1319

[gassmoeller]

Very nice! Please address the minor issues, and resolve the conflict in changes.h.

[gassmoeller]

We recently noticed that double underscores produce a warning in some modern compilers, and removed them from all other header files (they are reserved for compiler commands). Please replace them by single underscores (also in the following line, and in the other files).

[gassmoeller]

Please move this comment to the header file.

[gassmoeller]

Would it be easy to extend this to 3d, and would you be willing to extend it? If not that is totally fine, but then please add this limitation to the entry in changes.h, and to the description of the class, so that users are aware of it.

[gassmoeller]

Just saw that you mention it already in the description for the parameter file. Still, please also mention it in the class description of the header file, and the entry in changes.h.

[hlokavarapu]

For now, I'll have only implemented the 2D case. The 3D case is in the works and will be added on soon.

[gassmoeller]

I think there is an optimized version of this algorithm in interpolator/cell_average.cc, please use that one if possible.

[gassmoeller]

Same comment as above for the other class.

[gassmoeller]

Same as above.

[egpuckett]

There are several caveats to this PR.

1. The limiter for the bilinear interpolation algorithm can be improved by searching over adjacent cells. In a discussion with Harsha we (I?) decided it was a good idea to ask Rene how best to loop over adjacent cells, since my understanding is that there is a multimap containing particles sorted by cells and a mechanism for knowing nearest neighbor cells. We probably would be able to figure this out for ourselves, but my view is having Rene's guidance / suggestions would be better.

2. The biquadratic interpolation is not performing to design specifications; the error starts to increase at some point when other interpolation algorithms such as the bilinear algorithm will continue to decrease the error. Thus there is some debugging that needs to be done. Perhaps having more eyes on this problem will help.

[gassmoeller]

/run-tests

[bangerth]

About point 1: If you are sitting on cell, then you can ask for its (face) neighbors by writing cell->neighbor(f) where f is the number of the neighbor you care about (0..3 in 2d, 0..5 in 3d). This is again a cell, and you can access the particles on that cell (via the particles map) in exactly the same way as you access the particles on the current cell.

About point 2: I'm not sure how much anyone other than you guys can help with the debugging. Try to come up with as simple a problem as you can that shows the lack of convergence. Print output that shows things you think ought to be true (e.g., that the biquadratic interpolation's values should be within a certain range of values) with what actually is true in the implementation. Eventually you will come to a point where the implementation is doing something you know is not correct. Then figure out why this is so by going backward in the code from the point where something was already wrong.

[gassmoeller]

Very sorry for coming back so late to this PR. Could you rebase to a current master and resolve the conflict in chages.h?
About the caveats:

1. Sounds like an improvement that does not necessarily need to be included here, and can be added later. Let us not delay this PR, but rather open a new one when the feature is ready. About the implementation: Wolfgang's suggestion is a first step, but it will only catch face neighbors (not vertex or line neighbors), and you will need to check for children of the neighbors in case the cell is refined. A more complete approach is to use the deal.II function GridTools::vertex_to_cell_map to generate a map between all vertex indices and their adjacent cells, and check all cells that are adjacent to any of the vertices of the cell you are currently interested in. Be careful to only call that function once, and not once per cell (this would kill the scaling of the algorithm). A different application, but using this function is done in source/particle/world.cc:804-870
2. I agree with Wolfgang, debugging without a lot of information is hard, and you are probably in the best position to find the problem. Maybe some questions can help: Is the limitation dependent on the model setup (i.e. is it caused by something like SolCx benchmark with compositional fields is not as accurate as without compositional fields #1355)? At which point does the error start to increase (i.e. is it at 1e-2 at refinement level 3, or at 1e-13 at refinement level 8, which likely is just a numerical inaccuracy that might never appear in practice)? Following Wolfgang: Can you simplify the problem to a point, where only the interpolation is left, and you can compare results from each step of the algorithm to analytical solutions?

I am not too much concerned going forward with this PR as long as 2. does not seem to interfere with practical applications.

[hlokavarapu]

@gassmoeller , as you probably know, in both the bilinear and biquadratic interpolation scheme, the least square method involves the inversion of a matrix. I use the deal.ii algorithm gauss_jordan() to do the inversion. However, for a certain 2D model where I generate the particles using reference cell such that there are 16 evenly spaced particles per cell, I hit this assertion for the 10th or so call to properties_at_points

An error occurred in line <1810> of file </home/harsha_lv/CIG_codes/deal.II/include/deal.II/lac/full_matrix.templates.h> in function
    void dealii::FullMatrix<number>::gauss_jordan() [with number = double]
The violated condition was: 
    max > 1.e-16*typical_diagonal_element
Additional information: 
    The maximal pivot is 0.000800985, which is below the threshold. The matrix may be singular.

To circumvent this, I have realized that I can run this model in Release mode but should I be concerned? Is this the algorithm of choice in terms of inverting a matrix?

To reproduce this issue:

Particles_Ra_1e5_B_1_k_3.prm.txt

References:
@egpuckett , @yinghe616

[bangerth]

Yes, if something does not run in debug mode, you should always be very concerned. Ignoring errors in debug mode by just running in release mode is the same as just putting a piece of duct tape over the check engine light in your case because you don't like seeing it -- it doesn't make the problem go away.

Is the problem reproducible? If so, can you output the offending matrix? If I recall the algorithm correctly, then a small or zero pivot indicates that the matrix is (nearly) singular. That's of course a problem, but it would be good to check that by looking at the entries of the matrix and verifying the issue.

[egpuckett]

Ignoring errors in debug mode by just running in release mode is the same as just putting a piece of duct tape over the check engine light in your case because you don't like seeing it -- it doesn't make the problem go away.

Well said!

[hlokavarapu]

I believe that I may have found the cause or perhaps a second issue. Debugging this branch, I found the trail of a bug in my hastiness to keep up-to-date with #1384 . Interpolating a single particle property onto a single compositional field works as expected but interpolating two particle properties onto two compositional fields breaks down. I am using the aspect/tests/particle_interpolator_bilinear.prm added to this branch as my initial test case. After resolving this issue, if I then run into the issue of a pivot being close to zero during the inversion of the matrix, I will let you guys know ASAP.

[hlokavarapu]

After double checking that the returned vector of vectors for the function properties_at_points is consistent with the component mask, I fixed issue number 2. I then retraced my steps to the original issue.

For an active cell of level 0 and index 19, the constructed matrix A is

>>> A
array([[  2.37695312e+13,   6.09375000e+11,   9.75000000e+06],
       [  6.09375000e+11,   1.95312500e+10,   2.50000000e+05],
       [  9.75000000e+06,   2.50000000e+05,   4.00000000e+00]])

and taking the inverse,

A_inverse.gauss_jordan();

once again resulted in

--------------------------------------------------------
An error occurred in line <1810> of file </home/harsha_lv/CIG_codes/build/dealii_only_mpi/install/include/deal.II/lac/full_matrix.templates.h> in function
    void dealii::FullMatrix<number>::gauss_jordan() [with number = double]
The violated condition was: 
    max > 1.e-16*typical_diagonal_element
Additional information: 
    The maximal pivot is 0.000656922, which is below the threshold. The matrix may be singular.
-----------

Using python3 module numpy, the inverse of A is computed to be

>>> numpy.linalg.inv(A)
array([[  2.56000000e-10,  -7.15040699e-25,  -6.24000000e-04],
       [ -7.13626344e-25,   2.56000000e-10,  -1.60000000e-05],
       [ -6.24000000e-04,  -1.60000000e-05,   1.52225000e+03]])

and A*A_inv is

>>> numpy.matmul(A, numpy.linalg.inv(A))
array([[  1.00000000e+00,  -9.50048417e-15,   7.26882718e-07],
       [ -1.88841998e-14,   1.00000000e+00,   3.08646122e-08],
       [ -4.33680869e-19,   0.00000000e+00,   1.00000000e+00]])

In this problem, the dimensions of the 2D domain is 2e6 m by 6e6 m. After rescaling the domain such that it is 1 m by 3 m, for the active cell of level 0 and index of 19, the matrix A_rescale is then

>>> A_rescale
array([[  5.94238281e+00,   1.52343750e-01,   4.87500000e+00],
       [  1.52343750e-01,   4.88281250e-03,   1.25000000e-01],
       [  4.87500000e+00,   1.25000000e-01,   4.00000000e+00]])

and the inverse of A_rescale is

>>> numpy.linalg.inv(A_rescale)
array([[  1.02400000e+03,  -5.00289345e-12,  -1.24800000e+03],
       [ -5.00248517e-12,   1.02400000e+03,  -3.20000000e+01],
       [ -1.24800000e+03,  -3.20000000e+01,   1.52225000e+03]])

and A_rescale * A_rescale_inverse is

>>> numpy.matmul(A_rescale, numpy.linalg.inv(A_rescale))
array([[  1.00000000e+00,  -2.66453526e-15,  -1.13686838e-13],
       [ -2.84217094e-14,   1.00000000e+00,   0.00000000e+00],
       [ -9.09494702e-13,   0.00000000e+00,   1.00000000e+00]])

I am not sure what to make of these results. Any guidance would be much appreciated.

[bangerth]

If I understand you right, you are doing a least-squares fit onto a polynomial basis? Since you have a 3x3 matrix, can you say what basis you use? I'm slightly confused because you say you do a bilinear basis, but in 2d that would have 4 basis functions...

[hlokavarapu]

I am attempting to fit the function f(x,y) = c3 + c2*y + c1*x using least squares.

The original reference is from the following paper where they have also excluded the x*y term:

@Article{thielmann2014discretization,
title={Discretization errors in the hybrid finite element particle-in-cell method},
author={Thielmann, Marcel and May, DA and Kaus, BJP},
journal={Pure and Applied Geophysics},
volume={171},
number={9},
pages={2165--2184},
year={2014},
publisher={Springer}
}

[yinghe616]

@bangerth, @hlokavarapu currently is doing a P1 least squares fitting, and return the values to composition field at supporting points with a Q1 basis, so it's essentially a P1

[yinghe616]

@bangerth @hlokavarapu , if I remember correctly, A_inverse.gauss_jordan() only explicitly gives the inverse matrix for matrix size less than 2x2, so for larger matrix size, it's solved implicitly by not forming the inverse matrix? Therefore, for a least squares fitting problem, the matrix may not be invertable or nearly singular, but it can still return a least square solution of the matrix system, which may give a warning and cause the trouble in the debug mode?

[bangerth]

OK, so it's a linear, not bilinear function space. I think that by itself should be fixed -- since we are working on quadrilaterals, it would be more appropriate to use the bilinear space (or the trilinear space in 3d).

As far as the matrix is concerned: The monomial basis {1,x,y} leads to poorly conditioned problems if your variables are poorly scaled (as is the case here). I think a better approach would be to use the basis {1,x-x0,y-y0} where (x0,y0) are the coordinates of the cell center. But even that may be poorly suited; I would imagine that using the basis {1,(x-x0)/h,(y-y0)/h} is best, where h is the cell's diameter.

All of this is a problem because in your case, x and y are variables that are very large, and so the factors c2 and c1 in f(x,y) = c3 + c2*y + c1*x have a very different magnitude than c3. This results in the ill-conditionedness of the matrix. Choosing a different basis, as outlined above, will yield better results because then the coefficients will have roughly the same order of magnitude.

[yinghe616]

@bangerth Thanks for your explanation. It seems well explained the current issue. :-)