Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add support for interior face integrals #1223

Merged
merged 73 commits into from
Dec 5, 2024
Merged

Conversation

samuelpmishLLNL
Copy link
Contributor

@samuelpmishLLNL samuelpmishLLNL commented Aug 26, 2024

This PR aims to implement a feature requested by @tupek2 to be able to integrate over interior faces. The main challenge is in handling DG/L2 spaces for this type of integral, as they have separate values on each side.


The the new kind of integral behaves similarly to the existing boundary integral, except that for L2 spaces it represents:

Screenshot 2024-11-15 at 12 01 29 PM

where the domain of integration $\Omega$ represents the specified collection of interior faces, $\phi_1, \phi_2$ are the basis functions from side 1 and side 2, and $u_1, u_2$ are the values of the DG field evaluated from side 1 and side 2. For example, consider a mesh made of two adjacent quadrilateral elements, with a L2<1, 1> field defined on it. The red disks represent the nodes of the DG field (with the positions of the nodes shifted inward slightly to exaggerate which side of the interior face they belong to), and the blue edge represents the domain $\Omega$.

dg_example_diagram

In this example, the basis functions that appear in the integrand are:

Screenshot 2024-11-15 at 3 57 10 PM

where $N_0, N_1$ are the usual linear interpolation shape functions in 1D.


Here is an example of a residual calculation involving a vector-valued DG space and interior faces:

  ...

  constexpr int p = 2; // polynomial order
  using test_space  = L2<p, dim>;
  using trial_space = L2<p, dim>;

  constexpr int DERIVATIVE = 1;

  Domain interior_faces = InteriorFaces(*mesh);

  // Construct the new functional object using the specified test and trial spaces
  Functional<test_space(trial_space)> residual(&fespace, {&fespace});

  // register a new integral calculation with this residual object
  residual.AddInteriorFaceIntegral(
      Dimension<dim-1>{}, DependsOn<0>{},
      [=](double /*t*/, auto X, auto velocity) {

        // compute the (unit) surface normal
        auto dX_dxi = get<DERIVATIVE>(X);
        auto n = normalize(cross(dX_dxi));

        // extract the velocity values from each side of the interface
        // note: the orientation convention is such that the normal 
        //       computed as above will point from from side 1->2
        auto [u_1, u_2] = velocity; 

        // carry out whatever calculations are necessary to compute f_1, f_2
        auto a = dot(u_2 - u_1, n);
        auto f_1 = u_1 * a;
        auto f_2 = u_2 * a;

        // note: the values returned in the tuple are integrated 
        // against the basis functions from their respective sides
        return serac::tuple{f_1, f_2};

      }, interior_faces);

  ...

  // evaluate the residual with nodal values specified by U
  mfem::Vector U = ...; 
  mfem::Vector r = residual(t, U);

Here's a table showing which argument types are passed to qfunction and which weight functions are applied to the output of qfunctions, depending on the type of integral and trial/test space (respectively):

Screenshot 2024-11-19 at 8 34 25 AM

@samuelpmishLLNL samuelpmishLLNL added the WIP Work in progress label Aug 26, 2024
@samuelpmishLLNL
Copy link
Contributor Author

/style

cmake_dependent_option(SERAC_ENABLE_BENCHMARKS "Enable benchmark executables" ON "ENABLE_BENCHMARKS" OFF)
if (ENABLE_BENCHMARKS)
set(SERAC_ENABLE_BENCHMARKS ON)
endif()
Copy link
Contributor

@chapman39 chapman39 Dec 3, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While cmake_dependent_option is admittedly weird I think it is okay as is. SERAC_ENABLE_BENCHMARKS generally shouldn't be set explicitly unless you want to enable benchmarks yet disable Serac's. Simply set ENABLE_BENCHMARKS.

SERAC_ENABLE_BENCHMARKS can only be set if ENABLE_BENCHMARKS is ON. If ENABLE_BENCHMARKS is ON, it will automatically set SERAC_ENABLE_BENCHMARKS to ON. Otherwise SERAC_ENABLE_BENCHMARKS is forced to OFF.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please re-add cmake_dependent_option(SERAC_ENABLE_BENCHMARKS "Enable benchmark executables" ON "ENABLE_BENCHMARKS" OFF)

Copy link
Collaborator

@tupek2 tupek2 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm personally OK with disabling the .lua for now. We may be able to re-prioritize in the spring to revive, but we cannot keep this PR up forever, the scope is already massive.

@samuelpmishLLNL
Copy link
Contributor Author

/style


// the dofs mfem returns for Hcurl include information about
// dof orientation, but not for triangle faces on 3D elements.
// So, we need to manually
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this is unsupported for now, should we abort or warn?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good question, I'll defer that decision for now (it appears that this is new code in the PR, but it's just a copy of existing code)


ordering = fes->GetOrdering();
// note: this assumes that all the elements are the same polynomial order
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you have an idea how we can enforce this? Check on the Domain creation? Doesn't need to be in this PR, but maybe open as issue.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure. We could just iterate through the elements and check that they agree, but that doesn't catch the (unlikely) case that the ParMesh was partitioned into chunks that each have constant (but different) p.

src/serac/numerics/functional/element_restriction.hpp Outdated Show resolved Hide resolved
src/serac/numerics/functional/functional.hpp Outdated Show resolved Hide resolved
@samuelpmishLLNL
Copy link
Contributor Author

/style

@samuelpmishLLNL samuelpmishLLNL merged commit 96e272a into develop Dec 5, 2024
2 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request ready for review Ready for active inspection by reviewers refactoring Change the structure of the underlying data structures user-request Issue requested by user
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants