Add FEM based gradient reference (idaholab#401)

examples/spectral.i

@@ -12,7 +12,7 @@
 []
 
 [AuxVariables]
-  [./c_aux]
+  [./c_fft]
     order = CONSTANT
     family = MONOMIAL
     [./InitialCondition]
@@ -21,23 +21,23 @@
     [../]
   [../]
 
-  [./R0_aux]
+  [./R0_fft]
     order = CONSTANT
     family = MONOMIAL
     [./InitialCondition]
       type = FunctionIC
       function = 'cos(x/100*2*pi*2)*cos(y/100*2*pi*4)'
     [../]
   [../]
-  [./R1_aux]
+  [./R1_fft]
     order = CONSTANT
     family = MONOMIAL
     [./InitialCondition]
       type = FunctionIC
       function = 'cos(x/100*2*pi*3)*cos(y/100*2*pi*3)'
     [../]
   [../]
-  [./R2_aux]
+  [./R2_fft]
     order = CONSTANT
     family = MONOMIAL
     [./InitialCondition]
@@ -47,7 +47,7 @@
   [../]
 
 
-  [./u_aux]
+  [./u_fft]
     order = CONSTANT
     family = MONOMIAL
     [./InitialCondition]
@@ -60,11 +60,31 @@
       outvalue = 0
     [../]
   [../]
-  [./grad_u0_aux]
+  [./grad_u0_fft]
     order = CONSTANT
     family = MONOMIAL
   [../]
-  [./grad_u1_aux]
+  [./grad_u1_fft]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  
+  [./u_fem]
+    [./InitialCondition]
+      type = SmoothCircleIC
+      x1 = 50
+      y1 = 50
+      radius = 30
+      int_width = 20
+      invalue = 1
+      outvalue = 0
+    [../]
+  [../]
+  [./grad_u0_fem]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
+  [./grad_u1_fem]
     order = CONSTANT
     family = MONOMIAL
   [../]
@@ -82,73 +102,88 @@
   # Buffers
   [./c]
     type = RealFFTWBuffer
-    moose_variable = c_aux
+    moose_variable = c_fft
   [../]
   [./R]
     type = RealVectorValueFFTWBuffer
-    moose_variable = 'R0_aux R1_aux R2_aux'
+    moose_variable = 'R0_fft R1_fft R2_fft'
   [../]
 
   [./u]
     type = RealFFTWBuffer
-    moose_variable = u_aux
+    moose_variable = u_fft
   [../]
   [./grad_u]
     type = RealVectorValueFFTWBuffer
   [../]
 []
 
 [AuxKernels]
-  [./c_aux]
+  [./c_fft]
     type = FFTBufferAux
-    variable = c_aux
+    variable = c_fft
     fft_buffer = c
     execute_on = FINAL
   [../]
 
-  [./R0_aux]
+  [./R0_fft]
     type = FFTBufferAux
-    variable = R0_aux
+    variable = R0_fft
     fft_buffer = R
     component = 0
     execute_on = FINAL
   [../]
-  [./R1_aux]
+  [./R1_fft]
     type = FFTBufferAux
-    variable = R1_aux
+    variable = R1_fft
     fft_buffer = R
     component = 1
     execute_on = FINAL
   [../]
-  [./R2_aux]
+  [./R2_fft]
     type = FFTBufferAux
-    variable = R2_aux
+    variable = R2_fft
     fft_buffer = R
     component = 2
     execute_on = FINAL
   [../]
 
-  [./u_aux]
+  [./u_fft]
     type = FFTBufferAux
-    variable = u_aux
+    variable = u_fft
     fft_buffer = u
     execute_on = FINAL
   [../]
 
-  [./grad_u0_aux]
+  [./grad_u0_fft]
     type = FFTBufferAux
-    variable = grad_u0_aux
+    variable = grad_u0_fft
     fft_buffer = grad_u
     component = 0
     execute_on = FINAL
   [../]
-  [./grad_u1_aux]
+  [./grad_u1_fft]
     type = FFTBufferAux
-    variable = grad_u1_aux
+    variable = grad_u1_fft
     fft_buffer = grad_u
     component = 1
     execute_on = FINAL
   [../]
+
+  [./grad_u0_fem]
+    type = GradientComponentAux
+    variable = grad_u0_fem
+    v = u_fem
+    component = 0
+    execute_on = FINAL
+  [../]
+  [./grad_u1_fem]
+    type = GradientComponentAux
+    variable = grad_u1_fem
+    v = u_fem
+    component = 1
+    execute_on = FINAL
+  [../]
 []
 
 [Executioner]

include/auxkernels/GradientComponentAux.h

@@ -0,0 +1,31 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#pragma once
+
+#include "AuxKernel.h"
+
+/**
+ * Get a selected component of the gradient of a coupled variable
+ */
+class GradientComponentAux : public AuxKernel
+{
+public:
+  static InputParameters validParams();
+
+  GradientComponentAux(const InputParameters & parameters);
+
+protected:
+  virtual Real computeValue() override;
+
+  /// Gradient of the coupled variable
+  const VariableGradient & _grad_v;
+
+  /// Component of the gradient vector to match
+  const unsigned int _component;
+};

src/auxkernels/GradientComponentAux.C

@@ -0,0 +1,38 @@
+/**********************************************************************/
+/*                     DO NOT MODIFY THIS HEADER                      */
+/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
+/*                                                                    */
+/*            Copyright 2017 Battelle Energy Alliance, LLC            */
+/*                        ALL RIGHTS RESERVED                         */
+/**********************************************************************/
+
+#include "GradientComponentAux.h"
+
+registerMooseObject("MagpieApp", GradientComponentAux);
+
+InputParameters
+GradientComponentAux::validParams()
+{
+  InputParameters params = AuxKernel::validParams();
+  params.addClassDescription("Return specified component of the gradient of a coupled variable.");
+  params.addRequiredCoupledVar("v", "Coupled variable to extract gradient component of");
+
+  params.addRequiredParam<unsigned int>("component",
+                                        "Component of the gradient of the coupled variable v");
+  return params;
+}
+
+GradientComponentAux::GradientComponentAux(const InputParameters & parameters)
+  : AuxKernel(parameters),
+    _grad_v(coupledGradient("v")),
+    _component(getParam<unsigned int>("component"))
+{
+  if (_component >= LIBMESH_DIM)
+    paramError("component", "Component too large for LIBMESH_DIM");
+}
+
+Real
+GradientComponentAux::computeValue()
+{
+  return _grad_v[_qp](_component);
+}

src/executioners/SpectralExecutionerBase.C

@@ -72,7 +72,6 @@ SpectralExecutionerBase::execute()
   u_buffer.forward();
 
   auto & grad_u_buffer = getFFTBuffer<RealVectorValue>("grad_u");
-  auto grad_u = grad_u_buffer.realSpace();
   kVectorMultiply(u_buffer, grad_u_buffer);
 
   _time_step = 1;
@@ -92,6 +91,7 @@ SpectralExecutionerBase::execute()
   grad_u_buffer.backward();
 
   u /= 10000.0;
+  auto grad_u = grad_u_buffer.realSpace();
   grad_u /= 100.0;
 
   _time_step = 2;
@@ -111,14 +111,15 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase<Real> & in_buffer,
   mooseAssert(in.size() == out.size(), "Buffer sizes must be equal");
 
   const auto & grid = in_buffer.grid();
+  const Complex I(0.0, 1.0);
   switch (in_buffer.dim())
   {
     case 1:
     {
       const auto & ivec = in_buffer.kTable(0);
       const int ni = grid[0];
       for (int i = 0; i * 2 <= ni; ++i)
-        out[i](0) = in[i] * ivec[i];
+        out[i](0) = in[i] * ivec[i] * I;
       return;
     }
 
@@ -128,12 +129,12 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase<Real> & in_buffer,
       const auto & ivec = in_buffer.kTable(0);
       const auto & jvec = in_buffer.kTable(1);
       const int ni = grid[0];
-      const int nj = grid[1];
+      const int nj = (grid[1] >> 1) + 1;
       for (int i = 0; i < ni; ++i)
-        for (int j = 0; j * 2 <= nj; ++j)
+        for (int j = 0; j < nj; ++j)
         {
-          out[index](0) = in[index] * ivec[i];
-          out[index](1) = in[index] * jvec[j];
+          out[index](0) = in[index] * ivec[i] * I;
+          out[index](1) = in[index] * jvec[j] * I;
           index++;
         }
       return;
@@ -152,9 +153,9 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase<Real> & in_buffer,
         for (int j = 0; j < nj; ++j)
           for (int k = 0; k * 2 <= nk; ++k)
           {
-            out[index](0) = in[index] * ivec[i];
-            out[index](1) = in[index] * jvec[j];
-            out[index](2) = in[index] * kvec[k];
+            out[index](0) = in[index] * ivec[i] * I;
+            out[index](1) = in[index] * jvec[j] * I;
+            out[index](2) = in[index] * kvec[k] * I;
             index++;
           }
       return;

src/userobjects/FFTBufferBase.C

@@ -118,7 +118,7 @@ FFTBufferBase<T>::FFTBufferBase(const InputParameters & parameters)
     // precompute kvector components for current direction
     _k_table[i].resize(_grid[i]);
     for (int j = 0; j < _grid[i]; ++j)
-      _k_table[i][j] = (j * 2 > _grid[i] ? Real(_grid[i] - j) : Real(j)) / _box_size(i);
+      _k_table[i][j] = ((j < (_grid[i] >> 1) + 1) ? Real(j) : (Real(j) - _grid[i])) / _box_size(i);
   }
   _real_space_data.resize(real_space_buffer_size);
   _reciprocal_space_data.resize(reciprocal_space_buffer_size);