Skip to content

[WIP] Debugging check that variables of discontinuous family live on just one elem->dim() #4345

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.

Sign up for GitHub

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

Closed
lindsayad wants to merge 6 commits into from
Closed

[WIP] Debugging check that variables of discontinuous family live on just one elem->dim() #4345

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
1e949f5
Debugging check that variables of discontinuous family live on just o…
lindsayad Dec 5, 2025
cdcae75
Early return in sanity check if mesh isn't prepared
lindsayad Dec 6, 2025
ddbbb0c
Associate var in NodeElem test with 1D elements
lindsayad Dec 6, 2025
51dee79
Update sanity check message
lindsayad Dec 6, 2025
315f6eb
Move code to dof_map and compute on-the-fly
lindsayad Dec 6, 2025
138e739
Add variable name(s) in discontinuity error message
lindsayad Dec 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 49 additions & 0 deletions src/base/dof_map.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3139,6 +3139,47 @@ void DofMap::reinit_static_condensation()
_sc->reinit();
}

#ifdef DEBUG
namespace
{
unsigned char determine_elem_dims(const MeshBase & mesh,
const std::set<subdomain_id_type> & subdomains)
{
unsigned char elem_dims = 0;
for (const auto * const elem : mesh.active_subdomain_set_element_ptr_range(subdomains))
elem_dims |= static_cast<unsigned char>(1u << cast_int<unsigned char>(elem->dim()));
Copy link
Copy Markdown
Contributor

@GiudGiud GiudGiud Dec 9, 2025

Choose a reason for hiding this comment

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

is this chatgpt or does this make sense to you?
like the << and |= on chars

Copy link
Copy Markdown
Contributor

@GiudGiud GiudGiud Dec 9, 2025

Choose a reason for hiding this comment

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

well copilot explained it to me! Pretty creative

Copy link
Copy Markdown
Member

@roystgnr roystgnr Dec 9, 2025

Choose a reason for hiding this comment

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

I believe this is the sort of thing that's more clearly written with std::bitset in C++ code? But this was the standard idiom for bit sets in C code; I'm old enough that until I saw your comment it didn't even occur to me that the code here might be confusing.

Copy link
Copy Markdown
Member Author

@lindsayad lindsayad Dec 9, 2025

Choose a reason for hiding this comment

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

I really like it for being extremely lightweight compared to a std::set

mesh.comm().bitwise_or(elem_dims);
return elem_dims;
}

void discontinuity_sanity_check(const System & sys,
const FEType & fe_type,
const std::set<subdomain_id_type> * const active_subdomains)
{
const auto continuity = FEInterface::get_continuity(fe_type);
if (continuity != DISCONTINUOUS && continuity != SIDE_DISCONTINUOUS)
return;

const auto & mesh = sys.get_mesh();
// This check won't be meaninful if the mesh isn't prepared
if (!mesh.is_prepared())
return;

const auto & var_subdomains = (active_subdomains && !active_subdomains->empty())
? *active_subdomains
: mesh.get_mesh_subdomains();
const auto var_elem_dims = determine_elem_dims(mesh, var_subdomains);
// Power of two trick. Note that unsigned char automatically promoted to int for arithmetic and bitwise operations
const bool more_than_one_bit_set = (var_elem_dims & (var_elem_dims - 1)) != 0;
libmesh_assert_msg(
!more_than_one_bit_set,
"Discontinuous finite element families are typically associated with discontinuous Galerkin "
"methods. If degrees of freedom are associated with different values of element dimension, "
"then generally this will result in a singular system after application of the DG method.");
Copy link
Copy Markdown
Contributor

@GiudGiud GiudGiud Dec 9, 2025

Choose a reason for hiding this comment

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

is that true for a disconnected 1D domain and a disconnected 2D domain?

For example, bunch of TH channels solved with DG connected to a volume solve solved with DG, sharing p, v variables

Copy link
Copy Markdown
Member

@roystgnr roystgnr Dec 9, 2025

Choose a reason for hiding this comment

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

Hahaha - I was wondering if anyone else was going to consider that case. I thought pointing it out would be too nit-picky, since e.g. it would be natural to use different variables for e.g. v_parallel in a channel vs v_x/v_y/v_z in a volume. But now you've got me thinking about pressure, and it's actually not crazy to have a discontinuous pressure variable that's got the same meaning regardless of element dimension.

Copy link
Copy Markdown
Member Author

@lindsayad lindsayad Dec 9, 2025

Choose a reason for hiding this comment

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

I think we can add means to relax the check if they're needed. But this code in here entirely removes the need for isLowerD code in the framework which was only used for sanity checking, and the MOOSE sanity checking itself was overzealous; like in ProjectionAux there is not a conceptual need for excluding lower-d blocks. The only conceptual need is that all the blocks be of the same dimension which is what this check is doing

Copy link
Copy Markdown
Member Author

@lindsayad lindsayad Dec 10, 2025

Choose a reason for hiding this comment

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

Yea, thinking about it more, I think I'm just going to close this and prefer to error in the places that truly can't support a mix of dimensions. No need to an overzealous error check here

}
}
#endif

unsigned int DofMap::add_variable(System & sys,
std::string_view var,
const FEType & type,
Expand All @@ -3153,6 +3194,10 @@ unsigned int DofMap::add_variable(System & sys,
if (active_subdomains)
libmesh_assert(this->comm().verify(active_subdomains->size()));

#ifdef DEBUG
discontinuity_sanity_check(sys, type, active_subdomains);
#endif

// Make sure the variable isn't there already
// or if it is, that it's the type we want
for (auto v : make_range(this->n_vars()))
Expand Down Expand Up @@ -3262,6 +3307,10 @@ unsigned int DofMap::add_variables(System & sys,
if (active_subdomains)
libmesh_assert(this->comm().verify(active_subdomains->size()));

#ifdef DEBUG
discontinuity_sanity_check(sys, type, active_subdomains);
#endif

// Make sure the variable isn't there already
// or if it is, that it's the type we want
for (auto ovar : vars)
Expand Down
4 changes: 3 additions & 1 deletion tests/systems/equation_systems_test.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,11 +182,13 @@ public:
MeshTools::Generation::build_line (mesh, 10, 0., 1., EDGE2);
auto node_elem = mesh.add_elem(Elem::build(NODEELEM));
node_elem->set_node(0, mesh.node_ptr(0));
node_elem->subdomain_id() = 1;
mesh.prepare_for_use();

EquationSystems es(mesh);
System &sys = es.add_system<System> ("SimpleSystem");
sys.add_variable("u", CONSTANT, MONOMIAL);
const std::set<subdomain_id_type> u_subs = {0};
sys.add_variable("u", CONSTANT, MONOMIAL, &u_subs);
es.init();
es.reinit();
}
Expand Down