diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py index fb03e1b8..2ea99fa9 100644 --- a/src/underworld3/meshing/spherical.py +++ b/src/underworld3/meshing/spherical.py @@ -571,105 +571,118 @@ class regions(Enum): else: uw_filename = filename - # Check if r_i is greater than 0 if radiusInner <= 0: raise ValueError("The inner radius must be greater than 0.") + if not radiusInner < radiusInternal < radiusOuter: + raise ValueError( + "SphericalShellInternalBoundary requires " + "radiusInner < radiusInternal < radiusOuter." + ) if uw.mpi.rank == 0: gmsh.initialize() gmsh.option.setNumber("General.Verbosity", gmsh_verbosity) gmsh.model.add("SphereShell_with_Internal_Surface") - # Create three concentric spheres and use OCC fragment to split - # into two non-overlapping shell volumes sharing the internal surface - ball_outer = gmsh.model.occ.addSphere(0, 0, 0, radiusOuter) - ball_internal = gmsh.model.occ.addSphere(0, 0, 0, radiusInternal) - ball_inner = gmsh.model.occ.addSphere(0, 0, 0, radiusInner) + # Create the spherical shell volume. + outer = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusOuter) + inner = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInner) + gmsh.model.occ.cut( + [(3, outer)], + [(3, inner)], + removeObject=True, + removeTool=True, + ) - # Fragment creates non-overlapping pieces from the boolean intersection - out_dimtags, out_map = gmsh.model.occ.fragment( - [(3, ball_outer)], - [(3, ball_internal), (3, ball_inner)], + # Create an internal shell only to obtain a clean spherical surface at + # radiusInternal. That surface is embedded into the shell volume below; + # the duplicate volume and duplicate lower surface are removed before + # meshing. + internal = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInternal) + inner_copy = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInner) + gmsh.model.occ.cut( + [(3, internal)], + [(3, inner_copy)], + removeObject=True, + removeTool=True, ) gmsh.model.occ.synchronize() gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cellSize) - # Identify volumes and surfaces by bounding box - # For a sphere, bbox diagonal = sqrt(3) * radius - volumes = gmsh.model.getEntities(3) - surfaces = gmsh.model.getEntities(2) - def bbox_radius(dimtag): - """Estimate the sphere radius from a bounding box diagonal.""" + """Estimate a concentric sphere radius from the bounding box.""" bb = gmsh.model.get_bounding_box(dimtag[0], dimtag[1]) return np.sqrt(bb[3]**2 + bb[4]**2 + bb[5]**2) / np.sqrt(3.0) - inner_vols = [] - outer_vols = [] - solid_ball_vols = [] # r < radiusInner — to be removed - - for vol in volumes: - r_est = bbox_radius(vol) - if np.isclose(r_est, radiusInner, atol=cellSize): - solid_ball_vols.append(vol) - elif np.isclose(r_est, radiusInternal, atol=cellSize): - inner_vols.append(vol) - elif np.isclose(r_est, radiusOuter, atol=cellSize): - outer_vols.append(vol) - - # Remove the solid inner ball (r < radiusInner) - if solid_ball_vols: - gmsh.model.occ.remove(solid_ball_vols, recursive=True) - gmsh.model.occ.synchronize() - - # Re-query after removal volumes = gmsh.model.getEntities(3) - surfaces = gmsh.model.getEntities(2) + shell_vols = [ + vol + for vol in volumes + if np.isclose(bbox_radius(vol), radiusOuter, atol=cellSize * 0.5) + ] + duplicate_vols = [ + vol + for vol in volumes + if np.isclose(bbox_radius(vol), radiusInternal, atol=cellSize * 0.5) + ] + + if len(shell_vols) != 1: + raise RuntimeError( + "Could not identify the spherical-shell volume while building " + "SphericalShellInternalBoundary." + ) - # Classify surfaces by bounding box radius - for surface in surfaces: + shell_vol = shell_vols[0] + shell_boundary = { + dimtag[1] + for dimtag in gmsh.model.getBoundary([shell_vol], oriented=False, recursive=False) + if dimtag[0] == 2 + } + + lower_surface_tags = [] + internal_surface_tags = [] + upper_surface_tags = [] + duplicate_lower_tags = [] + + for surface in gmsh.model.getEntities(2): + surface_tag = surface[1] r_est = bbox_radius(surface) if np.isclose(r_est, radiusInner, atol=cellSize * 0.5): - gmsh.model.addPhysicalGroup( - surface[0], [surface[1]], - boundaries.Lower.value, name=boundaries.Lower.name, - ) + if surface_tag in shell_boundary: + lower_surface_tags.append(surface_tag) + else: + duplicate_lower_tags.append(surface_tag) elif np.isclose(r_est, radiusOuter, atol=cellSize * 0.5): - gmsh.model.addPhysicalGroup( - surface[0], [surface[1]], - boundaries.Upper.value, name=boundaries.Upper.name, - ) + upper_surface_tags.append(surface_tag) elif np.isclose(r_est, radiusInternal, atol=cellSize * 0.5): - gmsh.model.addPhysicalGroup( - surface[0], [surface[1]], - boundaries.Internal.value, name=boundaries.Internal.name, - ) - - # Classify remaining volumes into Inner and Outer - inner_vol_tags = [v[1] for v in inner_vols if v not in solid_ball_vols] - outer_vol_tags = [v[1] for v in outer_vols] - # Re-classify from current volumes in case tags changed after removal - inner_vol_tags = [] - outer_vol_tags = [] - for vol in volumes: - r_est = bbox_radius(vol) - if r_est < radiusInternal + cellSize * 0.5: - inner_vol_tags.append(vol[1]) - else: - outer_vol_tags.append(vol[1]) - - # Region physical groups - if inner_vol_tags: - gmsh.model.addPhysicalGroup(3, inner_vol_tags, - regions.Inner.value, name=regions.Inner.name) - if outer_vol_tags: - gmsh.model.addPhysicalGroup(3, outer_vol_tags, - regions.Outer.value, name=regions.Outer.name) - - # Combined elements group - all_vol_tags = inner_vol_tags + outer_vol_tags - gmsh.model.addPhysicalGroup(3, all_vol_tags, 99999, "Elements") + internal_surface_tags.append(surface_tag) + + if not lower_surface_tags or not upper_surface_tags or not internal_surface_tags: + raise RuntimeError( + "Could not identify Lower, Internal, and Upper spherical surfaces " + "while building SphericalShellInternalBoundary." + ) + + gmsh.model.mesh.embed(2, internal_surface_tags, shell_vol[0], shell_vol[1]) + + remove_dimtags = duplicate_vols + [(2, tag) for tag in duplicate_lower_tags] + if remove_dimtags: + gmsh.model.remove_entities(remove_dimtags, recursive=False) + gmsh.model.occ.remove(remove_dimtags, recursive=False) + gmsh.model.occ.synchronize() + + gmsh.model.addPhysicalGroup( + 2, lower_surface_tags, boundaries.Lower.value, name=boundaries.Lower.name + ) + gmsh.model.addPhysicalGroup( + 2, internal_surface_tags, boundaries.Internal.value, name=boundaries.Internal.name + ) + gmsh.model.addPhysicalGroup( + 2, upper_surface_tags, boundaries.Upper.value, name=boundaries.Upper.name + ) + + gmsh.model.addPhysicalGroup(shell_vol[0], [shell_vol[1]], 99999, "Elements") gmsh.model.mesh.generate(3) gmsh.write(uw_filename) diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py index 0cd9e18d..1721db7d 100644 --- a/tests/test_0502_boundary_integrals.py +++ b/tests/test_0502_boundary_integrals.py @@ -293,6 +293,52 @@ def test_bd_integral_annulus_internal_normal_tangential(): assert abs(value) < 0.05, f"Expected ~0, got {value}" +# --- Spherical shell internal boundary tests --- + +from underworld3.meshing import SphericalShellInternalBoundary + +_R_SHELL_INNER = 0.55 +_R_SHELL_INTERNAL = 0.775 +_R_SHELL_OUTER = 1.0 +_mesh_spherical_internal = None + + +def _get_spherical_internal_mesh(): + global _mesh_spherical_internal + if _mesh_spherical_internal is None: + _mesh_spherical_internal = SphericalShellInternalBoundary( + radiusOuter=_R_SHELL_OUTER, + radiusInternal=_R_SHELL_INTERNAL, + radiusInner=_R_SHELL_INNER, + cellSize=0.25, + degree=1, + qdegree=2, + ) + uw.discretisation.MeshVariable( + "T_spherical_internal", _mesh_spherical_internal, 1, degree=1 + ) + return _mesh_spherical_internal + + +def test_bd_integral_spherical_internal_boundary_areas(): + """SphericalShellInternalBoundary preserves Lower/Internal/Upper labels.""" + + mesh_spherical = _get_spherical_internal_mesh() + expected_areas = { + "Lower": 4.0 * np.pi * _R_SHELL_INNER**2, + "Internal": 4.0 * np.pi * _R_SHELL_INTERNAL**2, + "Upper": 4.0 * np.pi * _R_SHELL_OUTER**2, + } + + for boundary, expected in expected_areas.items(): + value = uw.maths.BdIntegral(mesh_spherical, fn=1.0, boundary=boundary).evaluate() + relative_error = abs(value - expected) / expected + assert relative_error < 0.06, ( + f"{boundary} area should be close to {expected:.4f}; " + f"got {value:.4f} (relative error {relative_error:.3f})" + ) + + def _build_spherical_shell_for_integrals(): from underworld3.meshing import SphericalShell diff --git a/tests/test_1064_constrained_spherical_shell_response.py b/tests/test_1064_constrained_spherical_shell_response.py new file mode 100644 index 00000000..3071047d --- /dev/null +++ b/tests/test_1064_constrained_spherical_shell_response.py @@ -0,0 +1,184 @@ +"""3-D spherical-shell constrained free-slip response regression. + +This test records two facts exposed by the Zhong et al. (2008)-style benchmark: + +* with monolithic LU, ``Stokes_Constrained`` and Nitsche free slip solve the same + internal-boundary Stokes response to within a tight tolerance; +* the practical fast grouped-Schur constrained path currently does not reproduce + the Zhong velocity response and remains an expected failure. + +Run: + pixi run -e amr-dev pytest -q tests/test_1064_constrained_spherical_shell_response.py +""" + +from functools import cache + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_3, pytest.mark.slow, pytest.mark.tier_c] + +RADIUS_INNER = 0.55 +RADIUS_INTERNAL = 0.775 +RADIUS_OUTER = 1.0 +CELL_SIZE = 1.0 / 8.0 +HARMONIC_DEGREE = 2 +NITSCHE_GAMMA = 10.0 + +ZHONG_SURFACE_VELOCITY = 1.006e-2 +ZHONG_CMB_VELOCITY = 1.186e-2 + + +@cache +def solve_response(method, solver_mode): + mesh = uw.meshing.SphericalShellInternalBoundary( + radiusOuter=RADIUS_OUTER, + radiusInternal=RADIUS_INTERNAL, + radiusInner=RADIUS_INNER, + cellSize=CELL_SIZE, + qdegree=2, + degree=1, + ) + + velocity = uw.discretisation.MeshVariable( + f"U_{method}_{solver_mode}", + mesh, + mesh.dim, + degree=2, + vtype=uw.VarType.VECTOR, + ) + pressure = uw.discretisation.MeshVariable( + f"P_{method}_{solver_mode}", + mesh, + 1, + degree=1, + continuous=True, + ) + + theta = mesh.CoordinateSystem.xR[1] + unit_r = mesh.CoordinateSystem.unit_e_0 + y_l0 = sympy.assoc_legendre(HARMONIC_DEGREE, 0, sympy.cos(theta)) + harmonic_norm = 4.0 * np.pi / (2 * HARMONIC_DEGREE + 1) + + if method == "constrained": + stokes = uw.systems.Stokes_Constrained( + mesh, + velocityField=velocity, + pressureField=pressure, + ) + elif method == "nitsche": + stokes = uw.systems.Stokes( + mesh, + velocityField=velocity, + pressureField=pressure, + ) + else: + raise ValueError(f"Unknown response method {method!r}") + + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, 0.0, 0.0]) + stokes.add_natural_bc(y_l0 * unit_r, mesh.boundaries.Internal.name) + + if method == "nitsche": + stokes.add_nitsche_bc("Upper", normal=unit_r, gamma=NITSCHE_GAMMA) + stokes.add_nitsche_bc("Lower", normal=-unit_r, gamma=NITSCHE_GAMMA) + else: + stokes.add_constraint_bc( + "Upper", + g=0.0, + normal=unit_r, + augmentation_base=1.0e4, + degree=2, + ) + stokes.add_constraint_bc( + "Lower", + g=0.0, + normal=-unit_r, + augmentation_base=1.0e4, + degree=2, + ) + + stokes.petsc_use_nullspace = True + stokes.tolerance = 1.0e-7 + stokes.petsc_options["snes_type"] = "ksponly" + + if solver_mode == "monolithic": + stokes.petsc_options["ksp_type"] = "preonly" + stokes.petsc_options["pc_type"] = "lu" + stokes.petsc_options["pc_factor_mat_solver_type"] = "mumps" + stokes.petsc_options["pc_use_amat"] = None + elif solver_mode == "fast_schur": + if method != "constrained": + raise ValueError("fast_schur mode is only defined for constrained runs") + stokes.petsc_options["pc_fieldsplit_schur_precondition"] = "selfp" + stokes.petsc_options["fieldsplit_1_ksp_type"] = "preonly" + stokes.petsc_options["fieldsplit_1_pc_type"] = "gasm" + else: + raise ValueError(f"Unknown solver mode {solver_mode!r}") + + stokes.solve() + + horizontal_v2 = velocity.sym.dot(velocity.sym) - velocity.sym.dot(unit_r) ** 2 + + surface_velocity = np.sqrt( + uw.maths.BdIntegral(mesh, horizontal_v2, boundary="Upper").evaluate() + / ( + RADIUS_OUTER**2 + * HARMONIC_DEGREE + * (HARMONIC_DEGREE + 1) + * harmonic_norm + ) + ) + cmb_velocity = np.sqrt( + uw.maths.BdIntegral(mesh, horizontal_v2, boundary="Lower").evaluate() + / ( + RADIUS_INNER**2 + * HARMONIC_DEGREE + * (HARMONIC_DEGREE + 1) + * harmonic_norm + ) + ) + + return ( + float(surface_velocity), + float(cmb_velocity), + int(stokes.snes.getConvergedReason()), + ) + + +def test_monolithic_constrained_matches_monolithic_nitsche_response(): + nitsche_surface, nitsche_cmb, nitsche_reason = solve_response( + "nitsche", + "monolithic", + ) + constrained_surface, constrained_cmb, constrained_reason = solve_response( + "constrained", + "monolithic", + ) + + assert nitsche_reason > 0 + assert constrained_reason > 0 + assert abs(constrained_surface - nitsche_surface) / nitsche_surface < 0.01 + assert abs(constrained_cmb - nitsche_cmb) / nitsche_cmb < 0.01 + + +@pytest.mark.xfail( + reason=( + "Known fast grouped-Schur constrained response failure for the " + "3-D SphericalShellInternalBoundary Zhong-style load." + ), + strict=True, +) +def test_fast_schur_constrained_matches_zhong_velocity_response(): + surface_velocity, cmb_velocity, snes_reason = solve_response( + "constrained", + "fast_schur", + ) + + assert snes_reason > 0 + assert abs(surface_velocity - ZHONG_SURFACE_VELOCITY) / ZHONG_SURFACE_VELOCITY < 0.05 + assert abs(cmb_velocity - ZHONG_CMB_VELOCITY) / ZHONG_CMB_VELOCITY < 0.05