Skip to content

Port SolCx analytic Stokes solution (uw.function.analytic.SolCx) - #223

Merged
lmoresi merged 2 commits into
developmentfrom
bugfix/analytic-solutions
Jun 8, 2026
Merged

Port SolCx analytic Stokes solution (uw.function.analytic.SolCx)#223
lmoresi merged 2 commits into
developmentfrom
bugfix/analytic-solutions

Conversation

@lmoresi

Copy link
Copy Markdown
Member

Summary

Ports the classic Velic SolCx exact Stokes solution into uw.function.analytic, giving UW3 an exact-solution reference for the canonical free-slip, discontinuous-viscosity benchmark (viscosity step eta_A|eta_B at x_c, density forcing (0, cos(pi x) sin(n pi z)), free-slip walls).

Mirrors the existing AnalyticSolNL port and establishes the reusable pattern for the rest of the Velic suite.

What's included

  • solCx.c / solCx.h — verbatim Velic kernel (_Velic_solCx), self-contained C
  • AnalyticSolCx.c / .h — thin vec2-shape adapters
  • analytic.pyx — raw wrappers (AnalyticSolCx_velocity_x/y, _pressure, _viscosity) plus a uw.function.analytic.SolCx(mesh, eta_A, eta_B, x_c, n) convenience exposing fn_velocity/pressure/viscosity/bodyforce and velocity_error()
  • setup.py — builds the new C sources into the analytic extension
  • tests/test_1015_analytic_solcx.py — binding + convergence-to-analytic (~4th order, rel 5.5e-7 at res 64)

Notes

  • The body-force sign matches UW3's momentum convention (the UW2 docs quote the opposite sign).
  • fn_velocity/fn_pressure wrap the compiled kernel and are evaluated point-wise via evalf (not lambdify-able) — use velocity_error(); fn_bodyforce/fn_viscosity are elementary and compile through the normal JIT path.

Follow-up

Finalise the full analytic suite (SolKx, SolKz, SolDB2d/3d, SolA/B/H, ...) by the same mechanical pattern, each with a convergence test. Tracked in planning.

Underworld development team with AI support from Claude Code

Adds the classic Velic "SolCx" exact solution (free-slip, discontinuous
viscosity step eta_A|eta_B at x_c, density forcing (0, cos(pi x) sin(n pi z)))
as an analytic validation reference. Mirrors the existing AnalyticSolNL port:
- solCx.c / solCx.h: verbatim Velic kernel (_Velic_solCx), pure self-contained C
- AnalyticSolCx.c / .h: thin vec2-shape adapters around the kernel
- analytic.pyx: AnalyticSolCx_velocity_x/y, _pressure, _viscosity wrappers
(JIT-printable + evalf), plus a uw.function.analytic.SolCx(mesh, eta_A,
eta_B, x_c, n) convenience exposing fn_velocity/pressure/viscosity/bodyforce
and velocity_error(). Body-force sign matches UW3's momentum convention.
- setup.py: build the new C sources into the analytic extension
- tests/test_1015_analytic_solcx.py: binding + convergence-to-analytic
(~4th order, rel 5.5e-7 at res 64)
This establishes the reusable pattern for porting the remaining Velic
solutions (SolKx, SolKz, SolDB, SolA/B/H, ...). It already caught a real
trap: a direct ksponly+lu solve silently mis-solves the singular saddle on
free-slip / pressure-nullspace problems; only the exact solution exposes it.
Underworld development team with AI support from Claude Code
CopilotAI review requested due to automatic review settings June 7, 2026 14:46

CopilotAI 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.

Pull request overview

Ports the classic Velic SolCx discontinuous-viscosity Stokes benchmark into underworld3.function.analytic by adding the upstream C kernel, thin UW3 adapters, Cython/SymPy wrappers, and a regression test that solves Stokes and checks convergence against the analytic solution.

Changes:

  • Adds the verbatim Velic SolCx kernel (solCx.c/.h) and UW3 adapter layer (AnalyticSolCx.c/.h).
  • Extends underworld3.function.analytic (Cython) with SolCx wrappers and a Python convenience class (SolCx) exposing fields and an error metric.
  • Introduces a new test validating binding non-triviality and solver convergence (tests/test_1015_analytic_solcx.py).

Reviewed changes

Copilot reviewed 7 out of 7 changed files in this pull request and generated 8 comments.

Show a summary per file
FileDescription
tests/test_1015_analytic_solcx.pyNew regression test for SolCx binding + Stokes convergence vs analytic
src/underworld3/function/solCx.hDeclares Velic SolCx kernel entry points
src/underworld3/function/solCx.cAdds the (verbatim) Velic SolCx analytic kernel implementation
src/underworld3/function/AnalyticSolCx.hDeclares UW3 adapter API for SolCx in vec2/scalar form
src/underworld3/function/AnalyticSolCx.cImplements thin adapters around _Velic_solCx plus viscosity step helper
src/underworld3/function/analytic.pyxAdds SolCx Cython bindings and Python convenience wrapper (SolCx)
setup.pyIncludes new SolCx C sources in the analytic extension build

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

self.fn_velocity = _sp.Matrix([AnalyticSolCx_velocity_x(*p, x, y),
AnalyticSolCx_velocity_y(*p, x, y)])
self.fn_pressure = AnalyticSolCx_pressure(*p, x, y)
self.fn_viscosity = _sp.Piecewise((self.eta_B, x > self.x_c), (self.eta_A, True))

# ----------------------------------------------------------------------------
# SolCx: viscosity step in x (eta_A | eta_B at x_c), trigonometric density
# forcing f = (0, -cos(pi x) sin(n pi z)) on a unit box with free-slip walls.
Comment on lines +11 to +13
isoviscous step in x (eta_A for x<x_c, eta_B for x>x_c), trigonometric
density forcing f = (0, -cos(pi x) sin(n pi z)) on a unit box with
free-slip walls. The heavy lifting is the verbatim Velic kernel in solCx.c
Comment threadsrc/underworld3/function/solCx.h Outdated
const double pos[],
double _eta_A, double _eta_B,
double _x_c, int _n,
double vel[], double* presssure,
Comment threadsrc/underworld3/function/solCx.h Outdated
const double pos[],
double _eta_A, double _eta_B, /* Input parameters: density, viscosity A, viscosity B */
double _x_c, int _n, /* Input parameters: viscosity jump location, wavenumber in x */
double vel[], double* presssure,
Comment threadsrc/underworld3/function/solCx.h Outdated
const double pos[],
double _eta_A, double _eta_B, /* Input parameters: density, viscosity A, viscosity B */
double _x_c, int _n, /* Input parameters: viscosity jump location, wavenumber in x */
double vel[], double* presssure,
self.eta_A = float(eta_A)
self.eta_B = float(eta_B)
self.x_c = float(x_c)
self.n = int(n)
self.eta_B = float(eta_B)
self.x_c = float(x_c)
self.n = int(n)
x, y = mesh.X
…n, typo
- fn_viscosity tie-break at x==x_c now matches the compiled SolCx_viscosity
(eta_A for x<x_c, eta_B for x>=x_c) — consistent PDE coefficients vs kernel
- correct the body-force sign in the comment block and AnalyticSolCx.h header
to the implemented +cos (UW3 convention; note UW2 docs quote -cos)
- SolCx.__init__: validate mesh.dim==2 and n a positive integer (fail fast;
n=0/negative is singular in the Velic kernel)
- solCx.h/.c: fix 'presssure' -> 'pressure' typo in the prototypes/defs
Underworld development team with AI support from Claude Code
@lmoresi
lmoresi merged commit d5afb4d into developmentJun 8, 2026
1 check passed
@lmoresi
lmoresi deleted the bugfix/analytic-solutions branch June 13, 2026 00:47
Sign up for freeto 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.

2 participants

@lmoresi