Assembling explicit advection term independently

[Corbin Foucart]

Hello all, 

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. 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>?

I suspect this can be done, but I'm unsure how to populate a Vector<double> directly with terms like these, since there's no global assembly to be done. Perhaps it can be filled directly with MeshWorker::mesh_loop?

Any advice would be appreciated!

Corbin

[Wolfgang Bangerth]

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: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Corbin Foucart]

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

[Corbin Foucart]

Upon reading the documentation closer, I think I've found the way to access u^k (+/-) across each interface

• 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.

Indeed,

        std::vector<Vector<double>> velocity_qpt_vals(n_face_qpts, Vector<double>(dim));

        fiv.get_fe_face_values(1).get_function_values(velocity_field, velocity_qpt_vals); 

accesses the values; the issue was I hadn't explicitly instantiated the Vector<double>s with length dim, which was causing a segfault.  However, I'm curious -- is there a better way to compute the jumps and averages of DG data other than requesting them via fe_face_values and computing them "by hand"?

Thank you,

Corbin

[Wolfgang Bangerth]

On 6/10/21 2:07 AM, Corbin Foucart wrote:
>
> Indeed,
>         std::vector<Vector<double>> velocity_qpt_vals(n_face_qpts,
> Vector<double>(dim));
>         fiv.get_fe_face_values(1).get_function_values(velocity_field,
> velocity_qpt_vals);
>
> accesses the values; the issue was I hadn't explicitly instantiated the
> Vector<double>s with length dim, which was causing a segfault. However, I'm
> curious -- is there a better way to compute the jumps and averages of DG data
> other than requesting them via fe_face_values and computing them "by hand"?

Not right now. The FEInterfaceValues class should eventually acquire functions
get_function_jumps
get_function_averages
get_function_gradient_jumps
get_function_gradient_averages
etc, similar to the FEValues::get_function_values function and friends. But
nobody has written these functions yet. It would of course be fantastic if you
could give this a try -- it would solve your problem, and would be very useful
to others as well! We'd love walking you through the steps to make that happen!

Best
W.

[Wolfgang Bangerth]

On 6/10/21 7:10 AM, Wolfgang Bangerth wrote:
>
> Indeed,
>          std::vector<Vector<double>> velocity_qpt_vals(n_face_qpts,
> Vector<double>(dim));
>          fiv.get_fe_face_values(1).get_function_values(velocity_field,
> velocity_qpt_vals);
>
> accesses the values; the issue was I hadn't explicitly instantiated the
> Vector<double>s with length dim, which was causing a segfault.

I forgot that earlier. Could you come up with a small program that shows the
segfault? We should have caught the wrong vector size via an assertion, and
I'd like to add that assertion.

[Corbin Foucart]

I forgot that earlier. Could you come up with a small program that shows the
segfault? We should have caught the wrong vector size via an assertion, and
I'd like to add that assertion.

Certainly--please find attached.

To clarify, deal.II *did* throw some type of exception, I just wasn't able to interpret it and opted to use GDB instead; I'm a Python programmer and not a C++ expert by any stretch of the imagination :)

The program terminates with:

An error occurred in line <3196> of file </home/foucartc/software/dealii-9.2.0/source/fe/fe_values.cc> in function
    void dealii::internal::do_function_values(const typename VectorType::value_type*, const dealii::Table<2, double>&, const dealii::FiniteElement<dim, spacedim>&, const std::vector<unsigned int>&, dealii::ArrayView<Number>, bool, unsigned int) [with int dim = 2; int spacedim = 2; VectorType = dealii::Vector<double>; typename VectorType::value_type = double]
The violated condition was:
    static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(values.size()), decltype(n_quadrature_points)>::type)>::type>(values.size()) == static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(values.size()), decltype(n_quadrature_points)>::type)>::type>(n_quadrature_points)
Additional information:
    Dimension 0 not equal to 3.

Stacktrace:
-----------
#0  /home/foucartc/CorbinFoucart_repos/dealii_applications/lib/libdeal_II.g.so.9.2.0: void dealii::internal::do_function_values<2, 2, dealii::Vector<double> >(dealii::Vector<double>::value_type const*, dealii::Table<2, double> const&, dealii::FiniteElement<2, 2> const&, std::vector<unsigned int, std::allocator<unsigned int> > const&, dealii::ArrayView<dealii::Vector<double>, dealii::MemorySpace::Host>, bool, unsigned int)
#1  /home/foucartc/CorbinFoucart_repos/dealii_applications/lib/libdeal_II.g.so.9.2.0: void dealii::FEValuesBase<2, 2>::get_function_values<dealii::Vector<double> >(dealii::Vector<double> const&, std::vector<dealii::Vector<dealii::Vector<double>::value_type>, std::allocator<dealii::Vector<dealii::Vector<double>::value_type> > >&) const
#2  ./step: DGVelocityField<2>::access_interface_values()
#3  ./step: main
--------------------------------------------------------
Aborted (core dumped)

As for implementing FEInterfaceValues::get_function_jumps et al: I have a deadline by which I'd really like these set of codes to be working, so I can't work on this in the next few weeks. However, I will have to implement that functionality anyway, so I wouldn't mind working on this sometime in July.

Best,

Corbin

[Wolfgang Bangerth]

On 6/10/21 12:58 PM, Corbin Foucart wrote:
>
> To clarify, deal.II *did* throw some type of exception, I just wasn't able to
> interpret it and opted to use GDB instead; I'm a Python programmer and not a
> C++ expert by any stretch of the imagination :)

Ah, I see. Yes, the error message is difficult to read, but at least it did
catch the issue.

> As for implementing FEInterfaceValues::get_function_jumps et al: I have a
> deadline by which I'd really like these set of codes to be working, so I can't
> work on this in the next few weeks. However, I will have to implement that
> functionality anyway, so I wouldn't mind working on this sometime in July.

That would be fantastic! There are some good resources linked to from here, if
you're unsure about the process:
https://dealii.org/participate.html

[Corbin Foucart]

That would be fantastic! There are some good resources linked to from here, if
you're unsure about the process:
https://dealii.org/participate.html

I'm ready to start working on this. I've read the participation FAQ, forked the repo, etc. This might be a basic question, but how should I go about re-compiling the library code once I've made changes to the source? Is there a makefile that makes and runs all the unit tests? I presume the new code that I use to test the functions should be added in dealii/tests/feinterface?

Thank you for the guidance in advance!

Corbin

 

[Wolfgang Bangerth]

On 7/19/21 1:56 PM, Corbin Foucart wrote:
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fdealii.org%2Fparticipate.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C96b24c1d9f0a4fbe746e08d94aef5069%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637623213971258194%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=NtX0aKBoVNC7XEJqVCSfSzLUlY4NU1j%2BIqje6JwMZ54%3D&reserved=0>

>
>
> I'm ready to start working on this. I've read the participation FAQ, forked
> the repo, etc. This might be a basic question, but how should I go about
> re-compiling the library code once I've made changes to the source? Is there a
> makefile that makes and runs all the unit tests? I presume the new code that I
> use to test the functions should be added in dealii/tests/feinterface?

You build deal.II from a build directory and call
cmake <SOURCE_DIR> && make
See here:
https://dealii.org/readme.html

To run the test suite, you'd then do the following steps in your build directory:
make setup_tests
ctest
or, if you only want to run tests in debug mode, then
ctest -R debug
That last step will take many hours on a single-processor machine. If you have
enough cores, you can do
ctest -j80 -R debug
and it will finish in about 30 minutes, but you really do need 80 cores for that.