Thanks for the detailed and helpful answer, Prof. Bangerth!
Yes, I thought about it in the past few days and arrived at the same
conclusion--that it does not make sense to represent the explicit advective
term as a single finite element field. I settled on a separate RHSAssembler
abstract base class with virtual methods to be overloaded for the volume
and face integrals on each cell. The main solvers simply request these
integrals from this class in the course of WorkStream::run.
The subtlety with respect to these advection body/face terms arises in the
following: as I'm sure you're well-aware, HDG methods involve an FESystem
(Q, U) on each element which is parametrized to a global system in the
skeleton space flux U_hat only (I used an FE_FaceQ to represent it). So all
of the right-hand side forcing is *not* communicated to the global linear
system directly, but indirectly through the U equation of the (Q, U) local
system (which is inverted to eliminate everything in terms of U_hat). So I
can't assemble the arbitrary RHS contributions directly into the linear
system like in the tutorials. I have to compute them and add them to the
RHS of the element-local system, like in step-51.
To that end, I'm having some trouble. I'd like to be able to compute
numerical fluxes of the known discontinuous field u^k (i.e., sum jumps and
averages of u^k) at the quadrature points on each face. However,
- I can't seem to be able to access u^k(+) and u^k(-) using a
FEInterfaceValues fiv, because it's vector-valued
- this same procedure works fine for scalars. On any face, I can use
something like fiv.get_face_values(1).get_function_values(scalar_field,
qpt_vector) to get scalar(-) at the face quadrature points (and use
get_face_values(0) for the scalar(+) value)
Am I missing something? There must be an easy way to access jumps and
averages of vector-valued DG data on the right-hand side. Step-47 almost
does this, but all the vector-valued jumps appear in the bilinear form
rather than the right-hand side.
If not, I could maybe extract the jumps and averages component by
component, using a storage scheme like step-35, but I feel like there must
be a better way.
Many thanks,
Corbin
On Wednesday, June 9, 2021 at 9:55:41 PM UTC-4 Wolfgang Bangerth wrote:
> On 6/8/21 3:33 AM, Corbin Foucart wrote:
> >
> > I've written generic, implicit HDG/DG solvers that solve a second-order
> > diffusion type equation $\frac{du}{dt} - \nabla^2 u = f$ (where f is a
> generic
> > DG finite element field). As I'm expanding the code to a set of
> projection
> > methods, I'd like this $f$ to represent an explicitly computed advection
> term
> > on the right-hand side $\nabla\cdot (u^k \otimes u^k)$ at time $k+1$.
> >
> > I'd like to assemble the weak form of the advection term, something like
> > $-(∇v, u^k\otimes u^k)_K +(v, F^*(u)\cdot n)_{\partial K}$ (where F^*(u)
> is
> > the convective numerical flux) on the right hand side for every element
> K and
> > store that as a DG field.
>
> This is conceptually not the right approach. A finite element field is a
> function of the form
> u_h(x) = \sum_j U_j \psi_j(x)
> where U_j is a vector of coefficients. But you don't have a field of this
> form: You have a vector
> F_j
> = sum_K -(∇psi_j, u^k\otimes u^k)_K +(\psi_j, F^*(u)\cdot n)_{\partial K}
> but this is a vector of (weighted) integrals of the previous solution, not
> a
> vector of expansion coefficients.
>
> In other words, your vector F_j is a vector alright, but it does *not*
> corresponding to a finite element field, whether continuous or
> discontinuous.
>
>
> > To avoid breaking the abstraction of the implicit
> > solver classes, is it possible to compute this term outside the main
> solver
> > classes and pass it in as a DG Field represented as a Vector<double>?
>
> What you plan to do is no different, really, than any of the other right
> hand
> side vectors we assemble. For example, in step-4, we assemble a right hand
> side vector ("system_rhs") in assemble_system(). The only difference is
> that
> in your case, this right hand side vector depends on the previous
> solution,
> but conceptually you're doing the same thing.
>
> Of course, it may be the case that you're not doing any other assembly
> operations in each time step, for example because you keep using the same
> matrix every time. But you can always write a function of the form
> assemble_rhs() that only assembles a rhs vector. As just one example, the
> new
> step-77 program has a function compute_residual() that does essentially
> this.
>
> Best
> Wolfgang
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/0b8c10c0-5991-4123-9154-73dabd5e4f4an%40googlegroups.com.
