/** # A proof of concept for the implicit-integral-equation solver We look for the field $S$, defined as the line integral of $s$ along the vector $\mathbf{n}$; $$S(\mathbf{x_0}) = \int_{\mathbf{x_0}} s\mathrm{d}\mathbf{n}$$ ![The field to integrate ($s$)](ti2/source.png) The code solves for $S$, by solving the implicit equation, $$\mathbf{n}\cdot\left(\nabla \left(\mathbf{n}\cdot\nabla S\right)\right) = -\mathbf{n}\cdot \nabla s,$$ in many directions: ![Looks good for all angles](ti2/S.mp4)(loop) ~~~gnuplot Convergence history as a function of the angle set polar set grid polar unset xtics unset ytics set xlabel 'MG cycles' set border 0 set style fill solid 0.5 set rrange [0.1 : 35] set size square set key right outside plot 'out' u 1:2 with filledcurve above r = 20 notitle , 'out' u 1:2 w l lw 2 t 'MG cycles' ,\ 'out' u 1:3 with filledcurve above r = 10 notitle , 'out' u 1:3 w l lw 2 t 'Relaxation sweeps' ~~~ The convergence for the near-grid-alliged integrations is poor. */ #include "integrator2.h" #include "utils.h" int main() { L0 = 6; X0 = Y0 = Z0 = -L0/2; init_grid (N); scalar s[], S[]; foreach() s[] = exp(-sq(x - 1) - sq(y - 1) - sq(z)) - exp(-sq(x + 1) - sq(y + 1) - sq(z)); boundary ({s}); output_ppm (s, file = "source.png", n = 300); for (double angle = 0; angle <= 2*pi + 1e-5; angle += pi/50) { coord n = (coord){cos(angle)- 0.01, sin(angle) + 0.01, 0}; mgstats lint = integrate_dn (S, s, n); printf ("%g %d %d\n", angle, lint.i, lint.nrelax);; output_ppm (S, file = "S.mp4", n = 300, min = -sqrt(pi), max = sqrt(pi)); } }