Hp adaptivity through local smoothness & error estimate #4207

skvarjun · 2025-07-01T19:05:23Z

hp-Adaptivity algorithm : logic flow

Projecting local solution in each element to Legendre expansion.
Added a function legepoly() to compute tensor-product Legendre basis functions evaluated at a point p.
Use the classical recurrence relation to evaluate 1D Legendre polynomials L_n(x), extended to multiple dimensions.
Construct the basis up to a specified total degree order of the element.
For each element, extract the local dof and compute the finite element solution u(h) at quadrature points.
Assemble the projection matrix K and right-hand side F using K_ij = integral (psi_i psi_i) dx and F_i = u_h psi_i dx
Solve the system K.c=F to obtain the Legendre coefficients c per element.
Then index-to-Degree Mapping: Here, a lookup vector mapping each index in the coefficient vector c to its corresponding total polynomial degree.
This mapping gives the basis ordering used during legepoly() construction
Coefficient Grouping and Filtering: Group coefficients c_i by their associated total degree.
From |k_i|, omitting all coefficients where ∣k∣=0 to exclude the constant mode. Then accumulate all coefficients (U_k) of the same degree into vectors keyed by degree.
Then for each degree (k >= 1), compute the L2 norm.
Finally, log-Log regression for smoothness esimation.
Construct vectors log(∣U_k∣) and log(k) using filtered norms.
Apply least-squares linear regression as log(∣U_k∣) = −s.log(k) + C
Compute the negtaive of slope s (call it regularity), and define the smoothness index as for the given element.
Finally calculates p corsen (p_error) and h corsen error (h_error) (as @roystgnr composed in hp_coarsentest.C) and then if the element's local solution is having highest regularity than mean and p_error > h_error then do p refine else h refine.

ps: Negative value of regularity indicates the singularity of solution for the given element. This means higher order terms in Legendre expansion contaminated by singularity behavior and low order spectral modes unable be capture the solution behaviour.

jwpeterson

I made a few minor comments, but hp-adaptivity is really under @roystgnr's purview, so he will hopefully have a lot more to say here.

examples/adaptivity/adaptivity_ex3/adaptivity_ex3.C

jwpeterson · 2025-07-01T20:05:15Z

include/error_estimation/hp_coarsentest.h

@@ -97,7 +98,7 @@ class HPCoarsenTest : public HPSelector
   * refinement and potentially change the desired
   * refinement type.
   */
-  virtual void select_refinement (System & system) override;
+  virtual void select_refinement (System & system, ErrorVector & smoothness);


This is changing a public HPCoarsenTest API which we typically try to avoid doing without some period of deprecation. I would be curious to know what @roystgnr thinks about this particular API, though, I'm not sure it's used widely.

I sure hope it's not used widely - this class has never really worked well.

And for that reason, the constructor has always been marked libmesh_experimental(), and for both those reasons, anyone crazy enough to use this class should be on the ball enough to handle an API change without a transition period.

I'm not sure about this specific change, though.

If this class had actually worked, the original plan had been for both it and HPSingularity (and any future selection options) to become subclasses of HPSelector. And although a smoothness vector sure appears to be what we need here, I don't think it's what we need in general. Wouldn't it make just as much sense to create the smoothness vector inside this implementation of select_refinement, rather than creating it separately and passing it in?

include/error_estimation/smoothness_estimator.h

include/utils/error_vector.h

src/error_estimation/smoothness_estimator.C

moosebuild · 2025-07-02T23:18:28Z

Job Coverage, step Generate coverage on ae0e884 wanted to post the following:

Coverage

	5868ff	#4207 ae0e88
	Total	Total	+/-	New
Rate	64.28%	64.31%	+0.03%	81.58%
Hits	75108	75269	+161	155
Misses	41736	41763	+27	35

Diff coverage report

Full coverage report

Warnings

New new line coverage rate 81.58% is less than the suggested 90.0%

This comment will be updated on new commits.

roystgnr · 2025-07-10T22:28:55Z

Is this ready to review? I'm in the middle of one but I just had a comment submission invalidated by a force-push.

skvarjun · 2025-07-10T22:32:48Z

@roystgnr It is ready. I just added a warning message. Please continue reviewing (accedentally force pushed)

.DS_Store

examples/adaptivity/adaptivity_ex3/adaptivity_ex3.C

examples/adaptivity/adaptivity_ex3/adaptivity_ex3.in

src/utils/error_vector.C

roystgnr · 2025-07-10T22:30:27Z

src/error_estimation/smoothness_estimator.C

+#include "libmesh/libmesh_common.h"
+#include "libmesh/smoothness_estimator.h"
+#include "libmesh/dof_map.h"
+#include "libmesh/fe_base.h"
+#include "libmesh/dense_matrix.h"
+#include "libmesh/dense_vector.h"
+#include "libmesh/error_vector.h"
+#include "libmesh/libmesh_logging.h"
+#include "libmesh/elem.h"
+#include "libmesh/quadrature_grid.h"
+#include "libmesh/system.h"
+#include "libmesh/mesh_base.h"
+#include "libmesh/numeric_vector.h"
+#include "libmesh/tensor_value.h"
+#include "libmesh/threads.h"
+#include "libmesh/tensor_tools.h"
+#include "libmesh/enum_error_estimator_type.h"
+#include "libmesh/enum_norm_type.h"
+#include "libmesh/int_range.h"
+#include "libmesh/enum_to_string.h"


This is just the full list from patch_recovery_error_estimator.C, right? We need to trim it down to just the necessary parts.

I forgot that. I just used all headers in patch_recovery_error_estimator.C :(. Cleaned

roystgnr · 2025-07-10T22:33:47Z

src/error_estimation/smoothness_estimator.C

+    const Order element_order = fe_type.order + elem->p_level();
+
+    if (element_order == 1)
+      libmesh_warning("Warning: Polynomial order in element is 1, smoothness cannot be stimated!\n"


Suggested change

libmesh_warning("Warning: Polynomial order in element is 1, smoothness cannot be stimated!\n"

libmesh_warning("Warning: Polynomial order in element is 1, smoothness cannot be estimated!\n"

Although, I'm not sure we want a warning here at all. Surely we expect to hit element_order=1 cases legitimately, and we don't need to warn the user because they're not doing anything wrong?

But we should document (not here, but in the header) what does end up in the output ErrorVector in those elements' entries in those cases.

I agree with this.
For element order 1, the linear fit would be absurd for log|order| vs log |coefficients| . This results in Nan on smoothness vector, and pure h refinement triggered.

Ooh - can we return std::numeric_limits<Real>::max() instead of triggering a NaN? We have a lot of users who turn on floating point exceptions.

roystgnr · 2025-07-10T22:50:26Z

include/enums/enum_error_estimator_type.h

@@ -40,7 +40,8 @@ enum ErrorEstimatorType : int {
                         LAPLACIAN               =  5,
                         PATCH_RECOVERY          =  6,
                         WEIGHTED_PATCH_RECOVERY =  7,
-                         UNIFORM_REFINEMENT      =  8};
+                         UNIFORM_REFINEMENT      =  8,
+                         SMOOTHNESS              =  9};


Calling the SmoothnessEstimator an ErrorEstimator feels like a category error to me. We're returning a smoothness, which should be non-zero for the exact solution, rather than error, which shouldn't be, right?

Yes, smoothness shall be non-zero for the exact solution. So do you think the ENUM for smoothness is better to remove?

Yeah - even if it's super useful for estimating error, it's still not a type of error.

src/error_estimation/smoothness_estimator.C

src/.DS_Store

roystgnr · 2025-07-10T23:05:44Z

Oh, BTW, that was just a question about readiness; I didn't mean to criticize the force push. Force push after a negative review bugs me (because then it's much harder to see what's changed between reviews), but force push while I'm in the middle of a review just means I'm not reviewing fast enough. ;-) Plus, I told you I'd be on vacation for 2 weeks, so you probably weren't expecting me back until Monday, not today.

roystgnr

A couple new bits along with the remaining previous discussion.

src/utils/error_vector.C

roystgnr · 2025-07-14T20:45:24Z

src/error_estimation/smoothness_estimator.C

+#include "libmesh/error_vector.h"
+#include "libmesh/libmesh_logging.h"
+#include "libmesh/elem.h"
+#include "libmesh/quadrature_grid.h"


I think you missed this one; PREE uses a QGrid but you don't.

Missed that

I need quadrature_grid.h as I'm using at at a qp_loop later #line250.

The qrule here isn't a QGrid. Just include "libmesh/quadrature.h".

roystgnr · 2025-07-14T20:48:15Z

src/error_estimation/smoothness_estimator.C

+    const Order element_order = fe_type.order + elem->p_level();
+
+    if (element_order == 1)
+      libmesh_warning("Warning: Polynomial order in element is 1, smoothness cannot be stimated!\n"


Ooh - can we return std::numeric_limits<Real>::max() instead of triggering a NaN? We have a lot of users who turn on floating point exceptions.

roystgnr · 2025-07-14T20:48:34Z

include/enums/enum_error_estimator_type.h

@@ -40,7 +40,8 @@ enum ErrorEstimatorType : int {
                         LAPLACIAN               =  5,
                         PATCH_RECOVERY          =  6,
                         WEIGHTED_PATCH_RECOVERY =  7,
-                         UNIFORM_REFINEMENT      =  8};
+                         UNIFORM_REFINEMENT      =  8,
+                         SMOOTHNESS              =  9};


Yeah - even if it's super useful for estimating error, it's still not a type of error.

roystgnr · 2025-07-16T03:15:33Z

How and why did you pull in those three completely-unrelated commits?

roystgnr

Can you do a force push just to remove the three unrelated commits? If we have to bisect anything later I don't want to be as baffled looking at git history then as I am right now.

That plus the header choice may be all we have left to fix here.

roystgnr · 2025-07-16T03:14:08Z

src/error_estimation/smoothness_estimator.C

+#include "libmesh/error_vector.h"
+#include "libmesh/libmesh_logging.h"
+#include "libmesh/elem.h"
+#include "libmesh/quadrature_grid.h"


The qrule here isn't a QGrid. Just include "libmesh/quadrature.h".

skvarjun · 2025-07-17T20:07:39Z

which commits do you want removed/deleted?
or did you mean squashing them?

…output file

hp-Adaptivity: Added documentation for few lines hp-Adaptivity: Few minor corrections - hp adaptivity

… from ErrorEstimator

…nement, thus override is safe hp-Adaptivity: fixed headers

jwpeterson reviewed Jul 1, 2025

View reviewed changes

skvarjun force-pushed the hp_adaptivity branch from e1c2dc0 to be66288 Compare July 2, 2025 15:37

skvarjun changed the title ~~Hp adaptivity~~ Hp adaptivity through smoothness estimate Jul 2, 2025

skvarjun force-pushed the hp_adaptivity branch from be66288 to c473bc1 Compare July 2, 2025 18:42

skvarjun changed the title ~~Hp adaptivity through smoothness estimate~~ Hp adaptivity through smoothness & error estimate Jul 2, 2025

skvarjun force-pushed the hp_adaptivity branch from c473bc1 to e4f8549 Compare July 2, 2025 20:26

skvarjun changed the title ~~Hp adaptivity through smoothness & error estimate~~ Hp adaptivity through local smoothness & error estimate Jul 2, 2025

skvarjun force-pushed the hp_adaptivity branch 2 times, most recently from 5d90b22 to d7f0482 Compare July 10, 2025 22:25

roystgnr requested changes Jul 10, 2025

View reviewed changes

roystgnr requested changes Jul 14, 2025

View reviewed changes

roystgnr requested changes Jul 16, 2025

View reviewed changes

skvarjun force-pushed the hp_adaptivity branch 2 times, most recently from 24308ac to 393a126 Compare July 17, 2025 17:38

skvarjun added 5 commits July 17, 2025 15:47

Fixed an issue with plot_error() which failed to update the error in …

cd2e7e6

…output file

hp-Adaptivity: Added smoothness marker and updated hp algorithm

1f1d2f6

hp-Adaptivity: Added documentation for few lines hp-Adaptivity: Few minor corrections - hp adaptivity

hp-Adaptivity: NaN triggering for 1st order poly. fixed

f7ef6ff

hp-adaptivity: SmoothnessEstimator class is no longer a derived class…

703f11b

… from ErrorEstimator

hp-adaptivity: smoothness vector inside implementation of select_refi…

ae0e884

…nement, thus override is safe hp-Adaptivity: fixed headers

skvarjun force-pushed the hp_adaptivity branch from 393a126 to ae0e884 Compare July 17, 2025 21:54

	libmesh_warning("Warning: Polynomial order in element is 1, smoothness cannot be stimated!\n"
	libmesh_warning("Warning: Polynomial order in element is 1, smoothness cannot be estimated!\n"

Hp adaptivity through local smoothness & error estimate #4207

Are you sure you want to change the base?

Hp adaptivity through local smoothness & error estimate #4207

Uh oh!

Conversation

skvarjun commented Jul 1, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jwpeterson left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

moosebuild commented Jul 2, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Coverage

Warnings

Uh oh!

roystgnr commented Jul 10, 2025

Uh oh!

skvarjun commented Jul 10, 2025

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

roystgnr commented Jul 10, 2025

Uh oh!

roystgnr left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

roystgnr commented Jul 16, 2025

Uh oh!

roystgnr left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

skvarjun commented Jul 17, 2025

Uh oh!

skvarjun commented Jul 1, 2025 •

edited

Loading

moosebuild commented Jul 2, 2025 •

edited

Loading