Array kernel #13314

YaqiWang · 2019-04-25T21:01:56Z

This is definitely still working in progress. I push it here to let you have an idea on what I am doing, particularly for @friedmud and @lindsayad .

The design here:

use Eigen::Array to hold local dofs for an array variable;
use variable groups in libMesh to group variables in an array variable together;
use template to avoid code duplication with standard and vector variables;

This map is useful for understanding the template:

 * OutputType          OutputShape           OutputData
 * ----------------------------------------------------
 * Real                Real                  Real
 * RealVectorValue     RealVectorValue       Real
 * RealArrayValue      Real                  RealArrayValue

So far, I have finished variables, coupeable, systems, initial condition. What I have not done is kernels, bcs, etc. but should be finished soon.

I also did not:

use Eigen::Map for faster solution vector access;
several switches based on isNodal in MooseVariableFE.C due to the ordering scheme difference between nodal variables and high order elemental variables.

So feel free to take a look. It is not mergeable yet.

YaqiWang · 2019-04-25T21:03:45Z

Hate conflicts...

YaqiWang · 2019-05-01T17:58:52Z

This could be easier than i thought. I do not need to use array type in local residual or jacobian but can instead resizing the dense vector or matrix differently. Just need to managing the index carefully.

YaqiWang · 2019-05-06T17:04:49Z

I made the first cut ;-) I added a simple test in e670cc1 and will be working on more tests. The final look of kernels, bcs are in b75aeab.

moosebuild · 2019-05-06T20:16:22Z

Job Documentation on 8d2e314 wanted to post the following:

View the site here

This comment will be updated on new commits.

YaqiWang · 2019-05-06T21:27:23Z

(Going to write down something here for the design, that should be useful for reviewing.)

An array variable is define as a set of standard field variables with the same finite element family and order. Each standard variable of an array variable is referred as a component of the array variable. An array kernel is a MOOSE kernel operating on an array variable and assemble the residuals and Jacobians for all the components of the array variable. The purpose of having array kernel is to reduce the number of kernels when the number of components of an array variable is large (potentially hundreds or thousands) and to avoid duplicated operations with lots of standard kernels. The array kernel can be useful for radiation transport where an extra independent direction variable can result into large number of standard variables. Similarly as array kernel, we will have array initial conditions, array boundary conditions, array DG kernels, array interface kernels, array constraints, etc.

The design:

to use variable groups in libMesh to group components in an array variable together;
to use template to avoid code duplication with standard and vector variables;
to use Eigen::Matrix to hold local dofs, solutions on quadrature points for an array variable;
to use dense matrices or vectors as standard or vector variables with proper sizes for holding the local Jacobians and residuals.

The following map is useful for understanding the template:

 * OutputType          OutputShape           OutputData
 * ----------------------------------------------------
 * Real                Real                  Real
 * RealVectorValue     RealVectorValue       Real
 * RealArrayValue      Real                  RealArrayValue

The three rows correspond to standard, vector and array variables.

The final form of an array kernel is quite simple. Using ArrayDiffusion as an example. The computeQpResidual function has

return _grad_u[_qp] * _array_grad_test[_i][_qp] * (*_d)[_qp];

where _grad_u[_qp] is an Eigen::Matrix with number of rows being equal to the number of components of the array variable and number of columns being LIBMESH_DIM. _array_grad_test[_i][_qp] is an Eigen::Map of classic _grad_test[_i][_qp] which is in type of Gradient. Thanks to the Eigen matrix arithmetic operators, we can have a simple multiplication expression here. _d is a pointer of a material property of Real type for scalar diffusion coefficient. Here we assume the diffusion coefficient is the same for all components. The returned value is an Eigen vector.

Of course, if we have different diffusion coefficients, we will have something like

    RealArrayValue v = _grad_u[_qp] * _array_grad_test[_i][_qp];
    for (unsigned int i = 0; i < _var.count(); ++i)
      v(i) *= (*_d_array)[_qp](i);

_var.count() gives the number of components of the array variable. RealArrayValue is a typedef of Eigen::Matrix in size of number of components of the array variable.
If we have coupled diffusion terms, i.e. diffusion coefficient is a matrix, we will have

  return (*_d_2d_array)[_qp] * (_grad_u[_qp] * _array_grad_test[_i][_qp]);

Correspondingly the computeQpJacobian has

    return RealArrayValue::Constant(_var.count(),
                                    _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d)[_qp]);

or

    return _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d_array)[_qp];

or

    return _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d_2d_array)[_qp].diagonal();

It is noted that only the diagonal entries of the diffusion coefficients are used in the fully-coupled (3rd) case because computeQpJacobian is supposed to only assemble the block-diagonal part of the Jacobian. The full local Jacobian is assembled in function computeQpOffDiagJacobian, where when the off-diagonal variable is the array variable, we have

  return _grad_phi[_j][_qp] * _grad_test[_i][_qp] * (*_d_2d_array)[_qp];

The retuned value is in type of an Eigen matrix with number of rows and columns equal to the number of components.

Future work:

to profile the code with large number of variables to find the hot spots and fix them if there are any.
to change the current dof ordering for elemental variables so that we can avoid bunch of if statements with isNodal() in MooseVariableFE.C, (refer to Dof Number Ordering libMesh/libmesh#2114).
to use Eigen::Map for faster solution vector access.

YaqiWang · 2019-05-06T21:53:15Z

I think this is ready for review. Tag @friedmud @lindsayad and @permcody . I am working on more array kernels, BCs, tests, but they should not affect the review.

permcody · 2019-05-06T22:38:47Z

Still a fair number of test problems, have you looked at all the build boxes yet? Bighorn for instance is failing nearly every test with the same error.

YaqiWang · 2019-05-06T22:50:23Z

I looked the error messages and am quite certain they can be easily fixed. For example, Bighorn is probably doing something similar like chemical reaction module deriving AddVariableAction, which I changed the type of one of its parameters. This is a new capability so the impact to the tests is actually quite small.

YaqiWang · 2019-05-07T20:23:07Z

Rattlesnake & MAMMOTH should have been fixed.

lindsayad · 2019-05-07T20:58:09Z

This is going to have merge conflicts with #13056 because I've refactored MooseVariableFE to store all its data in MooseVariableData objects. This eliminates the code duplication from having both element and neighbor calculations (and with the new mortar class, lower dimensional calculations).

YaqiWang · 2019-05-07T21:21:04Z

I am ok to resolve the conflict. Should be fairly simple although could look intimidating.

YaqiWang · 2019-05-07T21:27:57Z

I will need ArrayInterfaceKernel and ArrayConstraint for the multi-scheme stuff in Rattlesnake. But since ArrayConstaint will probably depend on #13056 and #13359 could be useful to ArrayInterfaceKernel, I guess I will not have these two in this PR. I plan to add ArrayDGKernel and few more tests coupling standard and array variables with various families. So within two more commits, I will possible be done on this PR.

YaqiWang · 2019-05-08T20:25:59Z

I have not add interface functions in NeighborCoupleable because so far I do not need them. I plan to added in a following-up PR for adding ArrayInterfaceKernel.

YaqiWang · 2019-05-08T20:28:59Z

I think now I have all the functionalities in my mind in this PR. I am still missing documentations. The tests under test/tests/kernels/array_kernels are verified with the inputs without array variables.

YaqiWang · 2019-05-08T22:46:01Z

The major pain on resolving the conflict is in variable and assembly, but looks may not be too bad...

…tended to array variable idaholab#6881

idaholab#6881

…nterface kernels or constraints

YaqiWang · 2019-05-10T09:52:11Z

This is the largest merge conflict I ever resolved, partly due to the bad arrangement of commits (it is very hard to keep a clean history of the commits for a large PR)... I am sure there are rough corners in this PR, but I think it is ready for review. The documentation failure is due to lack of SQA in some existing tests which I touched the tests file for correcting an error message. The failing apps should be fairly easy to fix.

lindsayad

I have some questions, but overall I think this looks relatively unobtrusive for adding a neat new functionality, and in general it's well documented.

lindsayad · 2019-05-13T16:36:46Z

framework/include/variables/MooseVariableData.h

@@ -503,7 +556,36 @@ class MooseVariableData
  void fetchADDoFValues();
  void zeroSizeDofValues();

-private:
+  void getArrayDoFValues(const NumericVector<Number> & sol,


This method is long. I'd like to see it moved to the .C file

lindsayad · 2019-05-13T16:42:36Z

framework/src/variables/MooseVariableFE.C

+void
+MooseVariableFE<OutputType>::addSolution(const DenseVector<Number> & v) const
+{
+  _element_data->addSolution(_sys.solution(), v);


I actually don't know how this compiles...This method is marked const so _sys should be const. That means it would call the const getter for solution, which returns a const reference to the NumericVector. And then we're trying to bind that const reference to a non-const reference...yea I don't know how this compiles!

Indeed. Let me search a little to see why.

Maybe the reference behaves like pointer? Although it is const, it just indicating the reference is not changed but not the content of the reference?

Found an answer: The constness of execute() only affects the this pointer of the class. It makes the type of this a const T instead of just T*. This is not a 'deep' const though - it only means the members themselves cannot be changed, but anything they point to or reference still can.* at https://stackoverflow.com/questions/2523516/why-can-i-call-a-non-const-member-function-pointer-from-a-const-method. Not sure it is the standard (although probably).

Yes, I'm aware of that, but what is the difference between our use case, and the program below, which emits the compiler error that I expect? (constness.cpp:25:62: error: binding value of type 'const int' to reference to type 'int' drops 'const' qualifier)

#include <memory> class SystemBase { public: SystemBase() : _solution(5) {} const int & solution() const { return _solution; } private: int _solution; }; class MooseVariableData { public: void addSolution(int & solution, int y) const { solution += y; } }; class MooseVariableFE { public: MooseVariableFE() { _element_data = std::make_unique<MooseVariableData>(); } void addSolution(int y) const { _element_data->addSolution(_sys.solution(), y); } private: std::unique_ptr<MooseVariableData> _element_data; SystemBase _sys; }; int main() { }

Sorry, yea just like the post said; the program below doesn't compile, but if we change to the reference then the program does compile, and the non-const solution getter is called. That's sad as far as ensuring const-correctness goes 😢

#include <memory> #include <iostream> class SystemBase { public: SystemBase() : _solution(5) {} SystemBase(const SystemBase & copy) : _solution(copy._solution) {} const int & solution() const { std::cout << "Const version\n"; return _solution; } int & solution() { std::cout << "Non-const version\n"; return _solution; } private: int _solution; }; class MooseVariableData { public: void addSolution(int & solution, int y) const { solution += y; } }; class MooseVariableFE { public: MooseVariableFE(SystemBase & sys) : _sys(sys) { _element_data = std::make_unique<MooseVariableData>(); } void addSolution(int y) const { _element_data->addSolution(_sys.solution(), y); } private: std::unique_ptr<MooseVariableData> _element_data; // SystemBase & _sys; SystemBase _sys; }; int main() { SystemBase sys; MooseVariableFE var(sys); var.addSolution(7); }

Do you want me to remove the const deco for this method?

Yes, since it is ultimately modifying a data member (through a reference)

lindsayad · 2019-05-13T16:48:20Z

framework/src/variables/MooseVariableData.C

@@ -654,11 +655,326 @@ MooseVariableData<OutputType>::computeValues()
    computeAD(num_dofs, nqp);
 }

+template <>
+void
+MooseVariableData<RealArrayValue>::computeValues()


Sigh...so much code duplication here, but there's enough eigen specific logic that I don't have a better suggestion. If there were just a few different lines, I would suggest adding some shim methods, but there's too many.

I thought of creating a derived Eigen class for ourselves and implement few more methods to avoid this duplication. But not sure if it is better or not.

lindsayad · 2019-05-13T16:49:09Z

framework/src/variables/MooseVariableData.C

 template <typename OutputType>
 void
-MooseVariableData<OutputType>::setNodalValue(OutputType value, unsigned int idx /* = 0*/)
+MooseVariableData<OutputType>::setNodalValue(const OutputType & value, unsigned int idx)


lindsayad · 2019-05-13T16:52:34Z

framework/src/variables/MooseVariableData.C

+  auto n = _dof_indices.size();
+  libmesh_assert(n);
+
+  /* Idea of using Eigen::map from Derek Gaston (Yaqi: we can reevaluate this later) */


So you may use these lines of code in the future?

Might although I think we may not have it in a short term. @friedmud mentioned some multi-threading issue of this. So I am ok to remove this comment completely.

lindsayad · 2019-05-13T20:00:45Z

framework/include/base/Assembly.h

@@ -1159,26 +1359,51 @@ class Assembly

  void computeCurrentNeighborVolume();

+  /// Appling scaling, constraints to the local residual block and populate the full DoF indices


Our convention is to use:

/** * */

for methods and /// for data.

lindsayad · 2019-05-13T20:07:05Z

framework/doc/content/source/variables/ArrayMooseVariable.md

+The three rows correspond to standard, vector and array variables.
+OutputType is the data type used for templating.
+RealArrayValue is a typedef in [MooseTypes.h] as *Eigen::Matrix<Real, Eigen::Dynamic, 1>*.
+OutputShape is for the type of shape functions and OutputData is the type of basis function expansion coefficients that are stored in the solution vector.


One solution vector only holds one type (e.g. Real), so I'm a little unclear on what the last part of this sentence means. Can you explain a little more?

Oh, this should be the that are stored in the moose array variable, which are grabbed from the solution vector.

lindsayad · 2019-05-13T20:07:16Z

framework/doc/content/source/variables/ArrayMooseVariable.md

@@ -0,0 +1,89 @@
+# ArrayMooseVariable


I appreciate the strong documentation!

lindsayad · 2019-05-13T20:09:46Z

framework/doc/content/source/variables/ArrayMooseVariable.md

+Future work:
+
+- To profile the code with large number of variables to find the hot spots and fix them if there are any.
+- To change the current dof ordering for elemental variables so that we can avoid bunch of if statements with `isNodal()` in `MooseVariableFE.C`, (refer to [libMesh Issue 2114](https://github.com/libMesh/libmesh/issues/2114).


That would be fantastic. Can you open a ticket for this after this goes in?

lindsayad · 2019-05-13T20:14:55Z

framework/src/kernels/ArrayDiffusion.C

+  else if (_dc_type == 1)
+    _d_array = &getMaterialProperty<RealArrayValue>("diffusion_coefficient");
+  else if (_dc_type == 2)
+    _d_2d_array = &getMaterialProperty<RealArray>("diffusion_coefficient");


I looked at these typedefs in MooseTypes and it's not immediately clear to me the motivation behind their names. I look at RealArray and presume that it's an array of Reals. Is that right? With the enumeration description of full, I'm guessing that it's not. And what is RealArrayValue supposed to mean?

RealArrayValue is actually an Eigen vector type, equivalent with VectorXd, where X stands for dynamic, d stands for double. RealArray is actually an Eigen matrix, MatrixXd. I have to admit that I came up this name without too much thinking, so I am willing to accept a better name. Or I can change it to RealMatrixValue, what do you think?

Hmm...we're starting to run out of names that don't clash with existing names. I personally like RealEigenVector and RealEigenMatrix. @friedmud likes names; maybe he has a suggestion.

I do not quite Eigen in the name because we also used it for eigenvalue. It can cause confusion although we know here Eigen is the package name. The underlining data structures for vector and matrix in Eigen are the same. I do not have good ideas now for naming these two ;-) I will need the final for changing the code.

Since I do not have better ideas, I will make this change, I like these names of differentiating vector and matrix.

YaqiWang · 2019-05-15T07:09:38Z

@lindsayad comments are addressed.

…d idea, codes are added to ensure the entries in the vectors are unique

moosebuild · 2019-05-15T19:00:23Z

Job Precheck on 1ba3315 wanted to post the following:

Your code requires style changes.

A patch was auto generated and copied here
You can directly apply the patch by running, in the top level of your repository:

curl -s https://mooseframework.inl.gov/docs/PRs/13314/style.patch | git apply -v

Alternatively, with your repository up to date and in the top level of your repository:

git clang-format 2438af0e91976255bca7b3d25eab43569a7a9bfb

YaqiWang · 2019-05-16T16:26:56Z

@lindsayad do you have any further comments? Looks like our HPC is online. I'd like this to be merged so I can finish fixing apps before I go for a vacation next week. Thanks.

permcody · 2019-05-16T16:31:00Z

We aren't going to merge anything until we get next moving along. We definitely have a huge backup of PRs needing to go in, but with all of the red boxes I'm not inclined to pile on even more just yet. It very well could take us a week to get caught up again so I wouldn't worry about it!

…

On Thu, May 16, 2019 at 10:27 AM Yaqi ***@***.***> wrote: @lindsayad <https://github.com/lindsayad> do you have any further comments? Looks like our HPC is online. I'd like this to be merged so I can finish fixing apps before I go for a vacation next week. Thanks. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#13314?email_source=notifications&email_token=AAXFOIASSRY5XMNKAGQ5CZLPVWDNHA5CNFSM4HIQUPYKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODVSLDLA#issuecomment-493138348>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAXFOIDYJFSGLW3254CV6TTPVWDNHANCNFSM4HIQUPYA> .

lindsayad · 2019-05-16T18:01:36Z

I'm not going to approve anything while the documentation target is failing...

YaqiWang · 2019-05-16T18:41:17Z

@permcody said that he will take care of it. The problem is that I touched a tests file due to fixing a wrong error message. That tests file contains lots of tests which I don’t know how to write the missing SQA terms. Other than this all doc recipes passed.

…y variable

YaqiWang · 2019-05-17T06:46:40Z

The added newsletter item can be found at https://mooseframework.inl.gov/docs/PRs/13314/site/newsletter/2019_06.html.

permcody · 2019-05-21T22:24:46Z

Superseded by #13459

YaqiWang force-pushed the array_kernel branch 2 times, most recently from 35a37b6 to f211bee Compare April 27, 2019 04:30

YaqiWang changed the title ~~Array kernel~~ WIP: Array kernel Apr 30, 2019

YaqiWang force-pushed the array_kernel branch from 7084259 to 9ee76d5 Compare April 30, 2019 22:43

YaqiWang force-pushed the array_kernel branch from 9ee76d5 to 3edba17 Compare May 5, 2019 00:13

YaqiWang changed the title ~~WIP: Array kernel~~ Array kernel May 6, 2019

YaqiWang force-pushed the array_kernel branch from 8b42313 to e1b69b9 Compare May 6, 2019 17:19

YaqiWang force-pushed the array_kernel branch from 49e27ef to bd344f7 Compare May 7, 2019 16:08

YaqiWang force-pushed the array_kernel branch from 34a29fd to 8885805 Compare May 8, 2019 20:24

YaqiWang force-pushed the array_kernel branch 4 times, most recently from 2a41873 to 223af90 Compare May 9, 2019 21:32

YaqiWang added 2 commits May 10, 2019 01:25

a compilable stage with variable, coupleable, assembly and systems ex…

88e4da6

…tended to array variable idaholab#6881

add array variable action and array initial condition and a simple test

adfa322

idaholab#6881

YaqiWang added 4 commits May 10, 2019 03:25

add capability to allow couple array variables in a standard kernel

4c512c5

documentations of array kernel and SQA requirements

9a8c164

final clang formatting after resolving conflict

5aaa6c8

add few neighbor coupleable functions that will be useful for array i…

cf2127c

…nterface kernels or constraints

YaqiWang force-pushed the array_kernel branch from 223af90 to cf2127c Compare May 10, 2019 09:49

add a missing function implementation and some cleanups

e92ff27

YaqiWang force-pushed the array_kernel branch from bce444d to e92ff27 Compare May 10, 2019 22:17

lindsayad requested changes May 13, 2019

View reviewed changes

YaqiWang force-pushed the array_kernel branch from 4bf04c2 to c151c65 Compare May 15, 2019 07:14

YaqiWang added 3 commits May 15, 2019 01:21

address Alex's comments

399c135

change back to vector, using set to store points could be a really ba…

2af76d8

…d idea, codes are added to ensure the entries in the vectors are unique

rename RealArrayValue to RealEigenVector

2663461

YaqiWang force-pushed the array_kernel branch from c151c65 to 1ba3315 Compare May 15, 2019 07:34

YaqiWang mentioned this pull request May 15, 2019

Local DoF ordering for elemental variables #13417

Open

rename RealArray to RealEigenMatrix

c73516f

YaqiWang force-pushed the array_kernel branch from 1ba3315 to c73516f Compare May 15, 2019 21:42

YaqiWang added 2 commits May 16, 2019 23:53

add array kernel in the coming newsletter

f5267d8

add a moose error to prevent users from using vector types in an arra…

8d2e314

…y variable

permcody mentioned this pull request May 21, 2019

Array kernel #13459

Closed

permcody closed this May 21, 2019

YaqiWang mentioned this pull request Jun 10, 2019

Array kernel #13528

Merged

		@@ -1159,26 +1359,51 @@ class Assembly

		void computeCurrentNeighborVolume();

		/// Appling scaling, constraints to the local residual block and populate the full DoF indices

Array kernel #13314

Array kernel #13314

Conversation

YaqiWang commented Apr 25, 2019 • edited

YaqiWang commented Apr 25, 2019

YaqiWang commented May 1, 2019

YaqiWang commented May 6, 2019

moosebuild commented May 6, 2019 • edited

YaqiWang commented May 6, 2019 • edited

YaqiWang commented May 6, 2019

permcody commented May 6, 2019

YaqiWang commented May 6, 2019 • edited

YaqiWang commented May 7, 2019

lindsayad commented May 7, 2019

YaqiWang commented May 7, 2019

YaqiWang commented May 7, 2019

YaqiWang commented May 8, 2019 • edited

YaqiWang commented May 8, 2019

YaqiWang commented May 8, 2019

YaqiWang commented May 10, 2019 • edited

lindsayad left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

YaqiWang commented May 15, 2019

moosebuild commented May 15, 2019

YaqiWang commented May 16, 2019

permcody commented May 16, 2019 via email

lindsayad commented May 16, 2019

YaqiWang commented May 16, 2019

YaqiWang commented May 17, 2019

permcody commented May 21, 2019

YaqiWang commented Apr 25, 2019 •

edited

moosebuild commented May 6, 2019 •

edited

YaqiWang commented May 6, 2019 •

edited

YaqiWang commented May 6, 2019 •

edited

YaqiWang commented May 8, 2019 •

edited

YaqiWang commented May 10, 2019 •

edited