1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495 |
- #include "ra/ra.hh"
- #include "ra/test.hh"
- #include "examples/cghs.hh"
- #include "ra/bench.hh"
- #include <numbers>
- using std::cout, std::endl, ra::TestRecorder;
- using ra::sqrm;
- constexpr auto PI = std::numbers::pi_v<double>;
- Benchmark::clock::duration tmul(0);
- struct StiffMat { double h; };
- template <class V, class W>
- void mult(StiffMat const & A, V const & v, W & w)
- {
- auto t0 = Benchmark::clock::now();
- auto i = ra::iota(v.len(0)-2, 1);
- auto j = ra::iota(v.len(1)-2, 1);
- w(i, j) = 4*v(i, j) -v(i-1, j) -v(i, j-1) -v(i+1, j) -v(i, j+1);
- tmul += (Benchmark::clock::now()-t0);
- }
- struct MassMat { double h; };
- template <class V, class W>
- void mult(MassMat const & M, V const & v, W & w)
- {
- auto t0 = Benchmark::clock::now();
- auto i = ra::iota(v.len(0)-2, 1);
- auto j = ra::iota(v.len(1)-2, 1);
- w(i, j) = sqrm(M.h) * (v(i, j)/2. + (v(i-1, j) + v(i, j-1) + v(i+1, j) + v(i, j+1) + v(i+1, j+1) + v(i-1, j-1))/12.);
- tmul += Benchmark::clock::now()-t0;
- }
- inline double f(double x, double y)
- {
- return 8.*PI*PI*sin(2.*PI*x)*sin(2.*PI*y);
- }
- inline double g(double x, double y)
- {
- return sin(2.*PI*x)*sin(2.*PI*y);
- }
- int main(int argc, char *argv[])
- {
- TestRecorder tr(std::cout);
- auto t0 = Benchmark::clock::now();
- int N = 50;
- double EPS = 1e-5;
- double h = 1./N;
- ra::Big<double, 2> v({N+1, N+1}, 0.), u({N+1, N+1}, 0.), w({N+1, N+1}, 0.), b({N+1, N+1}, 0.);
- auto i = ra::iota(N-1, 1);
- auto j = ra::iota(N-1, 1);
- auto ih = ra::iota(N-1, h, h);
- auto jh = ra::iota(N-1, h, h);
- v(i, j) = from(g, ih, jh);
- w(i, j) = from(f, ih, jh);
- StiffMat sm { h };
- MassMat mm { h };
- mult(mm, w, b);
- ra::Big<double, 3> work({3, N+1, N+1}, ra::none);
- int its = cghs(sm, b, u, work, EPS);
- double max = amax(abs(u-v));
- auto ttot = Benchmark::clock::now()-t0;
- cout << "total time " << Benchmark::toseconds(ttot)/1e-3 << " " << "ms" << endl;
- cout << "mul time " << Benchmark::toseconds(ttot)/1e-3 << " " << "ms" << endl;
- tr.info("max ", max).test_le(max, 0.00463);
- tr.info("its ", its).test_le(its, 67);
- return tr.summary();
- }
|