Linear elastic FFT solver #441

recuero · 2020-08-03T19:05:35Z

This PR aims at solving a periodic small deformation problem with multiple phases using a spectral approach.

Initial comparison with FE solution
Verify for various global strains
Establish a convergence criterion and include it in the test
Handle graciously the various strain modes, give option to user
Let the user choose a "homogeneized" material ratio for the solver

Refs. #401

To make this run, explicit instantiation of auto RankFourTensorTempl<T>::operator*(const Tensor<T2> & b) const is required in MOOSE's RankFourTensor.h

…verse strains

…l strain

moosebuild · 2020-08-07T21:31:52Z

Job Precheck on 786951b 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/magpie/docs/PRs/441/style.patch | git apply -v

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

git clang-format c19bd79fd44acf91a4276b80db1d1275daa319ee

moosebuild · 2020-08-07T22:18:52Z

Job Test on 4383ed8 wanted to post the following:

View the site here

This comment will be updated on new commits.

dschwen · 2020-08-10T16:07:51Z

examples/spectralLinearElastic2MaterialsTM.i

+    youngs_modulus = 1.0e4
+    poissons_ratio = 0.3
+    elasticity_tensor_prefactor = young_modulus_function
+    # Is the line above during the righ thing by scaling the young modulus?


Suggested change

# Is the line above during the righ thing by scaling the young modulus?

# Is the line above during the right thing by scaling the young modulus?

That comment shouldn't be there anymore...

include/executioners/SpectralExecutionerLinearElastic.h

src/executioners/SpectralExecutionerLinearElastic.C

dschwen · 2020-08-10T16:18:32Z

src/executioners/SpectralExecutionerLinearElastic.C

+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        epsilon_buffer.realSpace()[index] = _initial_strain_tensor;
+        index++;
+      }
+}


how about we add

FFTData<T> & operator=(const T & rhs);

which internally uses std::fill to set the data. This will be more compact and easier to optimize for the compiler.

The code in this file would then just be

epsilon_buffer.realSpace() = _initial_strain_tensor;

Makes it much cleaner. The meaning of the operation becomes more implicit at first sight though.

dschwen · 2020-08-10T16:20:04Z

src/executioners/SpectralExecutionerLinearElastic.C

+
+  const int ni = grid[0];
+  const int nj = grid[1];
+  const int nk = (grid[2] >> 1) + 1;


should probably be gamma_hat.kTable(2).size(), right?

Nope. We are mindlessly resizing and filling out the last dimension of _k_table (should we change this?). In addition, it should not make a difference for this problem as all variables and parameters are defined on the same grid.

// precompute kvector components for current direction _k_table[i].resize(_grid[i]); for (int j = 0; j < _grid[i]; ++j) _k_table[i][j] = 2.0 * libMesh::pi * ((j < (_grid[i] >> 1) + 1) ? Real(j) : (Real(j) - _grid[i])) / _box_size(i);

I guess I am changing it to address other comments

src/executioners/SpectralExecutionerLinearElastic.C

dschwen · 2020-08-10T16:25:21Z

src/executioners/SpectralExecutionerLinearElastic.C

+                Complex q_square = freq[0] * freq[0] + freq[1] * freq[1] + freq[2] * freq[2];
+                if (std::abs(q_square) > 1.0e-12)
+                {
+                  gamma_hat.reciprocalSpace()[index](i, j, k, l) =


it might be best to grab a reference to gamma_hat.reciprocalSpace() outside of the loop.

auto & gamma_hat_F = gamma_hat.reciprocalSpace()

(or whatever naming convention you want to adopt for recoprocal fields)

Would just be saving the function call, right?

A function call in a tight inner loop. As that function is implemented in the header file there is a very good chance the compiler will inline it, however we might as well be explicit about it.

dschwen · 2020-08-10T16:32:52Z

src/executioners/SpectralExecutionerLinearElastic.C

+  for (int freq_x = 0; freq_x < ni; ++freq_x)
+    for (int freq_y = 0; freq_y < nj; ++freq_y)
+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        stress_buffer.realSpace()[index] =
+            elastic_tensor.realSpace()[index] * epsilon_buffer.realSpace()[index];
+        index++;
+      }


We should have a templated helper function to do these types multiplications in FFTData

template <typename T> template <typename T1, typename T2> void FFTData<T>::setToProduct(const FFTData<T1> & m1, const FFTData<T2> & m2) { // assert that grid sizes are identical // loop over indices and multiply }

Grid size information is in FFTBuffer. In FFTData you have either a real or a complex buffer. I am adding a version of your suggestion.

dschwen · 2020-08-10T16:54:27Z

src/executioners/SpectralExecutionerLinearElastic.C

+          epsilon_buffer.reciprocalSpace()[index] = _initial_strain_tensor;
+        else
+          epsilon_buffer.reciprocalSpace()[index] -=
+              gamma_hat.reciprocalSpace()[index] * stress_buffer.reciprocalSpace()[index];


For this loop we could provide a helper function that takes a Lambda as an argument

template<typename T> template <typename T1> FFTData<T>::apply(T1 lambda) { size_t index = 0; for (auto ivec : kTable(0)) for (auto jvec : kTable(1)) for (auto kvec : kTable(2)) _buffer[index++] = lambda(ivec, jvec, kvec, index); }

(unfortunately FFTData does not have access to the kTable yet, but we can probably store a referece to the FFTBufferBase object that owns the FFTData.)

so the call would look something like

epsilon_buffer.reciprocalSpace().apply([](Real ivec, Real jvec, Real kvec, std::size_t index) { if (ivec == 0 && jvec == 0 && kvec == 0) return _initial_strain_tensor; else return gamma_hat.reciprocalSpace()[index] * stress_buffer.reciprocalSpace()[index]; } );

Since we are using templating to pass in the lambda the compiler will inline it and should produce very fast code.

Also you ultimately will want to do

epsilon_buffer.reciprocalSpace().apply([](Real ivec, Real jvec, Real kvec, std::size_t index) { return gamma_hat.reciprocalSpace()[index] * stress_buffer.reciprocalSpace()[index]; } ); epsilon_buffer.setKZero(_initial_strain_tensor);

where setKZero needs to be implemented (and will simply write to the index that corresponds to the (0,0,0) k-vector. That way you'll write twice to that location but avoid looping over the conditional. This is likely much faster.

dschwen · 2020-08-10T16:55:38Z

src/executioners/SpectralExecutionerLinearElastic.C

+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        stress_buffer.realSpace()[index] =
+            (elastic_tensor.realSpace()[index]) * epsilon_buffer.realSpace()[index];


see above for the setToProduct function

dschwen · 2020-08-10T16:56:35Z

src/executioners/SpectralExecutionerLinearElastic.C

+      for (int freq_z = 0; freq_z < nk; ++freq_z)
+      {
+        // Define elastic tensor
+        std::vector<Real> iso_const(2);


avoid creating short lived temporary vectors like the plague. They entail heap allocations and are thus very slow.

I.e. move it out of the loop and since the size is known at compile time we could also think about using arrays (and templating the fillFormVector routines).

dschwen · 2020-08-10T17:16:01Z

src/executioners/SpectralExecutionerLinearElastic.C

+                            stress.reciprocalSpace()[index](2, 2) * freq[2];
+
+        error_n += kvector_stress[0] * kvector_stress[0] + kvector_stress[1] * kvector_stress[1] +
+                   kvector_stress[1] * kvector_stress[1];


is this a typo? I'd expect some [2] here

dschwen · 2020-08-10T17:17:45Z

src/executioners/SpectralExecutionerLinearElastic.C

+      {
+        const std::array<Complex, 3> freq{ivec[freq_x] * I, jvec[freq_y] * I, kvec[freq_z] * I};
+
+        std::array<Complex, 3> kvector_stress;


If you used ComplexVectorValue this would be a lot more compact.

Mostly delegate some buffer operations to FFTData<T>. Use of lambdas.

include/utils/FFTData.h

dschwen · 2020-08-11T22:01:51Z

src/executioners/SpectralExecutionerLinearElastic.C

+        ComplexVectorValue kvector_stress;
+        kvector_stress(0) = stress.reciprocalSpace()[index](0, 0) * freq(0) +
+                            stress.reciprocalSpace()[index](1, 0) * freq(1) +
+                            stress.reciprocalSpace()[index](2, 0) * freq(2);
+        kvector_stress(1) = stress.reciprocalSpace()[index](0, 1) * freq(0) +
+                            stress.reciprocalSpace()[index](1, 1) * freq(1) +
+                            stress.reciprocalSpace()[index](2, 1) * freq(2);
+        kvector_stress(2) = stress.reciprocalSpace()[index](0, 2) * freq(0) +
+                            stress.reciprocalSpace()[index](1, 2) * freq(1) +
+                            stress.reciprocalSpace()[index](2, 2) * freq(2);


isn't that just

kvector_stress = stress.reciprocalSpace()[index] * freq;

??

dschwen · 2020-08-11T22:02:20Z

src/executioners/SpectralExecutionerLinearElastic.C

+        error_n += kvector_stress(0) * kvector_stress(0) + kvector_stress(1) * kvector_stress(1) +
+                   kvector_stress(2) * kvector_stress(2);


Suggested change

error_n += kvector_stress(0) * kvector_stress(0) + kvector_stress(1) * kvector_stress(1) +

kvector_stress(2) * kvector_stress(2);

error_n += kvector_stress.norm_sq();

.norm_sq() takes the norm of each component's complex number, not what we look for here.

Oh, it would be the component times its complex conjugate rather than component square? Thats really not what is wanted here?

norm_sq obtains the norm of each of the components and sums them up. The error formula I am basing this off takes the norm of the entire vector. Results are different.

There is more work to do in terms of convergence checks. The formula implemented here captures stress equilibrium. But we could add BC error and others. We can leave all this for next revision.

dschwen · 2020-08-11T22:03:56Z

src/executioners/SpectralExecutionerLinearElastic.C

-          for (int i = 0; i < ndim; i++)
-            for (int j = 0; j < ndim; j++)
-              error_0 += stress.reciprocalSpace()[0](i, j) * stress.reciprocalSpace()[0](i, j);
+          error_0 = stress.reciprocalSpace()[0].L2norm() * stress.reciprocalSpace()[0].L2norm();


not to self: add RankTwoTensorTempl<T>::L2norm_sq() analogous to norm_sq()

dschwen · 2020-08-11T22:05:09Z

src/utils/FFTData.C

+{
+  mooseAssert(grid.size() == 3, "Product defined for 3 dimensions");
+
+  for (unsigned int index = 0; index < grid[0] * grid[1] * grid[2]; index++)


passing in grid seems unnecessary if you can just evaluate _buffer.size()

Initial stab at linear elastic FFT solver Refs. idaholab#401

c0a274f

To make this run, explicit instantiation of auto RankFourTensorTempl<T>::operator*(const Tensor<T2> & b) const is required in MOOSE's RankFourTensor.h

recuero added the Do Not Merge label Aug 3, 2020

recuero added 6 commits August 6, 2020 17:04

Work in progress

77daaac

Converges to reasonable quantities

bf4c2f4

Minor changes. Does not compare well to tensor mechanics solution: in…

4705799

…verse strains

Clean up

aad41f1

Add a test on the linear elastic spectral executioner

39b76b9

More clean up

e3ba0a7

recuero self-assigned this Aug 6, 2020

recuero mentioned this pull request Aug 7, 2020

Add explicit instantiation for RankFourTensor complex operation idaholab/moose#15734

Closed

recuero added 3 commits August 7, 2020 13:40

Add tensor mechanics example for verification

9e9a200

Update example and test to reflect strain choice and solver tolerance

09e2838

Add convergence check and provide options to impose any type of globa…

e6d8c47

…l strain

recuero force-pushed the smallElasticFFT branch 2 times, most recently from 4e68f33 to 786951b Compare August 7, 2020 21:31

Add explicit instantiation of operator* for RankFourTensor

c799351

recuero force-pushed the smallElasticFFT branch from 786951b to c799351 Compare August 7, 2020 21:32

Some cleanup

4383ed8

recuero marked this pull request as ready for review August 7, 2020 22:44

recuero added ready for review and removed Do Not Merge labels Aug 7, 2020

recuero requested a review from dschwen August 7, 2020 22:45

Add const qualifiers

bb5ecbb

dschwen reviewed Aug 10, 2020

View reviewed changes

include/executioners/SpectralExecutionerLinearElastic.h Outdated Show resolved Hide resolved

dschwen reviewed Aug 10, 2020

View reviewed changes

src/executioners/SpectralExecutionerLinearElastic.C Outdated Show resolved Hide resolved

dschwen reviewed Aug 10, 2020

View reviewed changes

src/executioners/SpectralExecutionerLinearElastic.C Outdated Show resolved Hide resolved

dschwen reviewed Aug 10, 2020

View reviewed changes

src/executioners/SpectralExecutionerLinearElastic.C Outdated Show resolved Hide resolved

dschwen reviewed Aug 10, 2020

View reviewed changes

src/executioners/SpectralExecutionerLinearElastic.C Outdated Show resolved Hide resolved

dschwen reviewed Aug 10, 2020

View reviewed changes

src/executioners/SpectralExecutionerLinearElastic.C Outdated Show resolved Hide resolved

dschwen reviewed Aug 10, 2020

View reviewed changes

Adress reviewer's comments.

c49c47f

Mostly delegate some buffer operations to FFTData<T>. Use of lambdas.

dschwen reviewed Aug 11, 2020

View reviewed changes

include/utils/FFTData.h Show resolved Hide resolved

dschwen reviewed Aug 11, 2020

View reviewed changes

More changes

5e602e8

dschwen self-assigned this Aug 24, 2020

dschwen approved these changes Oct 6, 2020

View reviewed changes

dschwen merged commit 1c702b2 into idaholab:devel Oct 6, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Linear elastic FFT solver #441

Linear elastic FFT solver #441

recuero commented Aug 3, 2020 •

edited

Loading

moosebuild commented Aug 7, 2020

moosebuild commented Aug 7, 2020 •

edited

Loading

dschwen Aug 10, 2020

dschwen Aug 10, 2020

recuero Aug 10, 2020

dschwen Aug 10, 2020

dschwen Aug 10, 2020

recuero Aug 10, 2020

dschwen Aug 10, 2020

recuero Aug 11, 2020

recuero Aug 11, 2020

dschwen Aug 10, 2020

recuero Aug 11, 2020

dschwen Aug 11, 2020

dschwen Aug 10, 2020 •

edited

Loading

recuero Aug 11, 2020

dschwen Aug 10, 2020

dschwen Aug 10, 2020 •

edited

Loading

dschwen Aug 10, 2020

dschwen Aug 10, 2020

dschwen Aug 10, 2020

dschwen Aug 10, 2020

dschwen Aug 11, 2020

dschwen Aug 11, 2020

recuero Aug 11, 2020

dschwen Aug 11, 2020

recuero Aug 12, 2020 •

edited

Loading

dschwen Aug 11, 2020

dschwen Aug 11, 2020

	# Is the line above during the righ thing by scaling the young modulus?
	# Is the line above during the right thing by scaling the young modulus?

		error_n += kvector_stress(0) * kvector_stress(0) + kvector_stress(1) * kvector_stress(1) +
		kvector_stress(2) * kvector_stress(2);

	error_n += kvector_stress(0) * kvector_stress(0) + kvector_stress(1) * kvector_stress(1) +
	kvector_stress(2) * kvector_stress(2);
	error_n += kvector_stress.norm_sq();

Linear elastic FFT solver #441

Linear elastic FFT solver #441

Conversation

recuero commented Aug 3, 2020 • edited Loading

moosebuild commented Aug 7, 2020

moosebuild commented Aug 7, 2020 • edited Loading

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

dschwen Aug 10, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dschwen Aug 10, 2020 • edited Loading

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

recuero Aug 12, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

recuero commented Aug 3, 2020 •

edited

Loading

moosebuild commented Aug 7, 2020 •

edited

Loading

dschwen Aug 10, 2020 •

edited

Loading

dschwen Aug 10, 2020 •

edited

Loading

recuero Aug 12, 2020 •

edited

Loading