added test to gradient calculation logic described in #149 #150

bodduv · 2017-06-28T15:54:34Z

Edit: I will clean up the QC class with some updated documentation, with several test cases in the subsequent PRs.

davydden

i like it, but there is more to do

davydden · 2017-06-28T16:07:37Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+
+
+
+// Short test to check MoleculeHandler::get_neighbor_lists() function.


is description correct?

davydden · 2017-06-28T16:10:03Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+    auto my_pair_2 = GridTools::find_active_cell_around_point(MappingQ1<dim>(),
+                                                              dof_handler,
+                                                              p2);
+    std::vector<Point<dim>>points = {my_pair_1.second, my_pair_2.second};


space after >. p.s. did you run indentation on this?

davydden · 2017-06-28T16:11:35Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+                 ExcInternalError());
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    LinearAlgebraDealII::Vector local_gradient_I, local_gradient_II;


this should just be dealii::Vector<double>, it's not a parallel vector.

davydden · 2017-06-28T16:12:55Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+    std::vector<dealii::types::global_dof_index>
+    local_dof_indices(dofs_per_cell);
+
+    FEValues<dim>fe_values(mapping, fe, Quadrature<dim>(points), update_values);


you need to have separate FEValues for each cell. Currently you create a quadrature which has positions of both points inside each cell.

Also add space after >.

davydden · 2017-06-28T16:15:35Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+
+    for (unsigned int k = 0; k < dofs_per_cell; ++k)
+      {
+        local_gradient_I [k] += fe_values.shape_value(k, 0);


you don't do the correct gradient, it should be

n^{ab}*N_i(X^a)

and you need extractors.

I thought of doing just one dof at a mesh node first.

I thought of doing just one dof

but that's not what i mean and you don't do one dof. You calculate force due to this pair on all dofs currently (for k=0; k < dofs_per_cell; ++k)).

What I am saying is that as far as I can tell this line is wrong.

davydden · 2017-06-28T16:16:44Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+    constraints.distribute_local_to_global(local_gradient_II,
+                                           local_dof_indices,
+                                           gradient);
+


at this point we have answer from our favourite approach, now we need to compare it to what we are doing at the moment. Otherwise we don't know if it actually works. The fact that the test compiles and gives you some numbers does not mean yet that it is working 😄

I compared the values by manual calculation for each dof.

but you did not do scalar product n^{ab} * N_i(X^a). i wonder how/why does it work in its current form

I was just comparing the difference of the shape functions.

if you do it correctly, you should probably get the same result as for your cases below n^{ab}={1,0,0}.

davydden · 2017-06-28T16:33:37Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+                                                              p2);
+    std::vector<Point<dim>>points = {my_pair_1.second, my_pair_2.second};
+
+    run(points);


may i suggest to rework it into:

run (const Point<dim> &X1, auto &my_pair_1, const Point<dim> &X2, auto &my_pair_2)

and correct internals accordingly.

Then you should get the same result from your test regardless of the order of points, i.e.

Point<3> p1(0.8,0.8,0.8),p2(1.8,0.8,0.8); problem.check(p1,p2);

vs

Point<3> p2(0.8,0.8,0.8),p1(1.8,0.8,0.8);
problem.check(p1,p2);

bodduv · 2017-06-28T16:37:20Z

I will work on the test again. It appears that I only compared manually the first component. Regarding the older method we didn't finish it to compare to.

Edit: Should I finish the old method and then compare the both?

davydden · 2017-06-28T16:40:22Z

I will work on the test again. It appears that I only compared manually the first component.

May I suggest you use Maxima to produce MWE of the standard-QC (blessed output) and choose only 2 points which, however, have Xa-Xb not parallel to any of the axis?

Then you only need to test (pa,pb) vs (pb,pa) vs blessed to make sure it's correct. I don't think you need 3 different pairs to actually test it.

Regarding the older method we didn't finish it to compare to.

that's ok if you have MWE in Maxima or alike.

davydden · 2017-06-28T16:48:29Z

I suggest to do 2D as it's easier to test properly in Maxima and also easier to visualize

you can do vectors in Maxima easily:

xa_x: [-R*0.5, R];

depending on how you implement this, you could introduce DoFs for each element

DOFS: [[1,2,...],[2,3,...]];

you can get actual DoFs on the mesh in deal.II from extra output, i.e. https://github.com/dealii/dealii/blob/master/tests/dofs/dof_tools_extract_dofs_with_support_contained_within_01.cc#L90-L113

Here are some examples of how to do for loops in maxima:

N : [1-x,x];
H : [1,1];
nelem : 2;
u1: sin(W*t);
T1: c*t;
nlocaldofs : 2;
DOFS: [[1,2],[2,3]];
M : zeromatrix (3, 3);
K : zeromatrix (3, 3);
R : [0,0,0];
for el:1 thru nelem do
   for i:1 thru nlocaldofs do
     for j:1 thru nlocaldofs do
     (
       ii : DOFS[el][i],
       jj : DOFS[el][j],
       a  : H[el]*integrate(N[i]*N[j],x,0,1),
       b  : (1/H[el])*integrate(diff(N[i],x)*diff(N[j],x),x,0,1),
       M[ii][jj] : M[ii][jj] + a,
       K[ii][jj] : K[ii][jj] + b       
     );
for el:1 thru nelem do
   for i:1 thru nlocaldofs do
     (
       ii : DOFS[el][i],
       a  : H[el]*integrate(N[i],x,0,1),
       R[ii] : R[ii] + a       
     );

Of course, the implementation should be along the standard QC approach where you can evaluate each shape function on any element.

p.s. @vishalkenchan i edited this post to hopefully give enough Maxima syntax to implement the test.

bodduv · 2017-06-28T17:05:48Z

Thanks! I will have a look at it tomorrow.

davydden · 2017-06-29T09:40:33Z

thinking about this MWE, i would probably proceeded by defining global shape functions for two elements as scalar-valued functions (so that you can evaluate it on any point x,y).

Then, since the vector-valued shape functions in our case are primitive (i.e. F0=[f0,0]; F1=[0,f0], then the scalar products you are after to get forces is essentially: take either first or second component in n^{ab} and multiply with the difference of scalar valued shape functions.

At the end of the day you would have in Maxima x,y components of forces on each node (globally) of 2 element mesh. Knowing how deal.II enumerates DoFs (see above), you would be able to compare the results exactly.

davydden · 2017-06-29T09:54:27Z

here's how you can define such shape functions and plot them to make sure it's all good:

N1(x,y):= if x < 1 then (x*y) else (2-x)*y;
plot3d (N1(x,y), [x, 0, 2], [y, 0, 1], [grid, 10, 10],
[mesh_lines_color,false])$

bodduv · 2017-06-29T10:00:50Z

Thanks for the code snippets. This makes a sound way to write test.

bodduv · 2017-06-29T11:43:20Z

If I am using just FESystem<dim> fe (FE_Q<dim>(1),1);, this would give me scalar valued shape functions right? I could then directly compare the results from maxima. I would be comparing 6 x 3 doubles for 2D case. If I do this, am I still testing what we wanted to test?

davydden · 2017-06-29T11:51:54Z

Not really! In dealii u need to use vector-valued function as we use normally. In maxima u will have 2x6 values for forces, in deal.ii the vector will be of size 12. You can easily relate the two if u plot DoFs with the mesh (see above).

davydden · 2017-06-29T15:20:10Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+    for (unsigned int i = 0; i < dof_handler.n_dofs(); i+=dim)
+      {
+        for (int d = 0; d < dim; ++d)
+          std::cout << n[d]*gradient[i] <<  "\t";


you use vector-valued functions, why don't you just output the gradient itself? I was expecting that you use actual positions of the two points in your run().

see below on how i would rework it.

I was trying get to the difference of shape function values and comparing before multiplying with the unit normal. But I could do this as well.

I am just in favour to do it as close as the actual implementation so that you could just copy-paste this code and add partial_differentive_of_vector_wrt_dofs_02.cc which does exactly the same but using our QC class (with whatever potential and unit cluster weights).

Sure. I will change it.

davydden · 2017-06-29T15:20:34Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+      //problem.write_dofs_and_support_points_info();
+      {
+        Point<2> p1(0.23, 0.37), p2(1.73, 0.43);
+        problem.check(p1,p2);


please check

problem.check(p2,p1);

davydden · 2017-06-29T15:21:15Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+                                                              dof_handler,
+                                                              p2);
+
+    run (my_pair_1.second, my_pair_2.second);


in order not to assume which point is in which cell, i would provide here my_pair_1, my_pair_2, p1, p2

davydden · 2017-06-29T15:23:32Z

tests/core/partial_differentive_of_vector_wrt_dofs_01.cc

+
+      for (unsigned int k = 0; k < dofs_per_cell; ++k)
+        {
+          local_gradient_I [k]  =  fe_value_a.shape_value(k, 0);


Please use ~~extractors and~~ scalar products. Otherwise you may not be doing what you think you are doing.

p.s. i don't think this would change the result as shape_value returns

If the shape function is vector-valued, then this returns the only non- zero component

for primitive functions. So i guess it's fine to stick with it, but then at least get which component this corresponds to implement scalar product.

That's also how we should implement it in the code in order to avoid multiplications with zeroes.

to get what is the only non-zero component for a given shape functions, you would need fe.system_to_component_index(k).first, see dealii/dealii#3646 (comment)

Since in our case we have primitive shape functions, could I just do what we are doing in qc.cc:

for (unsigned int k = 0; k < dofs_per_cell; ++k) { local_gradient_I [k] = fe_value_a[u_fe].value(k, 0); local_gradient_II [k] = -fe_value_b[u_fe].value(k, 0); }

With scalar product.

no, you misunderstand my point! You need to do

for (unsigned int k = 0; k < dofs_per_cell; ++k) { const unsigned int comp = fe.system_to_component_index(k).first; local_gradient_I [k] = nab[comp] * fe_value_a.shape_value(k, 0); }

that would be the actual (TODO) implementation in QC class. Otherwise we will do 2 out of 3 multiplications with zeroes (in 3D case).

davydden

nice, minor changes to make it closer to the actual (TODO) implementation.

davydden

nice, i guess we are now fully convinced that this approach works! 😄

bodduv · 2017-06-29T18:28:20Z

👍

davydden · 2017-06-29T18:30:48Z

@vishalkenchan i would now add another test (together with MWE in Maxima) and use this logic in QC to calculate forces from two cluster atoms (say spring potential) for the same mesh.

The MWE in Maxima is the same up to dE/dr * (n1+n2)/2 multiplier which you can hard-code according to the chosen potential E and cluster weights n1 and n2.

bodduv · 2017-06-29T23:28:34Z

Are you going to touch compute_energy_gradient()? In this case, I will work on something else (probably address some of the issues in #118).

davydden · 2017-06-30T05:37:40Z

I meant that it's the next logical step for u to do ;-)

davydden · 2017-06-30T10:06:00Z

@vishalkenchan i will look at doing forces in QC, but not necessarily test it ;-)

davydden requested changes Jun 28, 2017

View reviewed changes

davydden reviewed Jun 28, 2017

View reviewed changes

bodduv force-pushed the fix/gradient_monoatomic branch from fb3bbdc to 8ae9fd5 Compare June 29, 2017 15:13

bodduv requested a review from davydden June 29, 2017 15:13

davydden reviewed Jun 29, 2017

View reviewed changes

davydden requested changes Jun 29, 2017

View reviewed changes

bodduv force-pushed the fix/gradient_monoatomic branch from 8ae9fd5 to 77f5824 Compare June 29, 2017 18:09

added test to gradient calculation logic described in #149

e6c95c0

bodduv force-pushed the fix/gradient_monoatomic branch from 77f5824 to e6c95c0 Compare June 29, 2017 18:15

bodduv requested a review from davydden June 29, 2017 18:16

davydden approved these changes Jun 29, 2017

View reviewed changes

davydden merged commit 3f7990f into ltm-erlangen:master Jun 29, 2017

bodduv deleted the fix/gradient_monoatomic branch June 29, 2017 18:28

davydden mentioned this pull request Jun 29, 2017

Gradient calculation without hanging nodes #149

Closed




		// Short test to check MoleculeHandler::get_neighbor_lists() function.

added test to gradient calculation logic described in #149 #150

added test to gradient calculation logic described in #149 #150

Conversation

bodduv commented Jun 28, 2017 • edited by davydden

davydden 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

davydden Jun 28, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden Jun 28, 2017 • edited

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

bodduv commented Jun 28, 2017 • edited

davydden commented Jun 28, 2017 • edited

davydden commented Jun 28, 2017 • edited

bodduv commented Jun 28, 2017

davydden commented Jun 29, 2017 • edited

davydden commented Jun 29, 2017 • edited

bodduv commented Jun 29, 2017

bodduv commented Jun 29, 2017 • edited

davydden commented Jun 29, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden Jun 29, 2017 • edited

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

davydden Jun 29, 2017 • edited

Choose a reason for hiding this comment

davydden Jun 29, 2017 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

davydden Jun 29, 2017 • edited

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

davydden Jun 29, 2017 • edited

Choose a reason for hiding this comment

davydden left a comment

Choose a reason for hiding this comment

davydden left a comment

Choose a reason for hiding this comment

bodduv commented Jun 29, 2017

davydden commented Jun 29, 2017 • edited

bodduv commented Jun 29, 2017

davydden commented Jun 30, 2017

davydden commented Jun 30, 2017

bodduv commented Jun 28, 2017 •

edited by davydden

davydden Jun 28, 2017 •

edited

davydden Jun 28, 2017 •

edited

bodduv commented Jun 28, 2017 •

edited

davydden commented Jun 28, 2017 •

edited

davydden commented Jun 28, 2017 •

edited

davydden commented Jun 29, 2017 •

edited

davydden commented Jun 29, 2017 •

edited

bodduv commented Jun 29, 2017 •

edited

davydden commented Jun 29, 2017 •

edited

davydden Jun 29, 2017 •

edited

davydden Jun 29, 2017 •

edited

davydden Jun 29, 2017 •

edited

davydden Jun 29, 2017 •

edited

davydden Jun 29, 2017 •

edited

davydden commented Jun 29, 2017 •

edited