Skip to content

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

Open
lmoresi wants to merge 1 commit into
developmentfrom
bugfix/analytic-solutions
Open

Port SolCx analytic Stokes solution (uw.function.analytic.SolCx)#223
lmoresi wants to merge 1 commit into
developmentfrom
bugfix/analytic-solutions

Conversation

@lmoresi
Copy link
Copy Markdown
Member

@lmoresi lmoresi commented Jun 7, 2026

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
Copilot AI review requested due to automatic review settings June 7, 2026 14:46
Copy link
Copy Markdown
Contributor

Copilot AI left a comment

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
File Description
tests/test_1015_analytic_solcx.py New regression test for SolCx binding + Stokes convergence vs analytic
src/underworld3/function/solCx.h Declares Velic SolCx kernel entry points
src/underworld3/function/solCx.c Adds the (verbatim) Velic SolCx analytic kernel implementation
src/underworld3/function/AnalyticSolCx.h Declares UW3 adapter API for SolCx in vec2/scalar form
src/underworld3/function/AnalyticSolCx.c Implements thin adapters around _Velic_solCx plus viscosity step helper
src/underworld3/function/analytic.pyx Adds SolCx Cython bindings and Python convenience wrapper (SolCx)
setup.py Includes 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
const double pos[],
double _eta_A, double _eta_B,
double _x_c, int _n,
double vel[], double* presssure,
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,
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
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.

2 participants