Skip to content

Fixes for BDDC#5152

Open
pbrubeck wants to merge 22 commits into
releasefrom
pbrubeck/fix/bddc
Open

Fixes for BDDC#5152
pbrubeck wants to merge 22 commits into
releasefrom
pbrubeck/fix/bddc

Conversation

@pbrubeck

@pbrubeck pbrubeck commented Jun 6, 2026

Copy link
Copy Markdown
Contributor

Description

Internal fixes required for BDDCPC to work on primal elasticity problems.

Adds new options -pc_bddc_use_divergence_mat and -pc_bddc_use_discrete_gradient. Their defaults are the correct ones for the Riesz map in the de Rham complex, H(curl) needs the gradient, while H(div) needs the divergence. The H1 Riesz map does not need an extra operator. However, elasiticity problems are posed on H1, but BDDC would need to be aware of the divergence, so in this case the user must prescribe it.

@pbrubeck pbrubeck force-pushed the pbrubeck/fix/bddc branch from e11e7da to 2f8b06c Compare June 6, 2026 12:25
Comment thread firedrake/preconditioners/bddc.py
Comment thread firedrake/preconditioners/bddc.py Outdated
Comment thread tests/firedrake/regression/test_bddc.py
@pbrubeck pbrubeck force-pushed the pbrubeck/fix/bddc branch from 804646e to 1bc6e36 Compare June 8, 2026 10:24
Comment thread tests/firedrake/regression/test_bddc.py Outdated
Comment thread tests/firedrake/regression/test_bddc.py
Comment thread firedrake/preconditioners/bddc.py Outdated
Comment thread firedrake/preconditioners/bddc.py Outdated
Comment thread firedrake/preconditioners/bddc.py Outdated
default values for options should guarantee robustness wrt H1,Hdiv and Hcurl problems
Comment thread tests/firedrake/regression/test_bddc.py Outdated
Comment thread tests/firedrake/regression/test_bddc.py Outdated
@@ -34,9 +34,12 @@ class BDDCPC(PCBase):
the options:
- ``'bddc_cellwise'`` to set up a MatIS on cellwise subdomains if P.type == python,
- ``'bddc_matfree'`` to set up a matrix-free MatIS if A.type == python,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Can you do matrix-free MatIS? Then the divergence matrix can be passed that way, I only need the action of its local part (in MATIS sense) on a given vector

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Yeah, we can

@pbrubeck pbrubeck Jun 15, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I think we need to fix support for matfree A in BDDC

@pbrubeck

Copy link
Copy Markdown
Contributor Author

@connorjward BDDC is younger than Submesh. Can this go into release? It does not break API.

@pbrubeck pbrubeck requested a review from connorjward June 24, 2026 08:32
num_entities = len(active_entities)

kernel_code = f"""
void compute_entity_coordinates(PetscScalar *coords, PetscScalar *cg1_coords) {{

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Sorry about the parloop. There is no better way of computing entity coordinates on extruded meshes.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

It's not an issue. Moving this to pyop3 is easy.

Comment on lines +423 to +424
cg1_offsets, cg1_flatmap = flatten_entity_ids(active_entities, cg1_closure_ids)
v_offsets, v_flatmap = flatten_entity_ids(active_entities, V_entity_ids)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

These data structures wouldn't be required if we use an op2.PermutedMap

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Why not use one then?

@pbrubeck pbrubeck Jun 24, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

It only occured to me just now. Only v_flatmap is permutation, cg1_flatmap will have repeated indices because it defines the dofs supported on the closure of an entity (and edge will include its vertices), so this case is not covered by a permutation. We would still need to flatten_entity_ids at least for the CG1 element, so it is better to keep the code as is.

@pbrubeck pbrubeck Jun 24, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Ideally we should be running a parloop on the mesh entities. I know, wait for pyop3...

@connorjward connorjward left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Seems more or less fine. I think release is OK because this is a sufficiently new feature and you're only touching code in one file.

Comment thread tests/firedrake/regression/test_bddc.py
Comment on lines +46 to +47
# On MacOSX the distributed right-hand side is bugged!
if DEFAULT_DIRECT_SOLVER == "mumps":

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

This comment does not explain the code changes. There's nothing here that's obviously specific to macOS

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

@stefanozampini did you report this?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

there's no way this can be reported officially to MUMPS people. It is a known issue on only some osx configurations

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

The comment would lead me to expect a check like: if os.system == "darwin" or similar.

Comment on lines +93 to +94
# "ksp_view": None,
# "ksp_monitor_singular_value": None,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Suggested change
# "ksp_view": None,
# "ksp_monitor_singular_value": None,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Why? it is a quick way to check if something goes wrong instead of having to retype the options

num_entities = len(active_entities)

kernel_code = f"""
void compute_entity_coordinates(PetscScalar *coords, PetscScalar *cg1_coords) {{

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

It's not an issue. Moving this to pyop3 is easy.

Comment thread firedrake/preconditioners/bddc.py
Comment on lines +46 to +49
# On MacOSX the distributed right-hand side is bugged!
if DEFAULT_DIRECT_SOLVER == "mumps":
sp.update({"bddc_pc_bddc_coarse_mat_mumps_icntl_20": 0})

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I think we can simply remove this from here

Suggested change
# On MacOSX the distributed right-hand side is bugged!
if DEFAULT_DIRECT_SOLVER == "mumps":
sp.update({"bddc_pc_bddc_coarse_mat_mumps_icntl_20": 0})

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

But shouldn't this be the default?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

What is the problem of leaving the option in? A test should just work, and this is a known MUMPS bug

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

But shouldn't this be the default?

No, the default is distributed right-hand side, which means no need to gather the rhs on the master process

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants