Skip to content

Free surface equation

free_surface_equation

FreeSurfaceTerm(test_space, trial_space, dx, ds, dS, **kwargs)

Bases: BaseTerm

Free Surface term: u \dot n

Source code in g-adopt/gadopt/equations.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
def __init__(
    self,
    test_space: firedrake.functionspaceimpl.WithGeometry,
    trial_space: firedrake.functionspaceimpl.WithGeometry,
    dx: firedrake.Measure,
    ds: firedrake.Measure,
    dS: firedrake.Measure,
    **kwargs,
):
    self.test_space = test_space
    self.trial_space = trial_space

    self.dx = dx
    self.ds = ds
    self.dS = dS

    self.mesh = test_space.mesh()
    self.dim = self.mesh.geometric_dimension()
    self.n = firedrake.FacetNormal(self.mesh)

    self.term_kwargs = kwargs

FreeSurfaceEquation(test_space, trial_space, quad_degree=None, prefactor=1, **kwargs)

Bases: BaseEquation

Free Surface Equation.

Source code in g-adopt/gadopt/equations.py 
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
84
85
86
87
88
def __init__(
    self,
    test_space: firedrake.functionspaceimpl.WithGeometry,
    trial_space: firedrake.functionspaceimpl.WithGeometry,
    quad_degree: Optional[int] = None,
    prefactor: Optional[Number | firedrake.Constant] = 1,
    **kwargs
):
    self.test_space = test_space
    self.trial_space = trial_space
    self.mesh = trial_space.mesh()

    p = trial_space.ufl_element().degree()
    if isinstance(p, int):  # isotropic element
        if quad_degree is None:
            quad_degree = 2*p + 1
    else:  # tensorproduct element
        p_h, p_v = p
        if quad_degree is None:
            quad_degree = 2*max(p_h, p_v) + 1

    if trial_space.extruded:
        # Create surface measures that treat the bottom and top boundaries similarly
        # to lateral boundaries. This way, integration using the ds and dS measures
        # occurs over both horizontal and vertical boundaries, and we can also use
        # "bottom" and "top" as surface identifiers, for example, ds("top").
        self.ds = CombinedSurfaceMeasure(self.mesh, quad_degree)
        self.dS = (
            firedrake.dS_v(domain=self.mesh, degree=quad_degree) +
            firedrake.dS_h(domain=self.mesh, degree=quad_degree)
        )
    else:
        self.ds = firedrake.ds(domain=self.mesh, degree=quad_degree)
        self.dS = firedrake.dS(domain=self.mesh, degree=quad_degree)

    self.dx = firedrake.dx(domain=self.mesh, degree=quad_degree)

    # General prefactor multiplying all terms in the equation
    # N.b. setting this to a firedrake constant (i.e. prefactor = Constant(1)) breaks
    # Drucker-Prager rheology test even though it is being multiplied by 1...
    self.prefactor = prefactor

    self.kwargs = kwargs

    # self._terms stores the actual instances of the BaseTerm-classes in self.terms
    self._terms = []
    for TermClass in self.terms:
        self._terms.append(TermClass(test_space, trial_space, self.dx, self.ds, self.dS, **kwargs))

mass_term(test, trial)

Return the UFL for the mass term \int test * trial * ds for the free surface time derivative term integrated over the free surface.

Source code in g-adopt/gadopt/free_surface_equation.py 
42
43
44
45
46
47
48
def mass_term(self, test, trial):
    r"""Return the UFL for the mass term \int test * trial * ds for the free surface time derivative term integrated over the free surface."""
    free_surface_id = self.kwargs['free_surface_id']
    k = self.kwargs['k']
    n = FacetNormal(self.mesh)
    # self.prefactor is a general prefactor from equations.py that multiplies all terms in an equation
    return self.prefactor * dot(n, k) * dot(test, trial) * self.ds(free_surface_id)