By default scale automatically (idaholab#401)

dschwen · Apr 29, 2020 · 58ccfcc · 58ccfcc
1 parent 31ec260
commit 58ccfcc
Show file tree

Hide file tree

Showing 6 changed files with 110 additions and 24 deletions.
diff --git a/examples/spectral.i b/examples/spectral.i
@@ -2,9 +2,10 @@
   type = MyTRIMMesh
   dim = 2
   xmax = 100
-  ymax = 100
+  ymin = -25
+  ymax = 25
   nx = 100
-  ny = 100
+  ny = 50
 []
 
 [Problem]
@@ -51,13 +52,23 @@
     order = CONSTANT
     family = MONOMIAL
     [./InitialCondition]
-      type = SmoothCircleIC
+      type = SmoothSuperellipsoidIC
       x1 = 50
-      y1 = 50
-      radius = 30
-      int_width = 20
+      y1 = 0
+      a = 30
+      b = 15
+      int_width = 15
       invalue = 1
       outvalue = 0
+      n = 2
+    [../]
+  [../]
+  [./v_fft]
+    order = CONSTANT
+    family = MONOMIAL
+    [./InitialCondition]
+      type = FunctionIC
+      function = if(y<0,y,-y)
     [../]
   [../]
   [./grad_u0_fft]
@@ -68,16 +79,30 @@
     order = CONSTANT
     family = MONOMIAL
   [../]
+  [./grad_v1_fft]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
   
   [./u_fem]
     [./InitialCondition]
-      type = SmoothCircleIC
+      type = SmoothSuperellipsoidIC
       x1 = 50
-      y1 = 50
-      radius = 30
-      int_width = 20
+      y1 = 0
+      a = 30
+      b = 15
+      int_width = 15
       invalue = 1
       outvalue = 0
+      n = 2
+    [../]
+  [../]
+  [./v_fem]
+    order = CONSTANT
+    family = MONOMIAL
+    [./InitialCondition]
+      type = FunctionIC
+      function = if(y<0,y,-y)
     [../]
   [../]
   [./grad_u0_fem]
@@ -88,6 +113,10 @@
     order = CONSTANT
     family = MONOMIAL
   [../]
+  [./grad_v1_fem]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
 []
 
 [Materials]
@@ -116,6 +145,13 @@
   [./grad_u]
     type = RealVectorValueFFTWBuffer
   [../]
+  [./v]
+    type = RealFFTWBuffer
+    moose_variable = v_fft
+  [../]
+  [./grad_v]
+    type = RealVectorValueFFTWBuffer
+  [../]
 []
 
 [AuxKernels]
@@ -154,6 +190,12 @@
     fft_buffer = u
     execute_on = FINAL
   [../]
+  [./v_fft]
+    type = FFTBufferAux
+    variable = v_fft
+    fft_buffer = v
+    execute_on = FINAL
+  [../]
 
   [./grad_u0_fft]
     type = FFTBufferAux
@@ -169,6 +211,13 @@
     component = 1
     execute_on = FINAL
   [../]
+  [./grad_v1_fft]
+    type = FFTBufferAux
+    variable = grad_v1_fft
+    fft_buffer = grad_v
+    component = 1
+    execute_on = FINAL
+  [../]
 
   [./grad_u0_fem]
     type = GradientComponentAux

diff --git a/include/userobjects/FFTBufferBase.h b/include/userobjects/FFTBufferBase.h
@@ -45,9 +45,19 @@ class FFTBufferBase : public ElementUserObject
   virtual void finalize() {}
   virtual void threadJoin(const UserObject &) {}
 
-  ///@{ transforms
-  virtual void forward() = 0;
-  virtual void backward() = 0;
+  ///@{ transforms with proper scaling guaranteed
+  virtual void forward();
+  virtual void backward();
+  ///@}
+
+  ///@{ transforms without proper scaling guaranteed
+  virtual void forwardRaw() = 0;
+  virtual void backwardRaw() = 0;
+  ///@}
+
+  ///@{ scaling required after respective transform
+  virtual Real forwardScale() { return 1.0; }
+  virtual Real backwardScale() { return forwardScale(); }
   ///@}
 
   ///@{ buffer access

diff --git a/include/userobjects/FFTWBufferBase.h b/include/userobjects/FFTWBufferBase.h
@@ -28,8 +28,11 @@ class FFTWBufferBase : public FFTBufferBase<T>
   ~FFTWBufferBase();
 
   // transforms
-  void forward() override;
-  void backward() override;
+  void forwardRaw() override;
+  void backwardRaw() override;
+
+  // scaling
+  Real forwardScale() override { return _scaling; }
 
 protected:
   ///@{ FFTW plans
@@ -42,6 +45,9 @@ class FFTWBufferBase : public FFTBufferBase<T>
   PerfID _perf_fft;
   ///@}
 
+  /// scale factor
+  Real _scaling;
+
   usingFFTBufferBaseMembers;
 };
 

diff --git a/src/executioners/SpectralExecutionerBase.C b/src/executioners/SpectralExecutionerBase.C
@@ -59,20 +59,25 @@ SpectralExecutionerBase::execute()
   // back and forth test
   auto & c_buffer = getFFTBuffer<Real>("c");
   auto c = c_buffer.realSpace();
-  auto c_tilde = c_buffer.reciprocalSpace();
+  // auto c_tilde = c_buffer.reciprocalSpace();
   c_buffer.forward();
 
   auto & R_buffer = getFFTBuffer<RealVectorValue>("R");
   auto R = R_buffer.realSpace();
   R_buffer.forward();
 
-  // gradient test
+  // gradient tests
   auto & u_buffer = getFFTBuffer<Real>("u");
   auto u = u_buffer.realSpace();
   u_buffer.forward();
+  auto & v_buffer = getFFTBuffer<Real>("v");
+  auto v = v_buffer.realSpace();
+  v_buffer.forward();
 
   auto & grad_u_buffer = getFFTBuffer<RealVectorValue>("grad_u");
   kVectorMultiply(u_buffer, grad_u_buffer);
+  auto & grad_v_buffer = getFFTBuffer<RealVectorValue>("grad_v");
+  kVectorMultiply(v_buffer, grad_v_buffer);
 
   _time_step = 1;
   _fe_problem.execute(EXEC_FINAL);
@@ -83,16 +88,14 @@ SpectralExecutionerBase::execute()
   // back and forth test
   c_buffer.backward();
   R_buffer.backward();
-  R /= 10000.0;
-  c /= 10000.0;
 
   // gradient test
   u_buffer.backward();
   grad_u_buffer.backward();
+  v_buffer.backward();
+  grad_v_buffer.backward();
 
-  u /= 10000.0;
-  auto grad_u = grad_u_buffer.realSpace();
-  grad_u /= 100.0;
+  // auto grad_u = grad_u_buffer.realSpace();
 
   _time_step = 2;
   _fe_problem.execute(EXEC_FINAL);

diff --git a/src/userobjects/FFTBufferBase.C b/src/userobjects/FFTBufferBase.C
@@ -156,6 +156,22 @@ FFTBufferBase<T>::execute()
     _real_space_data_start[a * _how_many + i] = (*_moose_variable[i])[0];
 }
 
+template <typename T>
+void
+FFTBufferBase<T>::forward()
+{
+  forwardRaw();
+  _reciprocal_space_data *= forwardScale();
+}
+
+template <typename T>
+void
+FFTBufferBase<T>::backward()
+{
+  backwardRaw();
+  _real_space_data *= backwardScale();
+}
+
 template <typename T>
 const T &
 FFTBufferBase<T>::operator()(const Point & p) const

diff --git a/src/userobjects/FFTWBufferBase.C b/src/userobjects/FFTWBufferBase.C
@@ -52,6 +52,8 @@ FFTWBufferBase<T>::FFTWBufferBase(const InputParameters & parameters)
                                1,
                                FFTW_ESTIMATE);
   }
+
+  _scaling = 1.0 / std::sqrt(_real_space_data.size());
 }
 
 template <typename T>
@@ -64,7 +66,7 @@ FFTWBufferBase<T>::~FFTWBufferBase()
 
 template <typename T>
 void
-FFTWBufferBase<T>::forward()
+FFTWBufferBase<T>::forwardRaw()
 {
   // execute plan
   {
@@ -75,7 +77,7 @@ FFTWBufferBase<T>::forward()
 
 template <typename T>
 void
-FFTWBufferBase<T>::backward()
+FFTWBufferBase<T>::backwardRaw()
 {
   // execute plan
   {