How to discretize equations

Now we demonstrate how to discretize a PDE object.

Pre-coded sample objects are stored in sample attribute of phyem. Invoke these samples by

import sys
ph_dir = '../../../../../'   # the path to dir that containing the phyem package.
sys.path.append(ph_dir)
import phyem as ph  # import the phyem package
ph.config._set_matplot_block(False)
samples = ph.samples

The partial differential equations of the canocical linear por-Hamiltonian are pre-coded. Call it through

oph = samples.pde_canonical_pH(n=3, p=3, periodic=True)[0]  # where o on `oph` means outer

Check oph, outer oriented port-Hamiltonian, with

oph.pr()

<Figure size 800x600 with 1 Axes>

We can take out the knowns and label them by

a3, b2 = oph.unknowns

We now test oph with test functions from the spaces where a3 and b2 come from, and label the test functions by \(v^3\) and \(u^2\),

wf = oph.test_with(oph.unknowns, sym_repr=[r'v^3', r'u^2'])

Now, we apply integration by parts to the term indexed by '1-1'.

wf = wf.derive.integration_by_parts('1-1')

We now apply a particular discretization to this weak formulation,

td = wf.td
td.set_time_sequence()  # initialize a time sequence
td.define_abstract_time_instants('k-1', 'k-1/2', 'k')
td.differentiate('0-0', 'k-1', 'k')
td.average('0-1', b2, ['k-1', 'k'])

td.differentiate('1-0', 'k-1', 'k')
td.average('1-1', a3, ['k-1', 'k'])
dt = td.time_sequence.make_time_interval('k-1', 'k')

wf = td()

wf.unknowns = [
    a3 @ td.time_sequence['k'],
    b2 @ td.time_sequence['k'],
]

wf = wf.derive.split(
    '0-0', 'f0',
    [a3 @ td.ts['k'], a3 @ td.ts['k-1']],
    ['+', '-'],
    factors=[1/dt, 1/dt],
)

wf = wf.derive.split(
    '0-2', 'f0',
    [ph.d(b2 @ td.ts['k-1']), ph.d(b2 @ td.ts['k'])],
    ['+', '+'],
    factors=[1/2, 1/2],
)

wf = wf.derive.split(
    '1-0', 'f0',
    [b2 @ td.ts['k'], b2 @ td.ts['k-1']],
    ['+', '-'],
    factors=[1/dt, 1/dt]
)

wf = wf.derive.split(
    '1-2', 'f0',
    [a3 @ td.ts['k-1'], a3 @ td.ts['k']],
    ['+', '+'],
    factors=[1/2, 1/2],
)

wf = wf.derive.rearrange(
    {
        0: '0, 3 = 2, 1',
        1: '3, 0 = 2, 1',
    }
)

We now can write the weak formulation with matrix proxies. Before doing that, we need to

ph.space.finite(3)
mp = wf.mp()

The matrix format of the weak formulation leads to a linear system,

ls = mp.ls()
ls.pr()

<Figure size 1200x600 with 1 Axes>

Note that, till this moment, everything is still abstract. To do the numerical simulation, we need to bring them to a particular implementation, for example, the msepy, standing for mimetic spectral elements python, by calling ph.fem.apply function, which will be shown in other notebooks.