Specify all discrete FEM spaces when discretizing forms

[yguclu]

Current Situation

At the moment, Psydac requires a certain number of spaces when calling discretize on a symbolic form:

• For bilinear forms, it wants (U_h, V_h) where U_h is the trial space and V_h is the test space;
• For linear forms, it wants the test space V_h;
• For functionals, it wants a unique space W_h as it assumes that the expression to be integrated depends on some field wh in W_h.

These expressions could depend on additional "free fields", which are passed to the assemble function later on. Only at that point the discrete spaces of the free fields are read, which is in contrast to the discrete FEM spaces U_h, V_h, and W_h above. It is not clear whether proper checks on these spaces are made at the assemble function call.

Proposed Improvement

As a first step towards providing a consistent and safe interface, we propose to have discretize require the positional parameter discrete_spaces_map : dict[FunctionSpace, FemSpace]. This consists of a dictionary where the keys are the SymPDE spaces of every function which appears in the symbolic expression, and the values are the corresponding discrete spaces to be used in Psydac.

The dictionary discrete_spaces_map will be required to contain all the necessary spaces, but it could contain more. This allows one to create a dictionary once and for all, and pass it to every call to discretize for bilinear forms, linear forms, and functionals. This approach minimizes the sources of user error, and it makes the call to discretize a unique safe point where consistency checks are made in Psydac.

As added bonus, the dictionary discrete_spaces_map appears very useful if one wants to create discrete objects with different levels of grid refinement; see the example below.

Example

domain = Square('Omega')
S = ScalarFunctionSpace('S', domain)
V = VectorFunctionSpace('V', domain)

u, v = elements_of(V, names='u, v')
w = element_of(S, name='h')

a = BilinearForm((u, v),  integral(domain, dot(curl(u), curl(v)) * w))
l = LinearForm(v, integral(domain, dot(f, v * w)))  # f is some analytical expression

l2norm = Norm(u - f, domain, kind='l2')

domain_h1 = discretize(domain, ncells=nc)
domain_h2 = discretize(domain, ncells=nc*2)

V_h1 = discretize(V, domain_h1, degree=p)
V_h2 = discretize(V, domain_h2, degree=p)

S_h1 = discretize(S, domain_h1, degree=p+1)
S_h2 = discretize(S, domain_h2, degree=p+1)

# NEW: Dictionary of spaces which maps SymPDE analytical spaces to Psydac discrete spaces
spaces_h1 = {V: V_h1, S: S_h1}

# OLD vs. NEW interface
# a_h1 = discretize(a, domain_h1, (V_h1, V_h1), nquads=p+1)
a_h1 = discretize(a, domain_h1, spaces_h1, nquads=p+1) # Check: all spaces have the same domain

# l_h1 = discretize(a, domain_h1, V_h1, nquads=p+2)
l_h1 = discretize(l, domain_h1, spaces_h1, nquads=p+2)

# l2norm_h1 = discretize(l2norm, domain_h1, V_h1, nquads=p+1)  #-- but why?
l2norm_h1 = discretize(l2norm, domain_h1, spaces_h1, nquads=p+1)