Test Hodge decomposition

Two cases are tested

• Case 1 from ‘Décomposition de Hodge-Helmholtz discrète’ A. Lemoine et al., 21ème Congrès Français de Mécanique, 2013. The field is composed of vortex and sinks/sources, the 2 dimensional vector U=(u,v) is \displaystyle u = \sin(\pi x) \cos(\pi y) + \cos(\pi x) \sin(\pi y) \displaystyle v = \cos(\pi x) \sin(\pi y) - \sin(\pi x) \cos(\pi y)
• Case 2 : the 2D field comes from the scalar fonction -e^{-(x^2+y^2)} then \displaystyle u = (2x - 2 y) e^{-(x^2+y^2)} \displaystyle v = (2y+ 2 x ) e^{-(x^2+y^2)}

We note that the Hodge decompostion needs a function with a compact support, then the boundary conditions are very important. In the Case 2 as -e^{-(x^2+y^2)} goes to zero rapidly we have not problem if L is large enough.

#include "run.h"

#include "hodge.h"

int LEVEL = 6;
int FLOW;

vector u[];
vector gam[];
scalar func[];
scalar alpha[];

int main()
{
size (10);

origin (-L0/2,-L0/2);

init_grid (1 << LEVEL);

for (FLOW = 1; FLOW <= 1; FLOW++)
{
run();
}

}

event init (i=0)
{
switch (FLOW) {
case 1:
foreach()
{
u.x[] = sin(pi*x)*cos(pi*y) + cos(pi*x)*sin(pi*y);
u.y[] = cos(pi*x)*sin(pi*y) - sin(pi*x)*cos(pi*y);
}
break ;
case 2:
foreach()
{
u.x[] = 2.*x*exp(-(sq(x)+sq(y)));
u.y[] = 2.*y*exp(-(sq(x)+sq(y)));
u.x[] -= 2.*y*exp(-(sq(x)+sq(y)));
u.y[] += 2.*x*exp(-(sq(x)+sq(y)));
}
break ;
}

boundary  ((scalar *){u});

}

event print (t=end)
{

char name[80];
sprintf (name, "log-case%d", FLOW);
static FILE * fp = fopen (name, "w");

foreach()
{
fprintf (fp,"%f %f %f %f %f %f %f %f %f %f \n",x,y,u.x[],u.y[],gradalpha.x[],
}
fclose(fp);

if (FLOW==2)
{
foreach()
printf("%f %f %f %f \n",x,y,func[],-exp(-(sq(x)+sq(y))));
}
}
set multiplot layout 2,2 rowsfirst
a = 2
scale = 10
file = "log-case1"
# --- GRAPH a
plot [-a:a]  [-a:a] file  u 1:2:($3/scale):($4/scale) with vect t "original"
# --- GRAPH b
plot [-a:a]  [-a:a]  file u 1:2:($5/scale):($6/scale) with vect  t "rot = 0"
# --- GRAPH c
plot [-a:a]  [-a:a] file u 1:2:($7/scale):($8/scale) with vect  t "div = 0"
# --- GRAPH d
plot [-a:a]  [-a:a] file u 1:2:(($5+$7)/scale):(($6+$8)/scale) with vect  t "reconstruction"
#plot [-a:a]  [-a:a] file u 1:2:($9/scale):($10/scale) with vect
unset multiplot