splines made but need to be checked, solver seems to make wrong results

blas.hpp

@@ -0,0 +1,54 @@
+#ifndef NUFI_BLAS_HPP
+#define NUFI_BLAS_HPP
+
+#include <cstddef>
+
+/*!
+ * \brief Convenience wrappers for BLAS, with overloads for single and double
+ *        precision.
+ */
+namespace blas
+{
+
+double dot( const size_t n, const double *x, size_t incx,
+                            const double *y, size_t incy );
+
+float  dot( const size_t n, const float  *x, size_t incx,
+                            const float  *y, size_t incy );
+
+void axpy( size_t n, double alpha, const double *x, size_t incx,
+                                         double *y, size_t incy );
+
+void axpy( size_t n, float  alpha, const float  *x, size_t incx,
+                                         float  *y, size_t incy );
+
+
+void scal( size_t n, double alpha, double *x, size_t incx );
+void scal( size_t n, float  alpha, float  *x, size_t incx );
+
+void copy( size_t n, const double *x, size_t incx, double *y, size_t incy );
+void copy( size_t n, const float  *x, size_t incx, float  *y, size_t incy );
+
+void ger( const size_t M, const size_t N, const double alpha,
+          const double *X, const size_t incX, const double *Y, const size_t incY,
+          double *A, const size_t lda);
+
+void ger( const size_t M, const size_t N, const float alpha,
+          const float *X, const size_t incX, const float *Y, const size_t incY,
+          float *A, const size_t lda);
+
+
+void gemv( const char trans, size_t m, size_t n,
+		   double alpha, const double *a, size_t lda,
+           const double *x, size_t incx, double beta,
+           double *y, size_t incy );
+
+void gemv( const char trans, size_t m, size_t n,
+		   float alpha, const float *a, size_t lda,
+           const float *x, size_t incx, float beta,
+           float *y, size_t incy );
+}
+
+
+#endif
+

fields.hpp

@@ -4,6 +4,8 @@
 #include <cmath>
 #include <deal.II/base/function.h>
 #include "parameters.hpp"
+#include "splines.hpp"
+#include "lsmr.hpp"
 
 using namespace dealii;
 
@@ -36,8 +38,8 @@ inline double compute_rho(const double x,
 }
 
 
-template <typename real, size_t order, size_t dx = 0>
+template <size_t dx = 0>
-real eval( real x, const real *coeffs) noexcept
+double eval(double x, const double *coeffs) noexcept
 {
     using std::floor;
 
@@ -48,7 +50,7 @@ real eval( real x, const real *coeffs) noexcept
     x = x - Parameters::LX * floor( x/Parameters::LX ); 
 
     // Knot number
-    real x_knot = floor( x/Parameters::SPLINE_DX); 
+    double x_knot = floor( x/Parameters::SPLINE_DX); 
 
     size_t ii = static_cast<size_t>(x_knot);
 
@@ -56,10 +58,10 @@ real eval( real x, const real *coeffs) noexcept
     x = x/Parameters::SPLINE_DX - x_knot;
 
     // Scale according to derivative.
-    real factor = 1;
+    double factor = 1;
     for ( size_t i = 0; i < dx; ++i ) factor *= 1/Parameters::SPLINE_DX;
 
-    return factor*splines1d::eval<real,order,dx>( x, coeffs + ii );
+    return factor*splines1d::eval<double,Parameters::SPLINE_ORDER,dx>(x, coeffs + ii);
 }
 
 template <typename real, size_t order>
@@ -72,17 +74,15 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be
 
     struct mat_t // STRUCT AND CONFIG NEEDS TO BE REVIEWED 
     {
-        const config_t<real> &config;
         real  N[ order ];
 
-        mat_t( const config_t<real> &conf ): config { conf }
+        mat_t()
         {
             splines1d::N<real,order>(0,N);
         }
 
         void operator()( const real *in, real *out ) const
         {
-            #pragma omp parallel for 
             for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
             {
                 real result = 0;
@@ -103,10 +103,9 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be
 
     struct transposed_mat_t // STRUCT AND CONFIG NEEDS TO BE REVIEWED 
     {
-        const config_t<real> &config;
         real  N[ order ];
 
-        transposed_mat_t( const config_t<real> &conf ): config { conf }
+        transposed_mat_t()
         {
             splines1d::N<real,order>(0,N);
         }
@@ -132,9 +131,9 @@ void interpolate( real *coeffs, const real *values) // Least Squares needs to be
         }
     };
 
-    mat_t M { config }; transposed_mat_t Mt { config };
+    mat_t M; transposed_mat_t Mt;
     lsmr_options<real> opt; opt.silent = true;
-    lsmr( config.Nx, config.Nx, M, Mt, values, tmp.get(), opt );
+    lsmr( Parameters::SPLINE_NX, Parameters::SPLINE_NX , M, Mt, values, tmp.get(), opt );
 
     if ( opt.iter == opt.max_iter )
         std::cerr << "Warning. LSMR did not converge.\n";

lsmr.hpp

@@ -0,0 +1,252 @@
+#ifndef LSMR_HPP
+#define LSMR_HPP
+
+#include <cmath>
+#include <limits>
+#include <iomanip>
+#include <iostream>
+#include "blas.hpp"
+
+template <typename real>
+struct lsmr_options
+{
+    ///////////
+    // INPUT //
+    ///////////
+
+    // Whether to print messages to std::cout.
+    bool silent            = true;
+
+    // Residual of normal equations AᵀAx = Aᵀb
+    bool relative_residual  = true;
+    real target_residual    = std::numeric_limits<real>::epsilon();
+    size_t      max_iter    = 1000;
+
+    // How many Lánczos vectors to keep for local reorthogonalisation.
+    // Choose zero for no reorthogonalisation, pure LSMR.
+    // Choose a large value for complete reorthognalisation.
+    //
+    // In an ideal world without roundoff errors, this would have no effect
+    // at all, as the Lánczos vectors would be perfectly orthogonal. In practice
+    // this property is lost rather quickly.  One may choose to store some of
+    // the most recent Lánczos vectors to enforce this property manually. This
+    // increase convergence speed at the cost of additional memory requirements.
+    size_t reorthogonalise_u = 50;
+    size_t reorthogonalise_v = 50;
+
+    ////////////
+    // OUTPUT //
+    ////////////
+
+    // Iteration count and reached residual.
+    // Estimates of ‖A‖ and cond(A)
+    size_t iter; real residual;
+    real norm_A_estimate, cond_estimate;
+};
+
+template <typename real, typename mat, typename transposed_mat>
+void lsmr( size_t m, size_t n, const mat& A, const transposed_mat& At,
+           const real *b, real *x, lsmr_options<real> &S );
+
+namespace lsmr_impl
+{
+
+template <typename real>
+real norm( size_t n, const real *x )
+{
+    using std::hypot;
+
+    real result = 0;
+    for ( size_t i = 0; i < n; ++i )
+        result = hypot(result,x[i]);
+
+    return result;
+}
+
+// Reorthognalise u with respect to the previous vectors in buffer,
+// using the modified Gram–Schmidt process. Overwrite the oldest vector
+// in buffer when full.
+template <typename real>
+void reorthogonalise( real *buf, size_t n, size_t buffer_max,  
+                      real *u,   size_t iter )
+{
+    using std::min;
+    using blas::dot;
+    using blas::axpy;
+    using blas::scal;
+    using blas::copy;
+
+    size_t n_buffered = min( iter+1, buffer_max );
+    for ( size_t i = 0; i < n_buffered; ++i )
+    {
+        real fac = -dot( n, u, 1, buf + i*n, 1 );
+        axpy( n, fac, buf + i*n, 1, u, 1 );
+    }
+
+    scal( n, 1/norm(n,u), u, 1 );
+    copy( n, u, 1, buf + ((iter+1)%buffer_max)*n, 1 );
+}
+
+}
+
+template <typename real, typename mat, typename transposed_mat>
+void lsmr( size_t m, size_t n, const mat& A, const transposed_mat& At,
+           const real *b, real *x, lsmr_options<real> &S )
+{
+    using std::min;
+    using std::max;
+    using std::abs;
+    using std::swap;
+    using std::hypot;
+    using blas::axpy; 
+    using blas::scal; 
+    using blas::copy;
+    using lsmr_impl::norm;
+    using lsmr_impl::reorthogonalise;
+
+
+    // Allocation of buffers.
+    size_t max_buf = min(n,m)-1;
+    S.reorthogonalise_u = min(S.reorthogonalise_u,max_buf);
+    S.reorthogonalise_v = min(S.reorthogonalise_v,max_buf);
+    size_t u_buffer_size = max( S.reorthogonalise_u, size_t(1) );
+    size_t v_buffer_size = max( S.reorthogonalise_v, size_t(1) );
+
+    std::unique_ptr<real[]> data { new real[ n*( 4 + v_buffer_size )  +
+                                             m*( 2 + u_buffer_size ) ] {} };
+
+    real *u     = data.get();
+    real *utmp  = u     + m;
+    real *ubuf  = utmp  + m;
+    real *v     = ubuf  + m*u_buffer_size;
+    real *vtmp  = v     + n;
+    real *h     = vtmp  + n;
+    real *h_bar = h     + n;
+    real *vbuf  = h_bar + n;
+
+    At(b,v);
+    const real norm_ATb = norm(n,v);
+
+
+    A(x,u); axpy(m,real(-1),b,1,u,1); 
+    scal(m, real(-1), u, 1 );         // u = b - Ax;
+   
+    real alpha = 0; 
+    real beta  = norm(m,u); 
+
+    if ( beta > real(0) )
+    {
+        scal(m, real(1)/beta, u, 1 ); // u = b - Ax / norm(b-Ax)
+        At(u,v);                      // v = At*u
+        alpha = norm(n,v);
+    }
+
+    if ( alpha > real(0) )
+        scal(n, real(1)/alpha, v, 1 ); // v = At*u/norm(At*u)
+
+    copy(n,u,1,ubuf,1);                // u_buf.col(0) = u_buf
+    copy(n,v,1,vbuf,1);                // v_buf.col(0) = v
+    copy(n,v,1,h,1);                   // h = v
+
+    if ( alpha * beta == real(0) ) return;
+
+    
+    real alpha_bar = alpha, zeta_bar = alpha*beta;
+    real rho = 1, rho_bar = 1, c_bar = 1, s_bar = 0;
+    real c, s, theta, zeta, theta_bar, rho_prev, rho_bar_prev;
+
+    // For estimating the condition number.
+    real   sigma_max = 0,   sigma_min = std::numeric_limits<real>::max();
+    real rho_bar_max = 0, rho_bar_min = std::numeric_limits<real>::max();
+
+    S.norm_A_estimate = 0;
+    for ( S.iter = 0; S.iter < S.max_iter; ++S.iter )
+    {
+        // Continue the bidiagonalisation.
+        A(v,utmp); axpy(m,-alpha,u,1,utmp,1); swap(u,utmp); // u = A*v - alpha*u
+        beta = norm(m,u);
+
+        if ( beta > 0 )
+        {
+            scal(m, real(1)/beta, u, 1 );
+            if ( S.reorthogonalise_u )
+                reorthogonalise( ubuf, m, u_buffer_size, u, S.iter );
+
+            S.norm_A_estimate = hypot( alpha, S.norm_A_estimate );
+            S.norm_A_estimate = hypot( beta , S.norm_A_estimate );
+
+            At(u,vtmp); axpy(n,-beta,v,1,vtmp,1); swap(v,vtmp); // v = At*u - beta*v
+            alpha = norm(n,v);
+
+            if ( alpha > 0 )
+            {
+                scal(n,real(1)/alpha, v, 1 );
+                if ( S.reorthogonalise_v )
+                    reorthogonalise( vbuf, n, v_buffer_size, v, S.iter );
+            }
+        }
+
+        // Construct and apply rotation P_k
+        rho_prev     = rho;
+        rho          = hypot(alpha_bar,beta);
+        c            = alpha_bar/rho;
+        s            = beta/rho;
+        theta        = s*alpha;
+        alpha_bar    = c*alpha;
+
+        // Construct and apply rotation \bar{P}_k
+        rho_bar_prev = rho_bar;
+        if ( S.iter )
+        {
+            rho_bar_max = max( rho_bar, rho_bar_max );
+            rho_bar_min = min( rho_bar, rho_bar_min );
+        }
+        theta_bar    = s_bar*rho;
+        rho_bar      = hypot( c_bar*rho, theta );
+        if ( S.iter )
+        {
+            sigma_max    = max( rho_bar_max, c_bar*rho );
+            sigma_min    = min( rho_bar_min, c_bar*rho );
+        }
+        c_bar        = c_bar * rho/rho_bar;
+        s_bar        = theta/rho_bar;
+        zeta         = c_bar * zeta_bar;
+        zeta_bar     = -s_bar*zeta_bar;
+
+
+        // Update h, h_bar, x
+        scal(n, -(theta_bar*rho)/(rho_prev*rho_bar_prev), h_bar, 1 ) ;
+        axpy(n, real(1), h, 1, h_bar, 1 );             // h_bar = h - factor*h_bar
+
+        axpy( n, zeta/(rho*rho_bar), h_bar, 1, x, 1 ); // x += factor * h_bar
+
+        scal(n, -theta/rho, h, 1 );
+        axpy(n, real(1), v, 1, h, 1 );                 // h = v - factor*h;
+
+        // Estimate quantities.
+        if ( S.relative_residual ) S.residual = abs(zeta_bar)/norm_ATb;
+        else                       S.residual = abs(zeta_bar);
+        S.cond_estimate = sigma_max / sigma_min;
+
+        if ( S.residual <= S.target_residual )
+        {
+            if ( S.silent == false )
+            {
+                std::cout << "LSMR: Iteration:  " << std::setw(4)  << S.iter << ", "
+                          << "Residual:  "        << std::setw(12) << std::scientific << S.residual << ", "
+                          << "cond estimate:  "   << std::setw(12) << std::scientific << S.cond_estimate << ".\n";
+            }
+            return;
+        }
+
+        if ( S.silent == false && (S.iter%10) == 0 )
+        {
+            std::cout << "LSMR: Iteration:  " << std::setw(4)  << S.iter << ", "
+                      << "Residual:  "        << std::setw(12) << std::scientific << S.residual << ", "
+                      << "cond estimate:  "   << std::setw(12) << std::scientific << S.cond_estimate << ".\n";
+        }
+    }
+}
+
+#endif
+

nufi_solver.hpp

@@ -8,15 +8,14 @@
 #include <deal.II/numerics/fe_field_function.h>
 
 #include <iostream>
+#include <memory>
 #include <ostream>
 #include <vector>
 #include <cstddef>
 
 #include "parameters.hpp"
 #include "poisson_problem.hpp"
-#include "fields.hpp" // holds f0(x,v), and compute_rho(x)
+#include "fields.hpp"
-#include "spline_field.hpp" // old GPT splines
-#include "splines.hpp"    //new splines                       
 
 using namespace dealii;
 
@@ -26,15 +25,15 @@ public:
   NuFISolver();
 
   void run();
-  double eval_rho(unsigned int n, double x, const std::vector<double> E_coeffs, unsigned int Nv = Parameters::NV);
+  double eval_rho(unsigned int n, double x, const double *E_coeffs, unsigned int Nv = Parameters::NV);
-  double eval_ftilda(unsigned int n, double x, double u, const std::vector<double> E_coeffs);
+  double eval_ftilda(unsigned int n, double x, double u, const double *E_coeffs);
-  void save_ftilda(unsigned int n, const std::vector<double> E_coeffs, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename);
+  void save_ftilda(unsigned int n, const double *E_coeffs, unsigned int Nx_out, unsigned int Nv_out, const std::string &filename);
 
 private:
 
 
   unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
-  [[maybe_unused]] unsigned int Nx = Parameters::SPLINE_NX;
+  unsigned int Nx = Parameters::SPLINE_NX;
 
   double Lx = Parameters::LX;
 
@@ -42,8 +41,6 @@ private:
 
   unsigned int order;
 
-  double dt = Parameters::DT;
-
   PoissonProblem<1> poisson;
 
 };
@@ -51,72 +48,82 @@ private:
 inline double NuFISolver::eval_ftilda(unsigned int n,
                                       double x,
                                       double u,
-                                      const std::vector<double> E_coeffs)
+                                      const double *E_coeffs)
+
 {
-  double Lu = std::abs(Parameters::V_DOMAIN_LEFT - Parameters::V_DOMAIN_RIGHT);
+  if ( n == 0 ) return f0(x,u);
 
-  if (n == 0)
+  const size_t stride_x = 1;
-      return f0(x, u);
+  const size_t stride_t = stride_x*(Parameters::SPLINE_NX + Parameters::SPLINE_ORDER - 1);
 
-  // Initial half-step.
+  double Ex;
-  u += 0.5*dt*E_spline.eval(x);
+  const double *c;
+
+  // We omit the initial half-step.
 
   while ( --n )
   {
-      x -= dt*u;
+      x  = x - Parameters::DT *u;
-      u += dt*E_spline.eval(x);
+      c  = E_coeffs + n*stride_t;
+      Ex = -eval<1>(x, c);
+      u  = u + Parameters::DT *Ex;
   }
 
-  // Final half-step.
+  // The final half-step.
-  x -= dt*u;
+  x -= Parameters::DT*u;
-  u += 0.5*dt*E_spline.eval(x);
+  c  = E_coeffs + n*stride_t;
+  Ex = -eval<1>(x, c);
+  u += 0.5*Parameters::DT*Ex;
 
-  double x_periodic = x - Lx * std::floor(x / Lx);
+  return f0(x,u);
-  double u_periodic = u - Lu * std::floor(u / Lu);
-
-  return f0(x_periodic, u_periodic);
 }
 
 inline double NuFISolver::eval_rho(const unsigned int n,
                           const double x,
-                          const std::vector<double> E_coeffs,
+                          const double *E_coeffs,
                           const unsigned int Nv)
 {
   const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
+  const double v_min = Parameters::V_DOMAIN_LEFT;
 
   double integral = 0.0;
   for (unsigned int i = 0; i < Nv; ++i)
-  {
+    integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs) * dv;
-    const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
-    AssertThrow(std::isfinite(E_spline.eval(x)), ExcMessage("NaN detected in E_spline.eval(x) inside NuFISolver::eval_rho integral loop"));
-    integral += eval_ftilda(n, x, v, E_spline) * dv;
-  }
   
   return 1.0 - integral;
 }
 
-class ChargeDensity_NuFI : public Function<1>
+template<unsigned int dim>
+class ChargeDensity_NuFI : public Function<dim>
 {
-public:
+  public:
-  ChargeDensity_NuFI(NuFISolver &solver, size_t n, const std::vector<double> E_coeffs)
+    ChargeDensity_NuFI(const double *rho_values, unsigned int Nx)
-      : solver(solver), n(n), E_coeffs(E_coeffs) {}
+      : Function<dim>(), rho(rho_values), Nx(Nx) {}
 
-  virtual double value(const Point<1> &p,
+    virtual double value(const Point<dim> &p,
                          [[maybe_unused]] const unsigned int component = 0) const override
     {
-      double x = p[0];
+      const double x = p[0];
 
-      return solver.eval_rho(n, x, E_coeffs); 
+      // Map x -> grid index
+      const double L = Parameters::LX;
+      const double dx = L / (Nx-1);
+
+      int i = static_cast<int>(std::floor((x - Parameters::X_DOMAIN_LEFT) / dx));
+
+      // periodic wrap
+      i = (i % Nx + Nx) % Nx;
+
+      return rho[i];
     }
 
-private:
+  private:
-  NuFISolver &solver;
+    const double *rho;
-  size_t n;
+    const unsigned int Nx;
-  const std::vector<double> E_coeffs;
 };
 
 inline void NuFISolver::save_ftilda(unsigned int n,
-                                    const std::vector<double> E_coeffs,
+                                    const double *E_coeffs,
                                     unsigned int Nx_out,
                                     unsigned int Nv_out,
                                     const std::string &filename)
@@ -144,7 +151,7 @@ inline void NuFISolver::save_ftilda(unsigned int n,
       {
           double v = vmin + (j + 0.5)*dv;
 
-          double val = eval_ftilda(n, x, v, E_spline);
+          double val = eval_ftilda(n, x, v, E_coeffs);
 
           file << val;
 
@@ -160,60 +167,39 @@ inline void NuFISolver::save_ftilda(unsigned int n,
 
 inline void NuFISolver::run()
 {
-  std::cout << "Start of NuFISolver::run()\n\n";
+  std::cout << "Building E_sline\n\n";
 
-  // init E_spline
+  using std::abs;
+  using std::max;
 
-  unsigned int Nx = Parameters::SPLINE_NX;
+  const size_t stride_t = Nx + order - 1;
-  // Nx grid points
-  double dx = Lx / (Nx-1);
 
-  std::vector<double> E_grid(Nx);
+  std::unique_ptr<double[]> coeffs { new double[ Nt*stride_t ] {} };
+  std::unique_ptr<double,decltype(std::free)*> rho { reinterpret_cast<double*>(std::aligned_alloc(64,sizeof(double)*Nx)), std::free };
 
-  //set initial E points
+  if ( rho == nullptr ) throw std::bad_alloc {};
-  for (unsigned int i=0; i<Nx; ++i)
-  {
-    [[maybe_unused]] double x = Parameters::X_DOMAIN_LEFT + i*dx;
-    E_grid[i] = 0; 
-  }
-
-  std::vector<double> E_coeffs(E_grid, Parameters::X_DOMAIN_LEFT, Parameters::X_DOMAIN_RIGHT); // Needs correction
 
   for (unsigned int it = 0; it < Nt; ++it)
   {
       std::cout << "Timestep " << it << " / " << Nt << std::endl << std::endl;
 
-      // Step 1: Evaluate rho^n(x) using current E_spline
+      // compute rho
-
+    	for(size_t i = 0; i<Nx; i++)
-      rho.resize(Nx);
-
-      for (unsigned int i = 0; i < (Nx); ++i)
     	{
-          double x = (i + 0.5) * dx;
+        double ith_rho = eval_rho(it, i, coeffs.get(), Parameters::NV); 
-          rho[i] = eval_rho(it, x, E_spline, Parameters::NV);
+        AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
-      }
+    		rho.get()[i] =  ith_rho;
-
-      ChargeDensity_NuFI rho_function(*this, it, E_spline);
-
-      poisson.set_rhs_function(rho_function);
-
-      for (unsigned int i=0; i< rho.size(); ++i) // check for bad rho[i]
-      {
-        AssertThrow(std::isfinite(rho[i]), ExcMessage("NaN detected in rho"));
     	}
 
+      poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Parameters::SPLINE_NX));
       poisson.solve_step();
 
       if (it % Parameters::PLOT_FREQUENCY == 0)
       {
         std::cout << "Saving results... \n\n";
-        save_ftilda(it, E_spline, 128, 128, "results/ftilda_" + std::to_string(it) + ".dat");
+        save_ftilda(it, coeffs.get(), 128, 128, "results/ftilda_" + std::to_string(it) + ".dat");
         poisson.output_results(it);
       }
-
-      E_grid = poisson.sample_electric_field(poisson, Nx, 0.0, Lx);
-
-      E_spline = std::vector<double> E_coeffs; // needs correction 
   }
 
   std::cout << "NuFI simulation finished.\n";

parameters.hpp

@@ -3,6 +3,7 @@
 
 #include <cmath>
 #include <cstdlib>
+
 namespace Parameters
 {
   constexpr unsigned int DIMENSION = 1;
@@ -27,9 +28,11 @@ namespace Parameters
   constexpr double DT=1./16.;
   constexpr unsigned int TMAX = 10;
 
+
   //spline options
   constexpr int SPLINE_NX = 512;
   constexpr double SPLINE_DX = LX/SPLINE_NX;
+  constexpr size_t SPLINE_ORDER = 4;
 
   //Plotting options
   constexpr int PLOT_FREQUENCY = 2;

poisson_problem.hpp

@@ -33,7 +33,9 @@
 #include <deal.II/numerics/vector_tools.h>
 #include <deal.II/numerics/fe_field_function.h>
 
+#include <memory>
 #include <string>
+#include <utility>
 #include <vector>
 
 #include "parameters.hpp"
@@ -53,7 +55,7 @@ public:
   void solve_step();
   void run();
 
-  void set_rhs_function(const Function<dim> &rhs);
+  void set_rhs_function(std::unique_ptr<Function<dim>> rhs_function);
 
   const Vector<double> &get_solution() const { return solution; }
   const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
@@ -83,7 +85,7 @@ private:
   Vector<double> solution;      // phi
   Vector<double> system_rhs;
 
-  const Function<dim> *rhs_function;
+  std::unique_ptr<const Function<dim>> rhs_function;
 
   MappingQ<dim> mapping;
 };
@@ -91,9 +93,9 @@ private:
 // Utilities
 
 template <int dim>
-void PoissonProblem<dim>::set_rhs_function(const Function<dim> &rhs)
+void PoissonProblem<dim>::set_rhs_function(std::unique_ptr<Function<dim>> rhs)
 {
-  rhs_function = &rhs;
+  rhs_function = std::move(rhs);
 }
 
 template <int dim>
@@ -214,9 +216,6 @@ public:
   }
 };
 
-
-// =-=-=-=-= Poisson equation solver =-=-=-=-=
-
 template <int dim>
 void PoissonProblem<dim>::assemble_system()
 {