vcbferreira/NuFI_deal.ii
1
0
Fork 0
mirror of https://codeberg.org/vcbferreira/NuFI_deal.ii synced 2026-08-12 14:33:18 +02:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3 from vascocbf/eval_with_fem

Browse Source 
Eval with fem
This commit is contained in:
Vasco C. B. Ferreira
2026-06-30 10:24:58 +02:00
committed by GitHub
parent 5a8622c7fe 96a91c06e9
commit e09c00cc02
18 changed files with 482 additions and 1200 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
+1 -1
View File
@@ -23,7 +23,7 @@ find_package(OpenMP REQUIRED)
add_library(nufi_lib
src/nufi_solver.cc
src/save_results.cc
src/blas.cc
# src/blas.cc
)
target_include_directories(nufi_lib PUBLIC
Makefile
-27
View File
@@ -142,30 +142,6 @@ nufi_poisson/fast:
$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_poisson.dir/build.make CMakeFiles/nufi_poisson.dir/build
.PHONY : nufi_poisson/fast
src/blas.o: src/blas.cc.o
.PHONY : src/blas.o
# target to build an object file
src/blas.cc.o:
$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.o
.PHONY : src/blas.cc.o
src/blas.i: src/blas.cc.i
.PHONY : src/blas.i
# target to preprocess a source file
src/blas.cc.i:
$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.i
.PHONY : src/blas.cc.i
src/blas.s: src/blas.cc.s
.PHONY : src/blas.s
# target to generate assembly for a file
src/blas.cc.s:
$(MAKE) $(MAKESILENT) -f CMakeFiles/nufi_lib.dir/build.make CMakeFiles/nufi_lib.dir/src/blas.cc.s
.PHONY : src/blas.cc.s
src/main.o: src/main.cc.o
.PHONY : src/main.o
@@ -248,9 +224,6 @@ help:
@echo "... rebuild_cache"
@echo "... nufi_lib"
@echo "... nufi_poisson"
@echo "... src/blas.o"
@echo "... src/blas.i"
@echo "... src/blas.s"
@echo "... src/main.o"
@echo "... src/main.i"
@echo "... src/main.s"
README.md
+2 -1
View File
@@ -6,4 +6,5 @@ This simulation of the Vlasov-Poisson system in 1x1v dimensions uses
---
Todo:
- Correct spline interpolator/evaluation for discontinuity
- fix eval and saving fields.
- something with periodicity or eval range in x
libnufi_lib.a
BIN
View File
Binary file not shown.
nufi/blas.h
-54
View File
@@ -1,54 +0,0 @@
#ifndef NUFI_BLAS_H
#define NUFI_BLAS_H
#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
nufi/fields.h
+82 -199
View File
@@ -1,229 +1,112 @@
#ifndef FIELDS_H
#define FIELDS_H
#include "nufi/parameters.h"
#include "nufi/poisson_problem.h"
#include <cmath>
#include <deal.II/base/function.h>
#include "nufi/parameters.h"
#include "nufi/splines.h"
#include "nufi/lsmr.h"
#include <deal.II/base/point.h>
using namespace dealii;
inline double f0(const double x,
const double v,
// inline std::vector<int> Indices_of_points(const std::vector<double> &points,
// double x_min, double x_max, double dx, int grid_type=0)
// {
// // grid type:
// // 0 => uniform
// // 1 => non uniform (TODO)
//
// if (dx <= 0.0) {
// throw std::invalid_argument("dx must be positive");
// }
// if (x_max <= x_min) {
// throw std::invalid_argument("x_max must be > x_min");
// }
//
// std::vector<int> indices;
// indices.reserve(points.size());
//
// switch (grid_type) {
// case 0:
// {
// const double L = x_max - x_min;
// const int N = std::floor(L/dx);
//
//
// for (double x : points) //GPT loop, to check
// {
// x-= x_min;
// x = x - L * std::floor(x/L);
//
// int i = static_cast<int>(std::floor(x / dx));
//
// // safety: handle rare edge case due to floating precision
// if (i == N) i = 0;
//
// indices.push_back(i);
// }
// }
// case 1:
// {
// throw std::invalid_argument("Case for non uniform grid is not
// completed");
// }
// default:
// throw std::invalid_argument("Invalid grid_type argument");
//
// }
// return indices;
// }
inline double f0(const double x, const double v,
const double eps = Parameters::EPS,
const double k = Parameters::WAVE_NR)
{
const double prefactor = Parameters::F0_FACTOR * (1.0 + eps * std::cos(k*x));
const double gaussian = v*v * std::exp(-0.5 * v*v);
const double k = Parameters::WAVE_NR) {
const double prefactor =
Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
const double gaussian = v * v * std::exp(-0.5 * v * v);
return prefactor * gaussian;
}
// wrapper for eval_point() { VectorTools::point_values() }
inline double eval(double x, const PoissonProblem<1> &poisson,
const Vector<double> &solution) noexcept {
inline double compute_rho(const double x,
const unsigned int Nv = Parameters::NV)
{
const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
x -= Parameters::X_DOMAIN_LEFT;
double integral = 0.0;
x = x - Parameters::LX * std::floor(x * Parameters::LX_INV); // in domain
for (unsigned int i = 0; i < Nv; ++i)
{
const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
integral += f0(x, v) * dv;
}
return 1.0 - integral;
return eval_point_grad<1>(poisson.get_mapping(), poisson.get_dof_handler(),
solution, Point<1>(x));
}
template <size_t dx = 0>
double eval(double x, const double *coeffs) noexcept
{
using std::floor;
// Shift to a box that starts at 0.
x -= Parameters::X_DOMAIN_LEFT;
// Get "periodic position" in box at origin.
x = x - Parameters::LX * floor( x*Parameters::LX_INV );
// Knot number
double x_knot = floor( x*Parameters::SPLINE_DX_INV);
size_t ii = static_cast<size_t>(x_knot);
// Convert x to reference coordinates.
x = x*Parameters::SPLINE_DX_INV - x_knot;
// Scale according to derivative.
double factor = 1;
for ( size_t i = 0; i < dx; ++i ) factor *= 1*Parameters::SPLINE_DX_INV;
return factor*splines1d::eval<double,Parameters::SPLINE_ORDER,dx>(x, coeffs + ii);
}
template <typename real, size_t order>
void interpolate( real *coeffs, const real *values)
{
std::unique_ptr<real[]> tmp { new real[ Parameters::SPLINE_NX ] };
for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
tmp[ i ] = coeffs[ i ];
struct mat_t
{
real N[ order ];
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;
if ( i + order <= Parameters::SPLINE_NX )
{
for ( size_t ii = 0; ii < order; ++ii )
result += N[ii] * in[ i + ii ];
}
else
{
for ( size_t ii = 0; ii < order; ++ii )
result += N[ii]*in[ (i+ii) % Parameters::SPLINE_NX];
}
out[ i ] = result;
}
}
};
struct transposed_mat_t
{
real N[ order ];
transposed_mat_t()
{
splines1d::N<real,order>(0,N);
}
void operator()( const real *in, real *out ) const
{
for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
out[ i ] = 0;
for ( size_t i = 0; i < Parameters::SPLINE_NX; ++i )
{
if ( i + order <= Parameters::SPLINE_NX )
{
for ( size_t ii = 0; ii < order; ++ii )
out[ i + ii ] += N[ii] * in[ i ];
}
else
{
for ( size_t ii = 0; ii < order; ++ii )
out[ (i+ii) % Parameters::SPLINE_NX ] += N[ii]*in[ i ];
}
}
}
};
mat_t M; transposed_mat_t Mt;
lsmr_options<real> opt; opt.silent = true;
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";
for ( size_t i = 0; i < Parameters::SPLINE_NX + order - 1; ++i )
coeffs[ i ] = tmp[ i % Parameters::SPLINE_NX ];
}
inline double integral_space_vector(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
{
inline double integral_space_vector(const PoissonProblem<1> &poisson,
const Vector<double> &solution,
double dx = Parameters::PLOT_DX,
size_t Nx = Parameters::PLOT_NX) {
double integral = 0.0;
double xmin = Parameters::X_DOMAIN_LEFT;
#pragma omp parallel for reduction (+:integral)
for (size_t i=0; i<Nx ; ++i) {
#pragma omp parallel for reduction(+ : integral)
for (size_t i = 0; i < Nx; ++i) {
double x = xmin + i * dx;
integral += eval<1>(x, current_coeffs);
integral += eval(x, poisson, solution);
}
return integral*dx;
return integral * dx;
};
inline double integral_space_vector_squared(const double *current_coeffs, double dx = Parameters::SPLINE_DX, size_t Nx = Parameters::SPLINE_NX)
{
inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
const Vector<double> &solution,
double dx = Parameters::PLOT_DX,
size_t Nx = Parameters::PLOT_NX) {
double integral = 0.0;
double xmin = Parameters::X_DOMAIN_LEFT;
#pragma omp parallel for reduction (+:integral)
for (size_t i=0; i<Nx ; ++i) {
double x = xmin + i*dx;
double val = eval<1>(x, current_coeffs);
integral += val*val;
#pragma omp parallel for reduction(+ : integral)
for (size_t i = 0; i < Nx; ++i) {
double x = xmin + i * dx;
double val = eval(x, poisson, solution);
integral += val * val;
}
return integral*dx;
return integral * dx;
};
class Gradient {
public:
Gradient(double xmin, double xmax, unsigned int Nx)
: xmin_(xmin), xmax_(xmax), Nx_(Nx)
{
if (xmax_ <= xmin_) {
throw std::invalid_argument("xmax must be greater than xmin");
}
}
std::vector<double> compute(const std::vector<double>& values) const {
size_t n = values.size();
if (n < 2) {
throw std::invalid_argument("Need at least 2 points");
}
std::vector<double> grad(n);
double dx = (xmax_ - xmin_) / (n-1);
// periodic boundaries
grad[0] = -(values[1] - values[n-1]) / (2.0 * dx);
grad[n-1] = -(values[0] - values[n-2]) / (2.0 * dx);
for (size_t i = 1; i < n-1; ++i) {
grad[i] = -(values[i+1] - values[i-1]) / (2.0 * dx);
}
return grad;
}
private:
double xmin_;
double xmax_;
[[maybe_unused]] unsigned int Nx_;
};
template <int dim>
class ChargeDensity : public Function<dim> // only uses f0
{
public:
ChargeDensity(double eps,
double k,
unsigned int Nv)
: Function<dim>(1), eps(eps), k(k), Nv(Nv) {}
virtual double value(const Point<dim> &p,
[[maybe_unused]] const unsigned int component = 0) const override
{
return compute_rho(p[0], Nv);
}
private:
const double eps;
const double k;
const unsigned int Nv;
};
#endif
nufi/lsmr.h
-252
View File
@@ -1,252 +0,0 @@
#ifndef LSMR_H
#define LSMR_H
#include <cmath>
#include <limits>
#include <iomanip>
#include <iostream>
#include "nufi/blas.h"
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/nufi_solver.h
+16 -43
View File
@@ -2,33 +2,37 @@
#define NUFI_SOLVER_H
#include <boost/qvm/mat_access.hpp>
#include <vector>
#include <cmath>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/vector_tools.h>
#include <vector>
#include "nufi/fields.h" //dont remove
#include "nufi/parameters.h"
#include "nufi/poisson_problem.h"
#include "nufi/fields.h" //dont remove
using namespace dealii;
class NuFISolver
{
class NuFISolver {
public:
NuFISolver();
void run();
double eval_rho(unsigned int n, double x, const double *E_coeffs, unsigned int Nv = Parameters::NV) const;
double eval_ftilda(unsigned int n, double x, double u, const double *E_coeffs) const;
double eval_f(unsigned int n, double x, double u, const double *E_coeffs) const;
double eval_rho(unsigned int n, const double x,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
const unsigned int Nv = Parameters::NV) const;
double eval_ftilda(unsigned int n, double x, double u,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history) const;
double eval_f(unsigned int n, double x, double u,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history) const;
private:
unsigned int Nt = std::floor(Parameters::TMAX/Parameters::DT);
unsigned int Nx = Parameters::SPLINE_NX;
unsigned int Nt = std::floor(Parameters::TMAX / Parameters::DT);
unsigned int Nx = Parameters::CALC_NX;
double Lx = Parameters::LX;
@@ -40,36 +44,5 @@ private:
unsigned int order;
PoissonProblem<1> poisson;
};
template<unsigned int dim>
class ChargeDensity_NuFI : public Function<dim>
{
public:
ChargeDensity_NuFI(const double *rho_values, unsigned int Nx)
: Function<dim>(), rho(rho_values), Nx(Nx) {}
virtual double value(const Point<dim> &p,
[[maybe_unused]] const unsigned int component = 0) const override
{
const double x = p[0];
// 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:
const double *rho;
const unsigned int Nx;
};
#endif
nufi/parameters.h
+28 -30
View File
@@ -4,43 +4,41 @@
#include <cmath>
#include <cstdlib>
namespace Parameters
{
constexpr unsigned int DIMENSION = 1;
namespace Parameters {
constexpr unsigned int DIMENSION = 1;
constexpr double X_DOMAIN_LEFT = 0.0;
constexpr double X_DOMAIN_RIGHT = 4*M_PI;
constexpr double LX = std::abs(X_DOMAIN_RIGHT- X_DOMAIN_LEFT);
constexpr double LX_INV = 1/LX;
constexpr double X_DOMAIN_LEFT = 0.0;
constexpr double X_DOMAIN_RIGHT = 4 * M_PI;
constexpr double LX = std::abs(X_DOMAIN_RIGHT - X_DOMAIN_LEFT);
constexpr double LX_INV = 1 / LX;
constexpr double V_DOMAIN_LEFT = -10.;
constexpr double V_DOMAIN_RIGHT = 10.;
constexpr size_t CALC_NX = 128;
constexpr double CALC_DX = LX / CALC_NX;
constexpr unsigned int NV = 512;
constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT)/NV;
constexpr double V_DOMAIN_LEFT = -10.;
constexpr double V_DOMAIN_RIGHT = 10.;
// deal.ii options
constexpr unsigned int GLOBAL_REFINEMENT = 8;
constexpr unsigned int FE_DEGREE = 4;
constexpr unsigned int CONVERGENCE_ITERATIONS = 10000;
constexpr double CONVERGENCE_LIMIT = 1e-12;
constexpr unsigned int NV = 128;
constexpr double DV = std::abs(V_DOMAIN_RIGHT - V_DOMAIN_LEFT) / NV;
constexpr double EPS = 0.01;
constexpr double WAVE_NR = 0.5;
constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
// deal.ii options
constexpr unsigned int GLOBAL_REFINEMENT = 6;
constexpr unsigned int FE_DEGREE = 2;
constexpr unsigned int CONVERGENCE_ITERATIONS = 5000;
constexpr double CONVERGENCE_LIMIT = 1e-8;
// NUFI options
constexpr double DT=1./16.;
constexpr unsigned int TMAX = 500;
constexpr double EPS = 0.01;
constexpr double WAVE_NR = 0.5;
constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi)
//spline options
constexpr int SPLINE_NX = 256;
constexpr double SPLINE_DX = LX/(SPLINE_NX);
constexpr double SPLINE_DX_INV = 1/SPLINE_DX;
constexpr size_t SPLINE_ORDER = 4;
// NUFI options
constexpr double DT = 1. / 8.;
constexpr unsigned int TMAX = 100;
//Plotting options
constexpr int PLOT_FREQUENCY = 10;
}
// Plotting options
constexpr int PLOT_FREQUENCY = 4;
constexpr size_t PLOT_NX = CALC_NX;
constexpr double PLOT_DX = LX / PLOT_NX;
} // namespace Parameters
#endif
nufi/poisson_problem.h
+140 -151
View File
@@ -3,41 +3,45 @@
#include <deal.II/base/function.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/mpi_remote_point_evaluation.h>
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/vector_tools_evaluate.h>
#include <deal.II/numerics/vector_tools_interpolate.h>
#include <deal.II/numerics/vector_tools_point_gradient.h>
#include <deal.II/numerics/vector_tools_point_value.h>
#include <functional>
#include <memory>
#include <string>
#include <utility>
@@ -49,9 +53,7 @@ using namespace dealii;
// =-=-=-=-= Poisson Solver =-=-=-=-=
template <int dim>
class PoissonProblem
{
template <int dim> class PoissonProblem {
public:
PoissonProblem(unsigned int degree);
@@ -59,13 +61,16 @@ public:
void solve_step();
void run();
void set_rhs_function(std::unique_ptr<Function<dim>> rhs_function);
void set_rhs_function(std::function<double(const Point<dim> &)> f);
const Vector<double> &get_solution() const { return solution; }
const MappingQ<dim> &get_mapping() const { return mapping; }
const DoFHandler<dim> &get_dof_handler() const { return dof_handler; }
std::vector<double> sample_electric_field(double x_min, double x_max, unsigned int Nx);
std::vector<double> sample_electric_potential(double x_min, double x_max, unsigned int Nx);
std::vector<double> sample_electric_field(double x_min, double x_max,
unsigned int Nx);
std::vector<double> sample_electric_potential(double x_min, double x_max,
unsigned int Nx);
private:
void create_mesh();
@@ -74,18 +79,18 @@ private:
void solve();
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
AffineConstraints<double> constraints;
SparsityPattern sparsity_pattern;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
Vector<double> solution; // phi
Vector<double> solution; // phi
Vector<double> system_rhs;
std::unique_ptr<const Function<dim>> rhs_function;
std::function<double(const Point<dim> &)> rhs_function;
MappingQ<dim> mapping;
};
@@ -93,81 +98,73 @@ private:
// Utilities
template <int dim>
void PoissonProblem<dim>::set_rhs_function(std::unique_ptr<Function<dim>> rhs)
{
rhs_function = std::move(rhs);
void PoissonProblem<dim>::set_rhs_function(
std::function<double(const Point<dim> &)> f) {
rhs_function = std::move(f);
}
template <int dim>
PoissonProblem<dim>::PoissonProblem(unsigned int degree)
: fe(degree)
, dof_handler(triangulation)
, mapping(degree)
{}
: fe(degree), dof_handler(triangulation), mapping(degree) {}
template <int dim>
std::vector<double> PoissonProblem<dim>::sample_electric_field(double x_min,double x_max,unsigned int Nx)
{
std::vector<double>
PoissonProblem<dim>::sample_electric_field(double x_min, double x_max,
unsigned int Nx) {
std::vector<double> E_values(Nx);
const double dx = (x_max - x_min) / (Nx - 1);
for (unsigned int i = 0; i < Nx; ++i)
{
const double x = x_min + i * dx;
const Point<dim> point(x);
for (unsigned int i = 0; i < Nx; ++i) {
const double x = x_min + i * dx;
const Point<dim> point(x);
// 1. Find the active cell containing x
const auto cell_point_pair =
GridTools::find_active_cell_around_point(mapping,
dof_handler,
point);
// 1. Find the active cell containing x
const auto cell_point_pair =
GridTools::find_active_cell_around_point(mapping, dof_handler, point);
const auto cell = cell_point_pair.first;
const Point<dim> &unit_point = cell_point_pair.second;
const auto cell = cell_point_pair.first;
const Point<dim> &unit_point = cell_point_pair.second;
// 2. FEPointEvaluation expects an ArrayView of points
std::vector<Point<dim>> points(1, unit_point);
ArrayView<const Point<dim>> point_view(points);
// 2. FEPointEvaluation expects an ArrayView of points
std::vector<Point<dim>> points(1, unit_point);
ArrayView<const Point<dim>> point_view(points);
FEPointEvaluation<1, dim> evaluator(mapping,
dof_handler.get_fe(),
update_gradients);
FEPointEvaluation<1, dim> evaluator(mapping, dof_handler.get_fe(),
update_gradients);
// reinit with ArrayView of points
evaluator.reinit(cell, point_view);
// reinit with ArrayView of points
evaluator.reinit(cell, point_view);
Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
cell->get_dof_values(solution, local_dofs);
Vector<double> local_dofs(dof_handler.get_fe().dofs_per_cell);
cell->get_dof_values(solution, local_dofs);
// 3. Evaluate gradient at this point
evaluator.evaluate(local_dofs, EvaluationFlags::gradients);
// 3. Evaluate gradient at this point
evaluator.evaluate(local_dofs, EvaluationFlags::gradients);
const Tensor<1, dim> grad_phi = evaluator.get_gradient(0);
const Tensor<1, dim> grad_phi = evaluator.get_gradient(0);
// 4. Compute E = -grad(phi)
E_values[i] = -grad_phi[0];
// 4. Compute E = -grad(phi)
E_values[i] = -grad_phi[0];
}
return E_values;
}
template <int dim>
std::vector<double> PoissonProblem<dim>::sample_electric_potential(
double x_min,
double x_max,
unsigned int Nx)
{
std::vector<double>
PoissonProblem<dim>::sample_electric_potential(double x_min, double x_max,
unsigned int Nx) {
std::vector<double> values(Nx);
std::vector<Point<dim>> eval_points(Nx);
double Lx = x_max - x_min;
double dx = Lx / Nx;
for(unsigned int i=0 ; i<Nx; ++i)
for (unsigned int i = 0; i < Nx; ++i)
eval_points[i] = Point<1, double>(x_min + i * dx);
Utilities::MPI::RemotePointEvaluation<dim,dim> cache;
Utilities::MPI::RemotePointEvaluation<dim, dim> cache;
cache.reinit(eval_points, triangulation, mapping);
values = VectorTools::point_values<dim>(cache, dof_handler, solution);
@@ -175,33 +172,45 @@ std::vector<double> PoissonProblem<dim>::sample_electric_potential(
return values;
}
template <int dim>
double
eval_point_grad(const Mapping<dim> &mapping, const DoFHandler<dim> &dof_handler,
const Vector<double> &solution, const Point<dim> &point) {
Tensor<1, dim> grad =
VectorTools::point_gradient<dim>(mapping, dof_handler, solution, point);
double Ex = grad[0];
return Ex;
}
template <int dim>
double eval_point_value(const Mapping<dim> &mapping,
const DoFHandler<dim> &dof_handler,
const Vector<double> &solution,
const Point<dim> &point) {
return VectorTools::point_value<dim>(mapping, dof_handler, solution, point);
}
// dealii Poisson
template <int dim> void PoissonProblem<dim>::create_mesh() {
template<int dim>
void PoissonProblem<dim>::create_mesh()
{
GridGenerator::hyper_cube(triangulation,
Parameters::X_DOMAIN_LEFT,
GridGenerator::hyper_cube(triangulation, Parameters::X_DOMAIN_LEFT,
Parameters::X_DOMAIN_RIGHT);
std::vector<
GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
periodic_faces;
std::vector<GridTools::PeriodicFacePair<
typename Triangulation<dim>::cell_iterator>> periodic_faces;
GridTools::collect_periodic_faces(triangulation,
0, 1, // boundary IDs
0,
periodic_faces);
GridTools::collect_periodic_faces(triangulation, 0, 1, // boundary IDs
0, periodic_faces);
triangulation.add_periodicity(periodic_faces);
triangulation.refine_global(Parameters::GLOBAL_REFINEMENT);
}
template <int dim>
void PoissonProblem<dim>::setup_system()
{
template <int dim> void PoissonProblem<dim>::setup_system() {
dof_handler.distribute_dofs(fe);
@@ -209,10 +218,25 @@ void PoissonProblem<dim>::setup_system()
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
DoFTools::make_periodicity_constraints(dof_handler,
0, 1,
0,
constraints);
DoFTools::make_periodicity_constraints(dof_handler, 0, 1, 0, constraints);
// Gauge fix for periodic Poisson:
// remove the constant nullspace by pinning one unconstrained DoF.
// (by Paul Wilhelm)
types::global_dof_index gauge_dof = numbers::invalid_dof_index;
for (types::global_dof_index i = 0; i < dof_handler.n_dofs(); ++i) {
if (!constraints.is_constrained(i)) {
gauge_dof = i;
break;
}
}
Assert(gauge_dof != numbers::invalid_dof_index,
ExcMessage("No unconstrained DoF found for gauge fixing."));
constraints.add_line(gauge_dof);
constraints.set_inhomogeneity(gauge_dof, 0.0);
constraints.close();
@@ -226,112 +250,77 @@ void PoissonProblem<dim>::setup_system()
system_rhs.reinit(dof_handler.n_dofs());
}
template <int dim>
void PoissonProblem<dim>::assemble_system()
{
QGauss<dim> quadrature_formula(fe.degree + 1);
// Paul
template <int dim> void PoissonProblem<dim>::assemble_system() {
Assert(system_matrix.m() == dof_handler.n_dofs(),
ExcMessage("Matrix not initialized correctly"));
system_matrix = 0;
system_rhs = 0;
QGauss<dim> quadrature_formula(fe.degree + 1);
FEValues<dim> fe_values(fe, quadrature_formula,
update_values |
update_gradients |
update_quadrature_points |
update_JxW_values);
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
Assert(rhs_function != nullptr, ExcMessage("RHS function not set"));
for (const auto &cell : dof_handler.active_cell_iterators())
{
for (const auto &cell : dof_handler.active_cell_iterators()) {
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
cell_rhs = 0;
for (const auto q : fe_values.quadrature_point_indices())
{
const double rho = rhs_function->value(fe_values.quadrature_point(q));
for (const auto q : fe_values.quadrature_point_indices()) {
const double rho = rhs_function(
fe_values.quadrature_point(q)); // Eval rhs_function at q points
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q) * // grad phi_i(x_q)
fe_values.shape_grad(j, q) * // grad phi_j(x_q)
fe_values.JxW(q)); // dx
cell_matrix(i, j) += fe_values.shape_grad(i, q) *
fe_values.shape_grad(j, q) * fe_values.JxW(q);
for (const unsigned int i : fe_values.dof_indices())
cell_rhs(i) += (fe_values.shape_value(i, q) * // phi_i(x_q)
rho * // f(x_q)
fe_values.JxW(q)); // dx
cell_rhs(i) += fe_values.shape_value(i, q) * rho * fe_values.JxW(q);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(cell_matrix,
cell_rhs,
local_dof_indices,
system_matrix,
system_rhs);
for (const unsigned int i : fe_values.dof_indices())
for (const unsigned int j : fe_values.dof_indices())
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
for (const unsigned int i : fe_values.dof_indices())
system_rhs(local_dof_indices[i]) += cell_rhs(i);
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
std::map<types::global_dof_index, double> boundary_values;
// VectorTools::interpolate_boundary_values(dof_handler,
// types::boundary_id(0),
// Functions::ZeroFunction<1>(),
// boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
template <int dim> void PoissonProblem<dim>::solve() {
template <int dim>
void PoissonProblem<dim>::solve()
{
SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS, Parameters::CONVERGENCE_LIMIT);
SolverControl solver_control(Parameters::CONVERGENCE_ITERATIONS,
Parameters::CONVERGENCE_LIMIT);
SolverCG<Vector<double>> solver(solver_control);
// PreconditionSSOR<SparseMatrix<double>> preconditioner;
// preconditioner.initialize(system_matrix, 1.2);
// solver.solve(system_matrix, solution, system_rhs, preconditioner);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
// constraints.distribute(solution);
constraints.distribute(solution);
}
template <int dim>
void PoissonProblem<dim>::initialize()
{
create_mesh(); // build grid
setup_system(); // distribute DoFs and matrices
template <int dim> void PoissonProblem<dim>::initialize() {
create_mesh(); // build grid
setup_system(); // distribute DoFs and matrices
}
template <int dim>
void PoissonProblem<dim>::solve_step()
{
template <int dim> void PoissonProblem<dim>::solve_step() {
assemble_system();
solve();
}
// NuFI doesnt use this, kept only for testing PoissonProblem
template <int dim>
void PoissonProblem<dim>::run()
{
template <int dim> void PoissonProblem<dim>::run() {
create_mesh();
setup_system();
assemble_system();
nufi/save_results.h
+16 -18
View File
@@ -1,28 +1,26 @@
#ifndef SAVE_RESULTS_H
#define SAVE_RESULTS_H
#include <string>
#include "nufi/nufi_solver.h"
#include "nufi/poisson_problem.h"
#include <deal.II/lac/vector.h>
#include <string>
void save_f(const NuFISolver &solver, unsigned int n,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history, unsigned int Nx_out,
unsigned int Nv_out, const std::string &filename);
void save_f( const NuFISolver &solver,
unsigned int n,
const double *E_coeffs,
unsigned int Nx_out,
unsigned int Nv_out,
const std::string &filename);
void save_rho(const NuFISolver &solver, unsigned int n,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
unsigned int Nx_out, const std::string &filename);
void save_rho(const NuFISolver &solver,
unsigned int n,
const double *E_coeffs,
unsigned int Nx_out,
const std::string &filename);
void save_Efield(unsigned int n, const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
unsigned int Nx_out, const std::string &filename);
void save_Efield(unsigned int n,
const double *E_coeffs,
unsigned int Nx_out,
const std::string &filename);
void save_space_vector(const std::vector<double>& vals, const std::string& filename, size_t it);
void save_space_vector(const std::vector<double> &vals,
const std::string &filename, size_t it);
#endif
nufi/splines.h
-90
View File
@@ -1,90 +0,0 @@
#ifndef SPLINES_HP
#define SPLINES_HP
#include <cstddef>
namespace splines1d
{
template <typename real>
constexpr real faculty( size_t n ) noexcept
{
return (n > 1) ? real(n)*faculty<real>(n-1) : real(1);
}
template <typename real, size_t order, size_t derivative = 0>
void N( real x, real *result, size_t stride = 1 ) noexcept
{
static_assert( order > 0, "Splines must have order greater than zero." );
constexpr int n { order };
constexpr int d { derivative };
if ( derivative >= order )
for ( size_t i = 0; i < order; ++i )
result[ i*stride ] = 0;
if ( n == 1 )
{
*result = 1;
return;
}
real v[n]; v[n-1] = 1;
for ( int k = 1; k < n - d; ++k )
{
v[n-k-1] = (1-x)*v[n-k];
for ( int i = 1-k; i < 0; ++i )
v[n-1+i] = (x-i)*v[n-1+i] + (k+1+i-x)*v[n+i];
v[n-1] *= x;
}
// Differentiate if necessary.
for ( size_t j = derivative; j-- > 0; )
{
v[j] = -v[j+1];
for ( size_t i = j + 1; i < order - 1; ++i )
v[i] = v[i] - v[i+1];
}
constexpr real factor = real(1) / faculty<real>(order-derivative-1);
for ( size_t i = 0; i < order; ++i )
result[i*stride] = v[i]*factor;
}
template <typename real, size_t order, size_t derivative = 0>
real eval( real x, const real *coefficients, size_t stride = 1 ) noexcept
{
static_assert( order > 0, "Splines must have order greater than zero." );
static_assert( order > derivative, "Too high derivative requested." );
constexpr size_t n { order };
constexpr size_t d { derivative };
if ( d >= n ) return 0;
if ( n == 1 ) return *coefficients;
// Gather coefficients.
real c[ order ];
for ( size_t j = 0; j < order; ++j )
c[j] = coefficients[ stride * j ];
// Differentiate if necessary.
for ( size_t j = 1; j <= d; ++j )
for ( size_t i = n; i-- > j; )
c[i] = c[i] - c[i-1];
// Evaluate using de Boor’s algorithm.
for ( size_t j = 1; j < n-d; ++j )
for ( size_t i = n-d; i-- > j; )
c[d+i] = (x+n-d-1-i)*c[d+i] + (i-j+1-x)*c[d+i-1];
constexpr real factor = real(1) / faculty<real>(order-derivative-1);
return factor*c[n-1];
}
}
#endif
nufi/stopwatch.h
+12 -20
View File
@@ -3,35 +3,27 @@
#include <chrono>
template <typename real>
class stopwatch
{
template <typename real> class stopwatch {
public:
void reset();
real elapsed();
void reset();
real elapsed();
private:
using clock = std::chrono::high_resolution_clock;
clock::time_point t0 { clock::now() };
using clock = std::chrono::high_resolution_clock;
clock::time_point t0{clock::now()};
};
template <typename real> inline
void stopwatch<real>::reset()
{
t0 = clock::now();
template <typename real> inline void stopwatch<real>::reset() {
t0 = clock::now();
}
template <typename real> inline
real stopwatch<real>::elapsed()
{
using seconds = std::chrono::duration<real,std::ratio<1,1>>;
template <typename real> inline real stopwatch<real>::elapsed() {
using seconds = std::chrono::duration<real, std::ratio<1, 1>>;
auto tnow = clock::now();
auto duration = std::chrono::duration_cast<seconds>( tnow - t0 );
auto tnow = clock::now();
auto duration = std::chrono::duration_cast<seconds>(tnow - t0);
return duration.count();
return duration.count();
}
#endif // STOPWATCH_H
rho_E.mp4
BIN
View File
Binary file not shown.
src/blas.cc
-95
View File
@@ -1,95 +0,0 @@
#include "nufi/blas.h"
#include <cblas.h>
namespace blas
{
double dot( const size_t n, const double *x, size_t incx,
const double *y, size_t incy )
{
return cblas_ddot(n,x,incx,y,incy);
}
float dot( const size_t n, const float *x, size_t incx,
const float *y, size_t incy )
{
return cblas_sdot(n,x,incx,y,incy);
}
void axpy( size_t n, double alpha, const double *x, size_t incx,
double *y, size_t incy )
{
cblas_daxpy(n,alpha,x,incx,y,incy);
}
void axpy( size_t n, float alpha, const float *x, size_t incx,
float *y, size_t incy )
{
cblas_saxpy(n,alpha,x,incx,y,incy);
}
void scal( size_t n, double alpha, double *x, size_t incx )
{
cblas_dscal(n,alpha,x,incx);
}
void scal( size_t n, float alpha, float *x, size_t incx )
{
cblas_sscal(n,alpha,x,incx);
}
void copy( size_t n, const double *x, size_t incx, double *y, size_t incy )
{
cblas_dcopy(n,x,incx,y,incy);
}
void copy( size_t n, const float *x, size_t incx, float *y, size_t incy )
{
cblas_scopy(n,x,incx,y,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)
{
cblas_dger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, 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)
{
cblas_sger( CblasColMajor, M, N, alpha, X, incX, Y, incY, A, 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 )
{
if ( trans == 'T' || trans == 'Y' )
{
cblas_dgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
}
else
{
cblas_dgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, 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 )
{
if ( trans == 'T' || trans == 'Y' )
{
cblas_sgemv( CblasColMajor, CblasTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
}
else
{
cblas_sgemv( CblasColMajor, CblasNoTrans, m, n, alpha, a, lda, x, incx, beta, y, incy );
}
}
}
src/main.cc
+25 -34
View File
@@ -1,44 +1,35 @@
#include <iostream>
#include <filesystem>
#include <iostream>
#include <nufi/nufi_solver.h>
#include <omp.h>
void clear_results_directory(const std::string &dir)
{
if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
return;
void clear_results_directory(const std::string &dir) {
if (!std::filesystem::exists(dir) || !std::filesystem::is_directory(dir))
return;
for (const auto &entry : std::filesystem::directory_iterator(dir))
{
if (std::filesystem::is_regular_file(entry))
{
std::filesystem::remove(entry.path());
std::cout << "Deleted: " << entry.path() << '\n';
}
for (const auto &entry : std::filesystem::directory_iterator(dir)) {
if (std::filesystem::is_regular_file(entry)) {
std::filesystem::remove(entry.path());
std::cout << "Deleted: " << entry.path() << '\n';
}
}
}
int main()
{
#include <omp.h>
std::cout << "Threads: " << omp_get_max_threads() << "\n";
try
{
clear_results_directory("results");
int main() {
std::cout << "Threads: " << omp_get_max_threads() << "\n";
try {
clear_results_directory("results");
NuFISolver solver;
solver.run();
}
catch (const std::exception &exc)
{
std::cerr << "\nException:\n" << exc.what() << "\n";
return 1;
}
catch (...)
{
std::cerr << "\nUnknown exception!\n";
return 1;
}
NuFISolver solver;
solver.run();
return 0;
} catch (const std::exception &exc) {
std::cerr << "\nException:\n" << exc.what() << "\n";
return 1;
} catch (...) {
std::cerr << "\nUnknown exception!\n";
return 1;
}
return 0;
}
src/nufi_solver.cc
+110 -125
View File
@@ -1,207 +1,192 @@
#include "nufi/nufi_solver.h"
#include <algorithm>
#include <boost/qvm/mat_access.hpp>
#include <cmath>
#include <cstdlib>
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/vector_tools.h>
#include <cstddef>
#include <iostream>
#include <memory>
#include <ostream>
#include <cstddef>
#include <vector>
#include "nufi/parameters.h"
#include "nufi/save_results.h"
#include "nufi/poisson_problem.h"
#include "nufi/fields.h"
#include "nufi/parameters.h"
#include "nufi/poisson_problem.h"
#include "nufi/save_results.h"
#include "nufi/stopwatch.h"
using namespace dealii;
double NuFISolver::eval_ftilda(unsigned int n,
double x,
double u,
const double *E_coeffs) const
{
if ( n == 0 ) return f0(x,u);
const size_t order = Parameters::SPLINE_ORDER;
const size_t stride_x = 1;
const size_t stride_t = stride_x*(Nx + order - 1);
double
NuFISolver::eval_ftilda(unsigned int n, double x, double u,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history) const {
if (n == 0)
return f0(x, u);
double Ex;
const double *c;
// We omit the initial half-step.
while ( --n )
{
x = x - Parameters::DT *u;
c = E_coeffs + n*stride_t;
Ex = -eval<1>(x, c);
u = u + Parameters::DT *Ex;
while (--n) {
x = x - Parameters::DT * u;
Ex = -eval(x, poisson, phi_history[n]);
u = u + Parameters::DT * Ex;
}
// The final half-step.
x -= Parameters::DT*u;
c = E_coeffs + n*stride_t;
Ex = -eval<1>(x, c);
u += 0.5*Parameters::DT*Ex;
x = x - Parameters::DT * u;
Ex = -eval(x, poisson, phi_history[n]);
u += 0.5 * Parameters::DT * Ex;
return f0(x,u);
return f0(x, u);
}
double NuFISolver::eval_f(unsigned int n,
double x,
double u,
const double *E_coeffs) const
{
if ( n == 0 ) return f0(x,u);
const size_t order = Parameters::SPLINE_ORDER;
const size_t stride_x = 1;
const size_t stride_t = stride_x*(Nx + order - 1);
double
NuFISolver::eval_f(unsigned int n, double x, double u,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history) const {
if (n == 0)
return f0(x, u);
double Ex;
const double *c;
// Initial half-step.
c = E_coeffs + n*stride_t;
Ex = -eval<1>(x, c);
u += 0.5*Parameters::DT * Ex;
Ex = -eval(x, poisson, phi_history[n]);
u += 0.5 * Parameters::DT * Ex;
while ( --n )
{
x = x - Parameters::DT *u;
c = E_coeffs + n*stride_t;
Ex = -eval<1>(x, c);
u = u + Parameters::DT *Ex;
while (--n) {
x = x - Parameters::DT * u;
Ex = -eval(x, poisson, phi_history[n]);
u = u + Parameters::DT * Ex;
}
// The final half-step.
x -= Parameters::DT*u;
c = E_coeffs + n*stride_t;
Ex = -eval<1>(x, c);
u += 0.5*Parameters::DT*Ex;
x = x - Parameters::DT * u;
Ex = -eval(x, poisson, phi_history[n]);
u += 0.5 * Parameters::DT * Ex;
return f0(x,u);
return f0(x, u);
}
double NuFISolver::eval_rho(const unsigned int n,
const double x,
const double *E_coeffs,
const unsigned int Nv) const
{
const double dv = (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
const double v_min = Parameters::V_DOMAIN_LEFT;
double NuFISolver::eval_rho(unsigned int n, const double x,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
const unsigned int Nv) const {
const double dv =
(Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;
const double v_min = Parameters::V_DOMAIN_LEFT + 0.5 * dv;
double integral = 0.0;
#pragma omp parallel for reduction (+ : integral)
#pragma omp parallel for reduction(+ : integral)
for (unsigned int i = 0; i < Nv; ++i)
integral += eval_ftilda(n, x, v_min + i * dv, E_coeffs);
return 1.0 - integral*dv;
integral += eval_ftilda(n, x, v_min + i * dv, poisson, phi_history);
return 1.0 - integral * dv;
}
void NuFISolver::run()
{
void NuFISolver::run() {
std::cout << "Building E_sline\n\n";
using std::abs;
using std::max;
const size_t stride_t = Nx + order - 1;
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 };
std::unique_ptr<double, decltype(std::free) *> rho{
reinterpret_cast<double *>(std::aligned_alloc(64, sizeof(double) * Nx)),
std::free};
std::vector<double> int_E_squared;
int_E_squared.reserve(Nt);
if ( rho == nullptr ) throw std::bad_alloc {};
std::vector<Vector<double>> phi_history;
Gradient grad(x_min, x_max, Nx);
if (rho == nullptr)
throw std::bad_alloc{};
double total_time = 0;
std::ofstream time_file("results/simulation_time.txt");
time_file << "# it step_time total_time" << "\n";
for (unsigned int it = 0; it < Nt; ++it)
{
stopwatch<double> timer;
const double x_min = Parameters::X_DOMAIN_LEFT;
double dx = Parameters::CALC_DX;
double time_elapsed_before = timer.elapsed();
for (unsigned int it = 0; it < Nt; ++it) {
stopwatch<double> timer;
std::cout << "Timestep " << it << " / " << Nt << " (simulation time = "<< it*Parameters::DT << ")"<< std::endl;
double time_elapsed_before = timer.elapsed();
// compute rho
std::cout << "Timestep " << it << " / " << Nt
<< " (simulation time = " << it * Parameters::DT << ")"
<< std::endl;
double dx = Parameters::SPLINE_DX;
// compute rho
#pragma omp parallel for
for(size_t i = 0; i<Nx; i++)
{
double x = Parameters::X_DOMAIN_LEFT + i*dx;
double ith_rho = eval_rho(it, x, coeffs.get(),Parameters::NV);
#pragma omp parallel for
for (size_t i = 0; i < Nx; i++) {
double x = Parameters::X_DOMAIN_LEFT + i * dx;
double ith_rho = eval_rho(it, x, poisson, phi_history, Parameters::NV);
AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
rho.get()[i] = ith_rho;
}
AssertThrow(std::isfinite(ith_rho), ExcMessage("NaN detected in rho"));
rho.get()[i] = ith_rho;
}
poisson.set_rhs_function(std::make_unique<ChargeDensity_NuFI<1>>(rho.get(), Nx));
poisson.solve_step();
poisson.set_rhs_function([&rho, x_min, dx, Nx = Nx](const Point<1> &p) {
double x = p[0];
int i = static_cast<int>(std::floor((x - x_min) / dx));
i = (i % Nx + Nx) % Nx;
return rho.get()[i];
});
std::vector<double> sampled_potential = poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of FE
poisson.solve_step();
// These have been tested to be equivalent
// ////////////////////////////////////////////////
// std::vector<double> E_vals = grad.compute(sampled_potential); // vector grad of FE solution
// save_space_vector(E_vals, "electric", it);
// std::vector<double> E_vals_deal = poisson.sample_electric_field(x_min, x_max, Nx); // FE grad of soution
// save_space_vector(E_vals_deal, "electricdeal", it);
// ////////////////////////////////////////////////
phi_history.push_back(poisson.get_solution());
// std::vector<double> sampled_potential =
// poisson.sample_electric_potential(x_min, x_max, Nx); // Solution of
// FE
// interpolate and save current field
double* current_coeffs = coeffs.get() + it*stride_t;
interpolate<double, Parameters::SPLINE_ORDER>(current_coeffs, sampled_potential.data());
double timer_elapsed = timer.elapsed();
double step_time = timer_elapsed - time_elapsed_before;
total_time += timer_elapsed;
time_file << it << " " << step_time << " " << total_time << "\n";
time_file.flush();
std::vector<double> E_x(Nx,0.0) ;
#pragma omp parallel for
for(size_t ix=0; ix<Nx; ++ix)
{
E_x[ix] = -eval<1>(Parameters::X_DOMAIN_LEFT+ix*dx, current_coeffs);
std::cout << "step made in " << step_time << " seconds\n\n";
if (it % Parameters::PLOT_FREQUENCY == 0) {
std::cout << "Saving results... ";
save_f(*this, it, poisson, phi_history, Parameters::PLOT_NX,
Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
save_rho(*this, it, poisson, phi_history, Parameters::PLOT_NX,
"results/rho_" + std::to_string(it) + ".dat");
// save_Efield(it, coeffs.get(), 128, "results/field_" +
// std::to_string(it) + ".dat");
std::vector<double> E_x(Nx, 0.0);
#pragma omp parallel for
for (size_t ix = 0; ix < Nx; ++ix) {
E_x[ix] = -eval(x_min + ix * dx, poisson, phi_history[it]);
}
double timer_elapsed = timer.elapsed();
total_time += timer_elapsed;
save_space_vector(E_x, "field", it);
std::cout << "step made in "<< timer_elapsed-time_elapsed_before <<" seconds\n\n";
if (it % Parameters::PLOT_FREQUENCY == 0)
{
std::cout << "Saving results... ";
save_f(*this, it, coeffs.get(), Parameters::SPLINE_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat");
save_rho(*this, it, coeffs.get(), Parameters::SPLINE_NX, "results/rho_" + std::to_string(it) + ".dat");
// save_Efield(it, coeffs.get(), 128, "results/field_" + std::to_string(it) + ".dat");
save_space_vector(E_x, "field", it);
double int_val = 0.5 * integral_space_vector_squared(current_coeffs);
int_E_squared.push_back(int_val);
save_space_vector(int_E_squared, "electricint", it);
std::cout << "Time since start = "<< total_time<<"\n\n";
}
double int_val =
0.5 * integral_space_vector_squared(poisson, phi_history[it]);
int_E_squared.push_back(int_val);
save_space_vector(int_E_squared, "electricint", it);
std::cout << "Time since start = " << total_time << "\n\n";
}
}
std::cout << "NuFI simulation finished in "<< total_time <<" seconds.\n";
std::cout << "NuFI simulation finished in " << total_time << " seconds.\n";
}
NuFISolver::NuFISolver()
: order(Parameters::FE_DEGREE),
poisson(order)
{
NuFISolver::NuFISolver() : order(Parameters::FE_DEGREE), poisson(order) {
std::cout << "Initializing dealii Poisson Solver\n";
poisson.initialize();
}
src/save_results.cc
+44 -54
View File
@@ -1,20 +1,18 @@
#include "nufi/save_results.h"
#include "nufi/fields.h"
#include "nufi/nufi_solver.h"
#include "nufi/parameters.h"
#include "nufi/poisson_problem.h"
#include <fstream>
#include <stdexcept>
#include <string>
#include <vector>
#include "nufi/nufi_solver.h"
void save_f( const NuFISolver &solver,
unsigned int n,
const double *E_coeffs,
unsigned int Nx_out,
unsigned int Nv_out,
const std::string &filename)
{
void save_f(const NuFISolver &solver, unsigned int n,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history, unsigned int Nx_out,
unsigned int Nv_out, const std::string &filename) {
std::ofstream file(filename);
double xmin = Parameters::X_DOMAIN_LEFT;
@@ -30,34 +28,30 @@ void save_f( const NuFISolver &solver,
file << xmin << " " << xmax << "\n";
file << vmin << " " << vmax << "\n";
for (unsigned int i = 0; i < Nx_out; ++i)
{
double x = xmin + (i + 0.5)*dx;
for (unsigned int i = 0; i < Nx_out; ++i) {
double x = xmin + (i + 0.5) * dx;
for (unsigned int j = 0; j < Nv_out; ++j)
{
double v = vmin + (j + 0.5)*dv;
for (unsigned int j = 0; j < Nv_out; ++j) {
double v = vmin + (j + 0.5) * dv;
double val = solver.eval_f(n, x, v, E_coeffs);
double val = solver.eval_f(n, x, v, poisson, phi_history);
file << val;
file << val;
if (j < Nv_out - 1)
file << " ";
}
if (j < Nv_out - 1)
file << " ";
}
file << "\n";
file << "\n";
}
file.close();
}
void save_rho(const NuFISolver &solver,
unsigned int n,
const double *E_coeffs,
unsigned int Nx_out,
const std::string &filename)
{
void save_rho(const NuFISolver &solver, unsigned int n,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
unsigned int Nx_out, const std::string &filename) {
std::ofstream file(filename);
double xmin = Parameters::X_DOMAIN_LEFT;
@@ -67,20 +61,18 @@ void save_rho(const NuFISolver &solver,
file << Nx_out << "\n";
file << xmin << " " << xmax << "\n";
for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx)
{
double val = solver.eval_rho(n, xmin, E_coeffs);
file << val;
file << "\n";
for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) {
double val = solver.eval_rho(n, xmin, poisson, phi_history);
file << val;
file << "\n";
}
file.close();
}
void save_Efield(unsigned int n,
const double *E_coeffs,
unsigned int Nx_out,
const std::string &filename)
{
void save_Efield([[maybe_unused]] unsigned int n,
const PoissonProblem<1> &poisson,
const std::vector<Vector<double>> &phi_history,
unsigned int Nx_out, const std::string &filename) {
std::ofstream file(filename);
double xmin = Parameters::X_DOMAIN_LEFT;
@@ -88,31 +80,29 @@ void save_Efield(unsigned int n,
double dx = (xmax - xmin) / Nx_out;
// select from E_coeffs
const size_t stride_x = 1;
const size_t stride_t = stride_x*(Parameters::SPLINE_NX + Parameters::SPLINE_ORDER - 1);
const double *c;
c = E_coeffs + n*stride_t;
file << Nx_out << "\n";
file << xmin << " " << xmax << "\n";
for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx)
{
double val = -eval<1>(xmin, c);
file << val;
file << "\n";
for (unsigned int i = 0; i < Nx_out; ++i, xmin += dx) {
double val = -eval(xmin, poisson, phi_history[n]);
file << val;
file << "\n";
}
file.close();
}
void save_space_vector(const std::vector<double>& vals, const std::string& filename, size_t it)
{
std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat");
void save_space_vector(const std::vector<double> &vals,
const std::string &filename, size_t it) {
std::ofstream file("results/" + filename + "_" + std::to_string(it) + ".dat");
if (!file) throw std::runtime_error("failed to start file in results/");
if (!file)
throw std::runtime_error("failed to start file in results/");
file << vals.size() << "\n";
file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT << "\n";
file << std::fixed << std::setprecision(8);
for (double val : vals) file << val << "\n";
file << vals.size() << "\n";
file << Parameters::X_DOMAIN_LEFT << " " << Parameters::X_DOMAIN_RIGHT
<< "\n";
file << std::fixed << std::setprecision(8);
for (double val : vals)
file << val << "\n";
}