Reading Time: 11 minutes

If your solver runs for hours and still gives you numbers you can’t trust, error control isn’t a luxury — it’s the difference between a result you can publish and one you have to retract.

Here’s what you need to know right away:

  • Error estimation tells you how wrong your solution is. Error control is what you do with that information to make the answer right.
  • The adaptive mesh refinement loop follows a concrete “Solve → Estimate → Mark → Refine” cycle — not theory, a real algorithm you can implement in FEniCS or COMSOL.
  • Tolerance tuning has practical rules: never set atol near machine epsilon (~10⁻¹⁵), scale it to your physical noise floor, and tighten incrementally.
  • Goal-oriented error control via adjoint/DWR methods targets specific quantities of interest (drag, max temperature, flux) rather than global error — a fundamentally different approach from standard energy-norm estimation.
  • h, p, and hp adaptivity are not interchangeable. Smoothness indicators and engineering reality checks (1–5% energy norm targets) provide concrete decision rules most textbooks skip.

The distinction matters because most researchers I talk to spend months tuning meshes manually without ever running an automated loop, and even more struggle with tolerance settings that either waste compute or hide real errors. This article covers practical error control strategies — adaptive mesh refinement, tolerance tuning, convergence criteria, and goal-oriented methods — with concrete examples and a clear decision framework for when to use each approach.

If you haven’t read our guide on error estimation and verification methods, start there. Post 520 answers “Are we solving the equations correctly?” This article answers “How do we actively manage the error to meet a target accuracy?”

Estimation vs Control: The Fundamental Distinction

The easiest way to think about this is through a simple analogy. Error estimation is like a dashboard warning light. Error control is when you actually pull over and fix the problem.

Mathematically, the distinction is clean:

  • Error estimation computes or bounds the error (or a proxy for it) without knowing the exact solution. It’s diagnostic.
  • Error control uses those estimates to modify the discretization — refining the mesh, raising polynomial degree, or adjusting solver tolerances — to bring the error within acceptable bounds.

You can estimate the error without ever controlling it. But you can’t control it without estimating it first.

This distinction becomes critical when you’re under real pressure. A student working on a conference deadline doesn’t need a certified global error bound — they need a practical estimate that tells them whether their mesh is fine enough. A journal reviewer, on the other hand, expects error control documentation that shows the convergence behavior was actively managed.

I recommend treating estimation and control as separate steps with different purposes. Estimate first to understand your error landscape. Control second to navigate it. Mixing them in your workflow leads to blind adjustments — tweaking tolerances or remeshing without knowing what’s actually causing the error.

The Adaptive Loop: Solve → Estimate → Mark → Refine

This is the workhorse of practical error control. The “Solve → Estimate → Mark → Refine” loop is not theoretical — it’s a concrete algorithm implemented in COMSOL, FEniCS, deal.II, MFEM, and several other production solvers. Understanding it gives you the framework to reason about any adaptive solver you encounter.

┌─────────┐    ┌────────────┐    ┌──────────┐    ┌──────────┐
│         │    │            │    │          │    │          │
│  Solve   │───→│  Estimate  │───→│  Mark    │───→│  Refine  │
│  (solve │    │  (compute  │    │  (select │    │  (split/ │
│  the PDE)│    │   error)   │    │  elements│    │  increase│
│         │    │            │    │  to refine)│  p-degree) │
└─────────┘    └────────────┘    └──────────┘    └──────────┘
     │                                   │
     │                                   │
     └───────────────────────────────────┘
               (repeat until tolerance met)

The adaptive loop — visuals will be added during final publication.

Here’s how each step actually works:

Step 1: Solve

Run your solver on the current discretization. In FEniCS, this is just solve(formulation, solution) or whatever form your solver takes. The key point: you need the solution on the current mesh before you can estimate the error on it.

Step 2: Estimate

Compute a local error indicator on each mesh element. This is where the two main approaches diverge:

  • Standard (energy-norm) estimation computes a global error proxy. It answers: “How wrong is my solution overall?”
  • Goal-oriented (adjoint/DWR) estimation computes an error proxy for a specific quantity of interest. It answers: “How wrong is the drag coefficient I care about?”

The estimate typically takes the form of a local error indicator η_K for each element K. Many implementations use a reconstructed gradient, a residual-based formula, or a dual-weight approach.

Step 3: Mark

This is the step most tutorials gloss over, but it’s where the real engineering happens. You don’t refine every element with significant error — that would be wasteful. Instead, you select which elements to refine using a marking strategy.

The most widely used strategy is Dörfler marking, named after its inventor. Here’s how it works:

  1. Compute all local error indicators η_K for every element K.
  2. Sort elements by error indicator value, descending.
  3. Accumulate the total error η_total = Σ η_K.
  4. Select elements in sorted order until the accumulated error exceeds a threshold: Σ_{marked} η_K ≥ θ · η_total
  5. The parameter θ typically ranges from 0.5 to 0.8.

Here’s what I’d choose: Start with θ = 0.6 for a balanced approach. If your simulation is tight on resources, use θ = 0.5 to refine fewer elements per cycle. If you’re pursuing certified accuracy for a journal, push θ toward 0.8 to catch more error early.

Did you know? Many researchers skip Dörfler marking entirely and use the simpler “maximum marking” strategy (refine every element with η_K > max(η_K)/10). It’s easy to code, but it often refines too aggressively, leading to wasted compute and potentially worse performance on large meshes.

If you don’t set θ at all and use automatic marking (default behavior in some solvers), expect unpredictable refinement cycles. The parameter controls your refinement budget.

Step 4: Refine

Split refined elements (h-adaptivity) or increase polynomial degree (p-adaptivity). In practice, this means:

  • h-adaptivity: Subdivide selected elements into smaller children. COMSOL does this automatically when you enable its adaptive solver. In FEniCS, you call mesh.resize() with refined subdivision counts.
  • p-adaptivity: Raise the polynomial order on elements with high error (useful for smooth solutions where h-refinement is expensive).
  • hp-adaptivity: Do both — refine in regions with singularities, raise polynomial degree in smooth regions.

Then loop back to Solve. Repeat until your error tolerance is met or you hit your compute budget.

The loop is iterative, not one-shot. Each cycle moves you closer to the target accuracy. Most production solvers run 5–15 cycles before convergence.

Tolerance Tuning: Setting rtol and atol Without Guesswork

Solver tolerances are where many researchers waste compute time — or worse, hide real errors behind artificially loose settings.

In FEniCS (and most PETSc-based solvers), you configure two primary tolerances:

from fenics import *

# Typical setup
problem = Problem(F, u)
solver = problem.solver()

# The two critical settings
solver.set_absolute_tolerance(1e-12)   # atol
solver.set_relative_tolerance(1e-6)    # rtol

Here’s what I wish every researcher understood about tolerance tuning:

The scale rule

atol should be scaled relative to your physical noise floor, not to abstract numbers. If your simulation models physical processes with inherent variability at 10⁻⁸ magnitude, an atol of 10⁻¹⁵ is meaningless — you’re solving for precision beyond the physics you’re modeling.

The machine epsilon trap

Never set tolerances near machine epsilon (~10⁻¹⁵). PETSc and most solvers use floating-point arithmetic, not exact arithmetic. Pushing tolerances this tight doesn’t give you mathematical precision — it gives you noise, divergence, or wasted iterations. I’ve seen papers report atol=1e-15 and wonder why their solver runs 10× longer than it should.

The incremental tightening rule

Tolerances should be tightened incrementally by single orders of magnitude. Start with loose settings (rtol=1e-3, atol=1e-8) to verify your solver runs. Then tighten gradually:

  1. rtol=1e-3 → verify convergence
  2. rtol=1e-6 → check iteration count stability
  3. rtol=1e-8 → verify results are physically reasonable

If iteration counts spike or diverge at a particular tolerance, you’ve found a practical ceiling. Don’t push past it.

The iteration count signal

Always monitor iteration counts alongside tolerance settings. A solver that takes 100 iterations at rtol=1e-3 and 500 iterations at rtol=1e-6 tells you the convergence is getting harder — not necessarily better. If iteration counts explode without corresponding accuracy gains, the tolerance is too tight for the problem structure.

Practical rule of thumb

rtol: 1e-6 is a reasonable default for most FEM problems
atol:  1e-10 to 1e-12, scaled to your problem magnitude

For stiff or poorly-conditioned problems (which is most of them), consider using a preconditioner alongside tighter tolerances. The FEniCS tutorial recommends solver.preconditioner('ilu') or solver.preconditioner('ml') for better scaling.

What to avoid: Setting tolerances once and never revisiting them. Your tolerance settings should be documented alongside your mesh settings, just as you would document the discretization scheme. If reviewers ask about accuracy, you should be able to say “I used rtol=1e-6 and confirmed that tightening to 1e-8 changed the result by less than 0.1%.”

Goal-Oriented Error Control: Adjoint Methods and DWR

This is where error control diverges fundamentally from the standard approach.

Standard error estimation asks: “How wrong is my solution overall?” Goal-oriented error estimation asks: “How wrong is the number I actually care about?”

The difference is enormous. You might be simulating fluid flow through a complex geometry. A global energy-norm error estimate might be tiny — your solution is “close” overall. But the drag coefficient on your body might still be 10% off because the error is concentrated in regions that matter for drag but not for the global norm.

Here’s what I’d choose: If your research output depends on a specific quantity (drag force, maximum temperature, heat flux through a boundary, concentration at a monitoring point), goal-oriented error control is almost always the right approach.

How adjoint methods work

The adjoint method (also called the Discrete Adjoint or Discrete Weighted Remainder — DWR — approach) works in two steps:

  1. Solve the primal problem (your original PDE). This gives you the solution u.
  2. Solve the adjoint problem (a related PDE derived from your quantity of interest). This gives you a sensitivity field z that tells you how changes in u affect your target.

The error in your quantity of interest Q is then approximated as:

Error(Q) ≈ z^T · R(u)

where R(u) is the residual of the primal solve and z is the adjoint solution. This formula is elegant because it turns the error in a single scalar (Q) into a computable quantity without knowing the exact solution.

The practical workflow:

  1. Define your quantity of interest (QoI) as a functional ℓ(u).
  2. Solve the primal problem to get u.
  3. Solve the adjoint problem with source term defined by ℓ.
  4. Compute the dual-weight error indicator η_K = z^T · R_K.
  5. Use Dörfler marking (θ=0.5–0.8) to select elements for refinement.
  6. Refine and loop.

The adjoint approach is implemented in dolfin-adjoint (FEniCS), in CART3D (Nemec & Aftosmis, 2007), and in several commercial solvers. Hartmann’s work on multitarget error estimation shows how this scales to multiple QoI (lift and drag simultaneously, for instance).

When to use adjoint methods

Use goal-oriented estimation when:

  • Your research question depends on a specific measurable quantity (not the full field)
  • You’re publishing simulation results that reviewers will scrutinize
  • The global error is small but your QoI is not (or vice versa)
  • You need to justify accuracy claims to a review board or journal

Use standard energy-norm estimation when:

  • You need certified global accuracy bounds
  • Your problem is inherently multi-target (no single QoI dominates)
  • The adjoint problem is too expensive or impossible to formulate

A practical note: The adjoint method doubles your solve workload (primal + adjoint). For expensive large-scale problems, this can be prohibitive. In those cases, consider surrogate-based estimation or simplified adjoint approximations.

h, p, and hp Adaptivity: Which Strategy When?

Not all mesh adaptation is equal. The choice between h, p, and hp adaptivity affects convergence rates, computational cost, and even the type of problems you can solve effectively.

Strategy What it does Best suited for Convergence rate Implementation complexity
h-adaptivity Refine mesh (split elements) Problems with singularities, discontinuities, sharp gradients Linear (O(N^(-1/3)) in 2D for energy norm) Low — most solvers support it natively
p-adaptivity Increase polynomial order Smooth solutions, well-resolved flows Exponential (O(e^(-bN))) for smooth problems Medium — requires hierarchical bases
hp-adaptivity Do both Mixed smoothness, complex geometries Optimal (combined exponential + local resolution) High — requires smoothness detection

The practical decision rule I use: Check your solution’s smoothness first. If the solution is smooth (no shocks, discontinuities, or boundary layers), use p-adaptivity — exponential convergence means fewer elements for the same accuracy. If you have singularities or discontinuities (which is common in real engineering problems), use h-adaptivity — polynomial refinement alone won’t help. If your problem has mixed characteristics (smooth in most regions but singular at a few points), hp-adaptivity is optimal but comes at the cost of implementation complexity.

Smoothness indicators

The Šolín et al. (2004) comparison shows how smoothness indicators guide the choice. A practical indicator is the geometric decay of hierarchical coefficients — if they decay rapidly, the solution is smooth and p-refinement will be efficient. If they plateau or grow, you need h-refinement.

Engineering reality check: Most practical problems fall somewhere between “smooth” and “singular.” The 1–5% energy norm target is a common engineering guideline — refine until the estimated error is below 1–5% of the solution magnitude. This avoids the theoretical trap of chasing infinite precision when your physical model already has inherent uncertainty.

When to Use Estimation-First vs Control-First

Not all error control workflows are created equal. Your choice between estimation-first and control-first depends on your stage in the research process.

Estimation-first approach

Use when:

  • Exploring a new problem or geometry
  • Tuning discretization schemes
  • Running parameter sweeps or sensitivity studies
  • Preparing for conference deadlines (need fast, reliable estimates)

Why: Estimation is cheap. It doesn’t modify your mesh or solver. You get diagnostic information immediately. This is the right approach when you’re still figuring out what you’re looking at.

What to do: Run a coarse mesh, estimate the error, note the patterns. If the error distribution makes sense (large near boundaries, small in the bulk), your discretization is working. If the error is concentrated where you don’t expect it, investigate the discretization scheme.

Control-first approach

Use when:

  • Producing publication-quality results
  • Submitting to journals or review boards that require documentation
  • Running production simulations where accuracy is non-negotiable
  • Communicating results to stakeholders who need certified accuracy

Why: Control modifies your mesh and solver settings. It’s expensive (more solves, more compute). But it guarantees your error is bounded. This is the right approach when you’re done exploring and ready to publish.

What to do: Run the full adaptive loop. Document the refinement cycles, the tolerance settings, and the final error estimate. This documentation is what reviewers ask for when they want to know “How accurate is your result?”

I recommend framing this as a decision tree: Are you exploring (estimation-first) or producing (control-first)? Your answer determines your workflow. Don’t run an adaptive loop for exploration — it wastes compute. Don’t use raw estimates for publication — reviewers expect documentation.

What to Avoid

Here are the most common mistakes I see researchers make when implementing error control:

1. Setting tolerances once and never revisiting them. Your tolerance settings should be part of your documented methodology. If you change meshes, discretizations, or solvers, re-check tolerances.

2. Confusing solver convergence with discretization accuracy. A solver that converges to rtol=1e-6 is not necessarily accurate — it’s just consistent. Accuracy depends on mesh resolution, discretization order, and physical model fidelity. Always run a convergence study (Post 520 covers this) alongside tolerance tuning.

3. Using the same θ for every problem. Dörfler marking’s θ parameter should be tuned per problem type. A θ of 0.5 for exploration, 0.6 for production, 0.8 for certification — not one-size-fits-all.

4. Running the adaptive loop too few times. Most solvers need 5–15 cycles before the error estimate stabilizes. Three cycles is almost never enough.

5. Ignoring iteration counts when tuning tolerances. If tightening rtol from 1e-3 to 1e-6 doubles your iteration count, ask whether the extra 0.1% accuracy is worth the 100% compute cost.

6. Using adjoint methods without understanding the QoI. The adjoint approach is powerful, but if your quantity of interest isn’t well-defined (or if the adjoint problem doesn’t exist), it won’t help. Verify the adjoint formulation before investing the doubled solve cost.

Summary + Next Steps

Error control is the bridge between “I solved the equations” and “I solved the equations correctly.” Here’s how to think about it practically:

Estimation first. Understand your error landscape with diagnostic estimates. Run a coarse mesh, compute error indicators, and look for patterns.

Control second. Use the “Solve → Estimate → Mark → Refine” loop to actively reduce error. Set θ=0.6 for Dörfler marking as a starting point. Tighten tolerances incrementally.

Choose h, p, or hp based on smoothness. Smooth solutions → p-adaptivity. Singularities → h-adaptivity. Mixed → hp-adaptivity (if you can handle the complexity).

Target specific quantities with adjoint methods. When your QoI matters more than the global field, the adjoint/DWR approach gives you targeted accuracy without wasted refinement.

Document everything. Tolerance settings, θ values, convergence cycles, and final error estimates should be part of your published methodology — not afterthoughts.

If you want to go deeper, explore our companion article on error estimation and verification methods for the foundation, and read our guide on adaptive time stepping for temporal error control alongside spatial adaptation.

Related Guides

External references used in this article: