From 8be3dadaa3a784b712f3ea225eddd47483fc12f7 Mon Sep 17 00:00:00 2001 From: "lisa.pizzo" Date: Sat, 10 Jan 2026 15:28:21 +0100 Subject: [PATCH] . --- Sheet7/E14/jacob_template/userset.h | 47 +++++ Sheet7/E14/jacob_template/vdop.cpp | 122 +++++++++++++ Sheet7/E14/jacob_template/vdop.h | 167 ++++++++++++++++++ .../jacob_template/visualize_par_results.m | 52 ++++++ Sheet7/E14/jacob_template/visualize_results.m | 20 +++ 5 files changed, 408 insertions(+) create mode 100644 Sheet7/E14/jacob_template/userset.h create mode 100644 Sheet7/E14/jacob_template/vdop.cpp create mode 100644 Sheet7/E14/jacob_template/vdop.h create mode 100644 Sheet7/E14/jacob_template/visualize_par_results.m create mode 100644 Sheet7/E14/jacob_template/visualize_results.m diff --git a/Sheet7/E14/jacob_template/userset.h b/Sheet7/E14/jacob_template/userset.h new file mode 100644 index 0000000..94a10b8 --- /dev/null +++ b/Sheet7/E14/jacob_template/userset.h @@ -0,0 +1,47 @@ +#ifndef USERSET_FILE +#define USERSET_FILE +#include +/** + * User function: f(@p x,@p y) + * @param[in] x x-coordinate of discretization point + * @param[in] y y-coordinate of discretization point + * @return value for right hand side f(@p x,@p y) + */ +double FunctF(double const x, double const y); + +/** + * User function: u(@p x,@p y) + * @param[in] x x-coordinate of discretization point + * @param[in] y y-coordinate of discretization point + * @return value for solution vector u(@p x,@p y) + */ +double FunctU(double const x, double const y); + + +/** + * User function: f(@p x,@p y) = @f$ x^2 \sin(2.5\pi y)@f$. + * @param[in] x x-coordinate of discretization point + * @param[in] y y-coordinate of discretization point + * @return value f(@p x,@p y) + */ +inline +double fNice(double const x, double const y) +{ + //return x * x * std::sin(2.5 * M_PI * y); // solution u + return std::sin(M_PI*2.5*y)*(M_PI*M_PI*2.5*2.5*x*x - 2); // -Laplacian(u) +} + +/** + * User function: f(@p x,@p y) = 0$. + * @param[in] x x-coordinate of discretization point + * @param[in] y y-coordinate of discretization point + * @return value 0 + */ +inline +double f_zero(double const x, double const y) +//double f_zero(double const /*x*/, double const /*y*/) +{ + return 0.0 + 0.0*(x+y); +} + +#endif diff --git a/Sheet7/E14/jacob_template/vdop.cpp b/Sheet7/E14/jacob_template/vdop.cpp new file mode 100644 index 0000000..a8b31e0 --- /dev/null +++ b/Sheet7/E14/jacob_template/vdop.cpp @@ -0,0 +1,122 @@ +#include "vdop.h" +#include // assert() +#include +#include +#include +using namespace std; + + +void vddiv(vector & x, vector const& y, + vector const& z) +{ + assert( x.size()==y.size() && y.size()==z.size() ); + size_t n = x.size(); + +#pragma omp parallel for + for (size_t k = 0; k < n; ++k) + { + x[k] = y[k] / z[k]; + } + return; +} + +//****************************************************************************** + +void vdaxpy(std::vector & x, std::vector const& y, + double alpha, std::vector const& z ) +{ + assert( x.size()==y.size() && y.size()==z.size() ); + size_t n = x.size(); + +#pragma omp parallel for + for (size_t k = 0; k < n; ++k) + { + x[k] = y[k] + alpha * z[k]; + } + return; +} +//****************************************************************************** + +double dscapr(std::vector const& x, std::vector const& y) +{ + assert( x.size()==y.size()); + size_t n = x.size(); + + double s = 0.0; +#pragma omp parallel for reduction(+:s) + for (size_t k = 0; k < n; ++k) + { + s += x[k] * y[k]; + } + + return s; +} + +//****************************************************************************** +bool CompareVectors(std::vector const& x, int const n, double const y[], double const eps) +{ + bool bn = (static_cast(x.size())==n); + if (!bn) + { + cout << "######### Error: " << "number of elements" << endl; + } + //bool bv = equal(x.cbegin(),x.cend(),y); + bool bv = equal(x.cbegin(),x.cend(),y, + [eps](double a, double b) -> bool + { return std::abs(a-b)(x.size())==n); + cout << "######### Error: " << "values" << endl; + } + return bn && bv; +} + +//****************************************************************************** +double par_scalar(vector const &x, vector const &y, MPI_Comm const& icomm) +{ + const double s = dscapr(x,y); + double sg; + MPI_Allreduce(&s,&sg,1,MPI_DOUBLE,MPI_SUM,icomm); + + return(sg); +} + +//****************************************************************************** +void ExchangeAll(vector const &xin, vector &yout, MPI_Comm const &icomm) +{ + int myrank, numprocs,ierr(-1); + MPI_Comm_rank(icomm, &myrank); // my MPI-rank + MPI_Comm_size(icomm, &numprocs); + int const N=xin.size(); + int const sendcount = N/numprocs; // equal sized junks + assert(sendcount*numprocs==N); // really all junk sized? + assert(xin.size()==yout.size()); + + auto sendbuf = xin.data(); + auto recvbuf = yout.data(); + ierr = MPI_Alltoall(sendbuf, sendcount, MPI_DOUBLE, + recvbuf, sendcount, MPI_DOUBLE, icomm); + assert(0==ierr); + + return; +} + +//****************************************************************************** +void ExchangeAllInPlace(vector &xin, MPI_Comm const &icomm) +{ + int myrank, numprocs,ierr(-1); + MPI_Comm_rank(icomm, &myrank); // my MPI-rank + MPI_Comm_size(icomm, &numprocs); + int const N=xin.size(); + int const sendcount = N/numprocs; // equal sized junks + assert(sendcount*numprocs==N); // really all junk sized? + + auto sendbuf = xin.data(); + ierr = MPI_Alltoall(MPI_IN_PLACE, sendcount, MPI_DOUBLE, + sendbuf, sendcount, MPI_DOUBLE, icomm); + assert(0==ierr); + + return; +} diff --git a/Sheet7/E14/jacob_template/vdop.h b/Sheet7/E14/jacob_template/vdop.h new file mode 100644 index 0000000..db6834b --- /dev/null +++ b/Sheet7/E14/jacob_template/vdop.h @@ -0,0 +1,167 @@ +#ifndef VDOP_FILE +#define VDOP_FILE +#include +#include // MPI +#include +#include + +/** @brief Element-wise vector divison x_k = y_k/z_k. + * + * @param[out] x target vector + * @param[in] y source vector + * @param[in] z source vector + * + */ +void vddiv(std::vector &x, std::vector const &y, + std::vector const &z); + +/** @brief Element-wise daxpy operation x(k) = y(k) + alpha*z(k). + * + * @param[out] x target vector + * @param[in] y source vector + * @param[in] alpha scalar + * @param[in] z source vector + * + */ +void vdaxpy(std::vector &x, std::vector const &y, + double alpha, std::vector const &z ); + + +/** @brief Calculates the Euclidean inner product of two vectors. + * + * @param[in] x vector + * @param[in] y vector + * @return Euclidean inner product @f$\langle x,y \rangle@f$ + * + */ +double dscapr(std::vector const &x, std::vector const &y); + + +inline +double L2_scapr(std::vector const &x, std::vector const &y) +{ + return dscapr(x, y) / x.size(); +} + + +/** Parallel inner product + @param[in] x vector + @param[in] y vector + @param[in] icomm MPI communicator + @return resulting Euclidian inner product +*/ +double par_scalar(std::vector const &x, std::vector const &y, + MPI_Comm const& icomm=MPI_COMM_WORLD); + + + +/* ReadId : Input and broadcast of an integer */ +inline +int ReadIn(std::string const &ss = std::string(), MPI_Comm const &icomm = MPI_COMM_WORLD) +{ + MPI_Barrier(icomm); + int myrank; /* my rank number */ + MPI_Comm_rank(icomm, &myrank); + int id; + + if (myrank == 0) { + std::cout << "\n\n " << ss << " : Which process do you want to debug ? \n"; + std::cin >> id; + } + MPI_Bcast(&id, 1, MPI_INT, 0, icomm); + + return id; +} + +/** + * Print entries of a vector to standard output. + * + * @param[in] v vector values + * @param[in] ss string containing the vector name + * @param[in] icomm communicator group for MPI + * +*/ +template +void DebugVector(std::vector const &v, std::string const &ss = std::string(), MPI_Comm const &icomm = MPI_COMM_WORLD) +{ + MPI_Barrier(icomm); + std::cout.flush(); + int numprocs; /* # processes */ + MPI_Comm_size(icomm, &numprocs); + int myrank; /* my rank number */ + MPI_Comm_rank(icomm, &myrank); + + int readid = ReadIn(ss); /* Read readid */ + + while ( (0 <= readid) && (readid < numprocs) ) { + if (myrank == readid) { + std::cout << "\n\n process " << readid; + std::cout << "\n .... " << ss << " (nnode = " << v.size() << ")\n"; + for (size_t j = 0; j < v.size(); ++j) { + std::cout.setf(std::ios::right, std::ios::adjustfield); + std::cout << "[" << j << "] "<< v[j] << " "; + } + std::cout << std::endl; + std::cout.flush(); + } + + readid = ReadIn(ss, icomm); /* Read readid */ + } + MPI_Barrier(icomm); + return; +} + +/** @brief Compares an STL vector with POD vector. + * + * The accuracy criteria @f$ |x_k-y_k| < \varepsilon \left({1+0.5(|x_k|+|y_k|)}\right) @f$ + * follows the book by + * Stoyan/Baran, p.8. + * + * @param[in] x STL vector + * @param[in] n length of POD vector + * @param[in] y POD vector + * @param[in] eps relative accuracy criteria (default := 0.0). + * @return true iff pairwise vector elements are relatively close to each other. + * + */ +bool CompareVectors(std::vector const &x, int n, double const y[], double const eps = 0.0); + + +/** Output operator for vector + * @param[in,out] s output stream, e.g. @p cout + * @param[in] v vector + * + * @return output stream +*/ +template +std::ostream& operator<<(std::ostream &s, std::vector const &v) +{ + for (auto vp: v) + { + s << vp << " "; + } + return s; +} + +/** Exchanges equal size partions of vector @p xin with all MPI processes. + * The received data are return in vector @p yout . + * + * @param[in] xin input vector + * @param[out] yout output vector + * @param[in] icomm MPI communicator + * +*/ +void ExchangeAll(std::vector const &xin, std::vector &yout, MPI_Comm const &icomm = MPI_COMM_WORLD); + +/** Exchanges equal size partions of vector @p xin with all MPI processes. + * The received data are return in vector @p xin . + * + * @param[in,out] xin input/output vector + * @param[in] icomm MPI communicator + * +*/ +void ExchangeAllInPlace(std::vector &xin, MPI_Comm const &icomm = MPI_COMM_WORLD); + + + +#endif diff --git a/Sheet7/E14/jacob_template/visualize_par_results.m b/Sheet7/E14/jacob_template/visualize_par_results.m new file mode 100644 index 0000000..9966382 --- /dev/null +++ b/Sheet7/E14/jacob_template/visualize_par_results.m @@ -0,0 +1,52 @@ +%% Visualize results +% +% flatpak run org.octave.Octave +% or +% octave --no-window-system --no-gui -qf +% +% or +% +% matlab -nosplash -nodesktop -r 'try visualize_par_results(4); catch; end; quit' +% +function visualize_par_results(nprocs) +%% +if nargin<1 + nprocs = 4; +end +fprintf('# procs = %d\n',nprocs) + +pre = 'uv_'; +post = '.txt'; + +xc = []; nnodes = []; +ia = []; nelems = []; +v = []; +node_offset = 0; +elem_offset = 0; +for rank=0:nprocs-1 + fname = [pre,num2str(rank,'%2u'),post]; + [lxc,lia,lv] = ascii_read_meshvector(fname); +% whos lxc lia lv + nnodes = [nnodes size(lxc,1)]; + nelems = [nelems size(lia,1)]; + %[xc,ia,v] + xc = [xc; lxc]; + v = [v ; lv ]; + ia = [ia; lia+node_offset]; +% node_offset +% lia = lia + node_offset +% ia = [ia; lia]; + % index offsets for next subdomain + node_offset = node_offset + nnodes(end); + elem_offset = elem_offset + nelems(end); +end + +% fname = 'uv.txt'; +% [xc,ia,v] = ascii_read_meshvector(fname); + +h = trisurf(ia, xc(:,1), xc(:,2), v); +xlabel('x'),ylabel('y'),zlabel('z') + +shading interp + +waitfor(h) % wait for closing the figure diff --git a/Sheet7/E14/jacob_template/visualize_results.m b/Sheet7/E14/jacob_template/visualize_results.m new file mode 100644 index 0000000..e40da61 --- /dev/null +++ b/Sheet7/E14/jacob_template/visualize_results.m @@ -0,0 +1,20 @@ +%% Visualize results +% +% flatpak run org.octave.Octave +% or +% octave --no-window-system --no-gui -qf +% +% or +% matlab -nosplash < + +clear all +clc + +%% +fname = 'uv.txt'; + +[xc,ia,v] = ascii_read_meshvector(fname); + +h = trisurf(ia, xc(:,1), xc(:,2), v); + +waitfor(h) % wait for closing the figure \ No newline at end of file