Basilisk source code (http://basilisk.fr/src/)

root / src / radial.h.page

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
/**
# Radial/cylindrical coordinates

This implements the radial coordinate mapping illustrated below.

![Radial coordinate mapping](figures/radial.svg)

It works in 1D, 2D and 3D. The 3D version corresponds to cylindrical
coordinates since the $z$-coordinate is unchanged. 

The only parameter is $d\theta$, the total angle of the sector. */

double dtheta = pi/3.;

/**
For convenience we add definitions for the radial and angular
coordinates $(r, \theta)$. */

map {
  double r = x, theta = y*dtheta/L0;
  NOT_UNUSED(r); NOT_UNUSED(theta);
}

event metric (i = 0) {

  /**
  We initialise the scale factors, taking care to first allocate the
  fields if they are still constant. */

  if (is_constant(cm)) {
    scalar * l = list_copy (all);
    cm = new scalar;
    free (all);
    all = list_concat ({cm}, l);
    free (l);
  }
  if (is_constant(fm.x)) {
    scalar * l = list_copy (all);
    fm = new face vector;
    free (all);
    all = list_concat ((scalar *){fm}, l);
    free (l);
  }

  /**
  The area (in 2D) of a mapped element is the area of an annulus
  defined by the two radii $r-\Delta/2$ and $r+\Delta/2$, divided by
  the total number of sectors $N=2\pi L0/(d\theta\Delta)$, this gives
  $$
  \frac{\pi [(r + \Delta / 2)^2 - (r - \Delta / 2)^2]}{2 \pi L 0 / (d \theta
  \Delta)} = \frac{2 \pi r \Delta}{2 \pi L 0 / (d \theta \Delta)} = \frac{rd
  \theta}{L 0} \Delta^2
  $$
  By definition, the (area) metric factor `cm` is the mapped area divided
  by the unmapped area $\Delta^2$. */

  scalar cmv = cm;
  foreach()
    cmv[] = r*dtheta/L0;

  /**
  It is important to set proper boundary conditions, in particular
  when refining the grid. */
  
  cm[left] = dirichlet (r*dtheta/L0);
  cm[right] = dirichlet (r*dtheta/L0);
  
  /**
  The (length) metric factor `fm` is the ratio of the mapped length of
  a face to its unmapped length $\Delta$. In the present case, it is
  unity for all dimensions except for the $x$ coordinates for which it
  is the ratio of the arclength by the unmapped length $\Delta$. We
  also set a small minimal value to avoid division by zero, in the
  case of a vanishing inner radius. */
  
  face vector fmv = fm;
  foreach_face()
    fmv.x[] = 1.;
  foreach_face(x)
    fmv.x[] = max(r*dtheta/L0, 1e-20);
  
  boundary ({cm, fm});
}