# GfsPoisson

GfsPoisson solves the Laplace equation ∇² P = 0, by default, for a scalar variable P; however, if a quantity Div is defined, typically in a GfsInit, this becomes the right-hand side Θ of the Poisson equation ∇² P = Θ.

Boundary conditions can be set on an embedded GfsSurface with GfsSurfaceBc, or on a GfsBoundary with a GfsBc; in either case the default is homogeneous Neumann.

### Examples

• Convergence of the Poisson solver
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine LEVEL

GModule hypre
ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }

Init {} {
Div = {
int k = 3, l = 3;
return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) print CYCLE, \$8;}' >> time
}
OutputProjectionStats { start = end } {
awk '{if (\$1 == "residual.infty:") print CYCLE, \$3, \$4;}' >> proj
}
OutputErrorNorm { start = end } {
awk '{print LEVEL, \$5, \$7, \$9}' >> error
} { v = P } {
s = (sin (M_PI*3.*x)*sin (M_PI*3.*y))
unbiased = 1
}
OutputSimulation { start = end } end-SOLVER.gfs
}
```

• Convergence with a refined circle
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine (x*x + y*y <= 0.25*0.25 ? LEVEL + 2 : LEVEL)

GModule hypre
ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }

Init {} {
Div = {
int k = 3, l = 3;
return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) print CYCLE, \$8;}' >> time
}
OutputProjectionStats { start = end } {
awk '{if (\$1 == "residual.infty:") print CYCLE, \$3, \$4;}' >> proj
}
OutputErrorNorm { start = end } {
awk '{print LEVEL, \$5, \$7, \$9}' >> error
} { v = P } {
s = (sin (M_PI*3.*x)*sin (M_PI*3.*y))
unbiased = 1
}
OutputSimulation { start = end } end-SOLVER.gfs
}
```

• Dirichlet boundary condition
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine LEVEL

Solid ({
double theta = atan2 (y, x);
double radius = 0.30 + 0.15*cos (6.*theta);
})
SurfaceBc P Dirichlet {
double theta = atan2 (y, x);
double r2 = x*x + y*y;
return r2*r2*cos (3.*theta);
}

GModule hypre
ApproxProjectionParams {
tolerance = 1e-30
erelax = 2
nitermin = CYCLE nitermax = CYCLE
}

Init {} {
Div = {
double theta = atan2 (y, x);
double r2 = x*x + y*y;
return 7.*r2*cos (3.*theta);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) print CYCLE, \$8;}' >> time
}
OutputProjectionStats { start = end } {
awk '{if (\$1 == "residual.infty:") print CYCLE, \$3, \$4;}' >> proj
}
OutputErrorNorm { start = end } {
awk '{print LEVEL, \$5, \$7, \$9}' >> error
} { v = P } {
s = {
double theta = atan2 (y, x);
double r2 = x*x + y*y;
return r2*r2*cos (3.*theta);
}
}
OutputSimulation { start = end } end-SOLVER.gfs
}
```

• Convergence of the Poisson solver with solid boundaries
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine LEVEL
Solid (ellipse (0, 0, 0.25, 0.25))

GModule hypre
ApproxProjectionParams { tolerance = 1e-30 erelax = 2 nitermin = CYCLE nitermax = CYCLE }

Init {} {
Div = {
int k = 3, l = 3;
return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) {print "CYCLE " \$8}}' >> time
}
OutputProjectionStats { start = end } {
awk '{
if (\$1 == "residual.infty:") print "CYCLE "\$3 " " \$4;
}' >> proj
}
OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
}
GfsBox {}
```

• Star-shaped solid boundary with refinement
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine LEVEL
RefineSolid (LEVEL + 2)
Solid ({
double dr = 0.1;
double theta = atan2 (y, x);
double radius = 0.79*(0.45 - dr + dr*cos (6.*theta));
})
GModule hypre
ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
Init {} {
Div = {
int k = 3, l = 3;
return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) {print "CYCLE " \$8}}' >> time
}
OutputProjectionStats { start = end } {
awk '{
if (\$1 == "residual.infty:") print "CYCLE "\$3 " " \$4;
}' >> proj
}
OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
}
GfsBox {}
```

• Star-shaped solid boundary
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine LEVEL
Solid ({
double dr = 0.1;
double theta = atan2 (y, x);
double radius = 0.79*(0.45 - dr + dr*cos (6.*theta));
})
GModule hypre
ApproxProjectionParams { erelax = 2 tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
Init {} {
Div = {
int k = 3, l = 3;
return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) {print "CYCLE " \$8}}' >> time
}
OutputProjectionStats { start = end } {
awk '{
if (\$1 == "residual.infty:") print "CYCLE "\$3 " " \$4;
}' >> proj
}
OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
}
GfsBox {}
```

• Thin wall at box boundary
• ```4 3 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
Refine LEVEL
GModule hypre
ApproxProjectionParams { erelax = 2 tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }
Init {} {
Div = {
int k = 3, l = 3;
x = (x - 0.5)/2.;
y = (y + 0.5)/2.;
return -M_PI*M_PI*(k*k + l*l)*sin (M_PI*k*x)*sin (M_PI*l*y);
}
}
OutputTime { istep = 1 } {
awk '{if (\$2 == 1) {print "CYCLE " \$8}}' >> time
}
OutputProjectionStats { start = end } {
awk '{
if (\$1 == "residual.infty:") print "CYCLE "\$3 " " \$4;
}' >> proj
}
OutputSimulation { start = end } sim-LEVEL-SOLVER { variables = P }
}
GfsBox {}
GfsBox {}
GfsBox {}
GfsBox {}
1 2 right
2 3 bottom
3 4 left
```

• Poisson solution in a dumbbell-shaped domain
• ```1 0 GfsPoisson GfsBox GfsGEdge {} {
Refine 3
ApproxProjectionParams { nitermax = 1000 minlevel = 1 tolerance = 1e-30 }
Solid dumbell.gts
Init {} { Div = y }
OutputProjectionStats { istep = 1 } stdout
}
```

• Poisson problem with a pure spherical harmonics solution
• ```6 12 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
PhysicalParams { L = 2.*M_PI/4. }
MetricCubed M LEVEL

# Use alternative implementation
# MetricCubed1 M
# MapFunction {
#     x = atan2 (X, Z)*180./M_PI
#     y = asin (CLAMP(Y,-1.,1.))*180./M_PI
# }

Refine LEVEL

GModule hypre
ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }

Global {
#include <gsl/gsl_sf.h>

double fact (double n) {
if (n <= 1)
return 1.;
else
return n*fact(n - 1.);
}

double legendre (int l, int m, double x) {
if (m < 0)
return pow(-1.,fabs(m))*fact(l-fabs(m))/fact(l+fabs(m))*
gsl_sf_legendre_Plm (l, fabs(m), x);
else
return gsl_sf_legendre_Plm (l, m, x);
}

double spherical_harmonic_re (int l, int m, double X, double Y) {
return sqrt((2*l+1)/(4*M_PI)*fact(l-m)/fact(l+m))*
legendre (l, m, sin(Y/180*M_PI))*cos(m*X/180*M_PI);
}
}

Init { } {
Div = -4*(4+1)*spherical_harmonic_re (4, 2, x, y)
Sol = spherical_harmonic_re (4, 2, x, y)
}

OutputTime { istep = 1 } {
awk '{if (\$2 == 1) print CYCLE, \$8;}' >> time
}
OutputProjectionStats { start = end } {
awk '{if (\$1 == "residual.infty:") print CYCLE, \$3, \$4;}' >> proj
}
OutputErrorNorm { start = end } {
awk '{print LEVEL, \$5, \$7, \$9}' >> error
} { v = P } {
s = Sol
v = E
unbiased = 1
}
OutputSimulation { start = end } end-SOLVER-LEVEL.gfs
}
```

• Spherical harmonics with longitude-latitude coordinates
• ```2 2 GfsPoisson GfsBox GfsGEdge {} {
Time { iend = 1 }
PhysicalParams { L = M_PI }
MetricLonLat M 1.
Refine LEVEL

GModule hypre
ApproxProjectionParams { tolerance = 1e-30 nitermin = CYCLE nitermax = CYCLE }

Global {
#include <gsl/gsl_sf.h>

double fact (double n) {
if (n <= 1)
return 1.;
else
return n*fact(n - 1.);
}

double legendre (int l, int m, double x) {
if (m < 0)
return pow(-1.,fabs(m))*fact(l-fabs(m))/fact(l+fabs(m))*
gsl_sf_legendre_Plm (l, fabs(m), x);
else
return gsl_sf_legendre_Plm (l, m, x);
}

double spherical_harmonic_re (int l, int m, double X, double Y) {
return sqrt((2*l+1)/(4*M_PI)*fact(l-m)/fact(l+m))*
legendre (l, m, sin(Y/180*M_PI))*cos(m*X/180*M_PI);
}
}

Init { } {
Div = -4*(4+1)*spherical_harmonic_re (4, 2, x, y)
Sol = spherical_harmonic_re (4, 2, x, y)
}

OutputTime { istep = 1 } {
awk '{if (\$2 == 1) print CYCLE, \$8;}' >> time
}
OutputProjectionStats { start = end } {
awk '{if (\$1 == "residual.infty:") print CYCLE, \$3, \$4;}' >> proj
}
OutputErrorNorm { start = end } {
awk '{print LEVEL, \$5, \$7, \$9}' >> error
} { v = P } {
s = Sol
v = E
unbiased = 1
}
OutputSimulation { start = end } end-SOLVER-LEVEL.gfs
}
```

• Poisson equation on a sphere with Gaussian forcing
• ```6 12 GfsPoisson GfsBox GfsGEdge {} {
PhysicalParams { L = 2.*M_PI/4. }
MetricCubed M LEVEL
Time { iend = 1 }
Refine LEVEL

ProjectionParams { tolerance = 1e-12 }
ApproxProjectionParams { tolerance = 1e-12 }

Init {} { Div = exp(-2.*5.*5.*(1.-cos((y + 90.)/180.*M_PI))) }

OutputLocation { start = end } prof.dat profile

OutputSimulation { start = end } end.gfs
OutputProjectionStats { istep = 1 } stderr

EventScript { start = end } {
gnuplot <<EOF
set term pos enhanced eps color solid 20 lw 3
set out 'profile.eps'

set key bottom right
set xl "Latitude"
set yl "{/Symbol F}"
set xr [-90:90]

plot './prof.dat' u 3:5 every 5 ps 2 t "Gerris",\
'prof.ref'u 1:2 w l t "(Boyd and Zhou, 2009)"
EOF
}

}
```

• Gaussian forcing using longitude-latitude coordinates
• ```2 2 GfsPoisson GfsBox GfsGEdge {} {
PhysicalParams { L = M_PI }
MetricLonLat M 1.
Time { iend = 1 }
Refine LEVEL

ProjectionParams { tolerance = 1e-12 }
ApproxProjectionParams { tolerance = 1e-12 }

Init {} { Div = exp(-2.*5.*5.*(1.-cos((y + 90.)/180.*M_PI))) }

OutputLocation { start = end } prof.dat profile

OutputSimulation { start = end } end.gfs
OutputProjectionStats { istep = 1 } stderr

EventScript { start = end } {
gnuplot <<EOF
set term pos enhanced eps color solid 20 lw 3
set out 'profile.eps'

set key bottom right
set xl "Latitude"
set yl "{/Symbol F}"
set xr [-90:90]

plot './prof.dat' u 3:(\\$5-0.0078011) every 5 ps 2 t "Gerris",\
'prof.ref'u 1:2 w l t "(Boyd and Zhou, 2009)"
EOF
}

}
```