diff --git a/src/underworld3/function/AnalyticSolCx.c b/src/underworld3/function/AnalyticSolCx.c index 6d59b1a14..a17560fc4 100644 --- a/src/underworld3/function/AnalyticSolCx.c +++ b/src/underworld3/function/AnalyticSolCx.c @@ -43,3 +43,17 @@ double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, const dou (void) n; (void) z; return ( x < x_c ) ? eta_A : eta_B; } + +tensor2 SolCx_stress( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) +{ + double pos[2]; + double total_stress[3]; /* [0]=sigma_xx, [1]=sigma_zz, [2]=sigma_xz (full Cauchy) */ + tensor2 out; + pos[0] = x; + pos[1] = z; + _Velic_solCx( pos, eta_A, eta_B, x_c, n, NULL, NULL, total_stress, NULL ); + out.xx = total_stress[0]; + out.zz = total_stress[1]; + out.xz = total_stress[2]; + return out; +} diff --git a/src/underworld3/function/AnalyticSolCx.h b/src/underworld3/function/AnalyticSolCx.h index aafd0c5b3..841e9d6d1 100644 --- a/src/underworld3/function/AnalyticSolCx.h +++ b/src/underworld3/function/AnalyticSolCx.h @@ -24,8 +24,9 @@ the result depends solely on the arguments, so redundant calls (e.g. the x and z components evaluated separately) can be collapsed by the compiler. */ -vec2 SolCx_velocity( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); -double SolCx_pressure( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); -double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); +vec2 SolCx_velocity( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); +double SolCx_pressure( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); +double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); +tensor2 SolCx_stress( double eta_A, double eta_B, double x_c, int n, const double x, const double z ) __attribute__((const)); #endif diff --git a/src/underworld3/function/analytic.pyx b/src/underworld3/function/analytic.pyx index dbc237516..18ad20018 100644 --- a/src/underworld3/function/analytic.pyx +++ b/src/underworld3/function/analytic.pyx @@ -12,15 +12,20 @@ cdef extern from "AnalyticSolNL.h" nogil: ctypedef struct vec2: double x double z + ctypedef struct tensor2: + double xx + double zz + double xz vec2 SolNL_velocity( double eta0, unsigned n, double r, double x, double y ) vec2 SolNL_bodyforce( double eta0, unsigned n, double r, double x, double y ) double SolNL_viscosity( double eta0, unsigned n, double r, double x, double y ) cdef extern from "AnalyticSolCx.h" nogil: - # vec2 already declared above (shared typedef via AnalyticSolNL.h) - vec2 SolCx_velocity( double eta_A, double eta_B, double x_c, int n, double x, double z ) - double SolCx_pressure( double eta_A, double eta_B, double x_c, int n, double x, double z ) - double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, double x, double z ) + # vec2 / tensor2 already declared above (shared typedefs via AnalyticSolNL.h) + vec2 SolCx_velocity( double eta_A, double eta_B, double x_c, int n, double x, double z ) + double SolCx_pressure( double eta_A, double eta_B, double x_c, int n, double x, double z ) + double SolCx_viscosity( double eta_A, double eta_B, double x_c, int n, double x, double z ) + tensor2 SolCx_stress( double eta_A, double eta_B, double x_c, int n, double x, double z ) class sympy_function_printable(sympy.Function): """ @@ -131,6 +136,25 @@ class AnalyticSolCx_viscosity(AnalyticSolCx_base): from sympy import sympify return sympify(SolCx_viscosity( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] )) +# Exact total (Cauchy) stress components sigma_ij from the Velic kernel's +# total_stress array: xx, zz (= yy in the x-z plane), xz. Dynamic topography on +# the top boundary (n = y_hat) is -n.sigma.n = -sigma_zz. +class AnalyticSolCx_stress_xx(AnalyticSolCx_base): + _printstr = "SolCx_stress({}).xx" + def _eval_evalf(self,prec): + from sympy import sympify + return sympify(SolCx_stress( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ).xx) +class AnalyticSolCx_stress_yy(AnalyticSolCx_base): + _printstr = "SolCx_stress({}).zz" + def _eval_evalf(self,prec): + from sympy import sympify + return sympify(SolCx_stress( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ).zz) +class AnalyticSolCx_stress_xy(AnalyticSolCx_base): + _printstr = "SolCx_stress({}).xz" + def _eval_evalf(self,prec): + from sympy import sympify + return sympify(SolCx_stress( self.args[0],self.args[1],self.args[2],self.args[3],self.args[4],self.args[5] ).xz) + class SolCx: r"""Velic *SolCx* analytic Stokes solution — the canonical free-slip, @@ -175,6 +199,10 @@ class SolCx: 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) + # exact total (Cauchy) stress components (sigma_xx, sigma_yy, sigma_xy) + self.fn_stress_xx = AnalyticSolCx_stress_xx(*p, x, y) + self.fn_stress_yy = AnalyticSolCx_stress_yy(*p, x, y) + self.fn_stress_xy = AnalyticSolCx_stress_xy(*p, x, y) # tie-break at x==x_c matches the compiled SolCx_viscosity (eta_B for x>=x_c) self.fn_viscosity = _sp.Piecewise((self.eta_A, x < self.x_c), (self.eta_B, True)) self.fn_bodyforce = _sp.Matrix([_sp.Integer(0), @@ -195,4 +223,26 @@ class SolCx: """Relative L2 error of a velocity MeshVariable against the analytic.""" import numpy as _np ua = self.evaluate_velocity(velocity_var.coords) - return float(_np.linalg.norm(velocity_var.data - ua) / _np.linalg.norm(ua)) \ No newline at end of file + return float(_np.linalg.norm(velocity_var.data - ua) / _np.linalg.norm(ua)) + + def evaluate_stress(self, coords): + """Exact total (Cauchy) stress at ``coords`` (N×2 array), via the + compiled kernel. Returns an (N, 3) array of (sigma_xx, sigma_yy, + sigma_xy).""" + import numpy as _np + cdef Py_ssize_t i + cdef double xi, yi + cdef tensor2 s + out = _np.empty((len(coords), 3), dtype=_np.float64) + for i in range(len(coords)): + xi = float(coords[i, 0]); yi = float(coords[i, 1]) + s = SolCx_stress(self.eta_A, self.eta_B, self.x_c, self.n, xi, yi) + out[i, 0] = s.xx + out[i, 1] = s.zz + out[i, 2] = s.xz + return out + + def topography_top(self, coords): + """Exact dynamic topography -n.sigma.n on the top boundary (n = y_hat), + i.e. -sigma_yy, at ``coords`` (N×2 array). Returns a length-N array.""" + return -self.evaluate_stress(coords)[:, 1] \ No newline at end of file diff --git a/tests/test_1015_analytic_solcx.py b/tests/test_1015_analytic_solcx.py index b9d623d63..4ba67c102 100644 --- a/tests/test_1015_analytic_solcx.py +++ b/tests/test_1015_analytic_solcx.py @@ -56,3 +56,32 @@ def test_solcx_stokes_converges_to_analytic(): e32 = _solve_solcx(32) assert e32 < e16, f"no convergence: e16={e16:.2e} e32={e32:.2e}" assert e32 < 1.0e-3, f"poor match to analytic: e32={e32:.2e}" + + +def test_solcx_stress_binding(): + """The exact total-stress kernel binding returns finite, non-trivial values.""" + pts = [(0.25, 1.0), (0.7, 0.2), (0.9, 0.6), (0.15, 0.4)] + sxx = [float(A.AnalyticSolCx_stress_xx(1.0, 1.0e6, 0.5, 1, x, y).evalf()) for x, y in pts] + syy = [float(A.AnalyticSolCx_stress_yy(1.0, 1.0e6, 0.5, 1, x, y).evalf()) for x, y in pts] + sxy = [float(A.AnalyticSolCx_stress_xy(1.0, 1.0e6, 0.5, 1, x, y).evalf()) for x, y in pts] + for s in sxx + syy + sxy: + assert np.isfinite(s) + # the stress field is non-trivial (the dynamic-topography signal on the top is -syy) + assert max(abs(s) for s in syy) > 1.0e-3 + assert max(abs(s) for s in sxx) > 1.0e-3 + + +def test_solcx_stress_pressure_consistency(): + """Convention-independent check: in 2D sigma = -p I + tau with tau traceless, + so the mean normal stress (sigma_xx + sigma_yy)/2 must equal -p at every + point. This validates the new stress wrapper against the kernel's own + pressure without assuming any UW3 sign convention.""" + pts = [(0.2, 0.3), (0.7, 0.8), (0.5, 0.5), (0.9, 0.1), (0.35, 0.65)] + for xi, yi in pts: + sxx = float(A.AnalyticSolCx_stress_xx(1.0, 1.0e6, 0.5, 1, xi, yi).evalf()) + syy = float(A.AnalyticSolCx_stress_yy(1.0, 1.0e6, 0.5, 1, xi, yi).evalf()) + p = float(A.AnalyticSolCx_pressure(1.0, 1.0e6, 0.5, 1, xi, yi).evalf()) + mean_normal = 0.5 * (sxx + syy) + assert abs(mean_normal + p) < 1.0e-10, ( + f"(sxx+syy)/2 != -p at ({xi},{yi}): {mean_normal:.3e} vs {-p:.3e}" + )