{ "cells": [ { "cell_type": "markdown", "id": "d0fcf68c", "metadata": {}, "source": [ "# Oseen MHM: write momentum, incompressibility and interface forms\n", "\n", "This notebook follows **meshes → local equations → global equations → assemble\n", "→ solve → velocity and pressure**. A user writes mixed UFL forms directly;\n", "PyMHM supplies mesh-associated trace numbering, orientation and condensation.\n", "\n", "We first solve an independently manufactured affine problem exactly in stable\n", "Taylor–Hood P2/P1 local spaces. Its polynomial exactness is a correctness check,\n", "not a convergence rate. The final section reads the separately attributed smooth\n", "stabilized P3/P3 refinement family and plots its measured asymptotic rates.\n", "\n", "Run with the locked `introduction` profile. No optional native cell is silently\n", "skipped.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "companion-bootstrap", "metadata": { "tags": [ "pymhm-companion-bootstrap" ] }, "outputs": [], "source": [ "from pathlib import Path\n", "import os\n", "import sys\n", "from pymhm.io.workspace import workspace_from_archive\n", "\n", "# Download checksum-verified support files; extraction does not execute them.\n", "COMPANION_URL = \"https://ipes-lncc.github.io/pymhm/downloads/b0dc88575f8a4f7dec3cee89030544a52c7e03ce7a302aaad9b86191d4c8f29f/oseen_variational-companion.zip\"\n", "COMPANION_SHA256 = \"b0dc88575f8a4f7dec3cee89030544a52c7e03ce7a302aaad9b86191d4c8f29f\"\n", "WORKSPACE = Path(\n", " os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd() / \".pymhm-companions\" / COMPANION_SHA256)\n", ")\n", "ROOT = workspace_from_archive(COMPANION_URL, sha256=COMPANION_SHA256, directory=WORKSPACE)\n", "os.environ[\"PYMHM_WORKSPACE\"] = str(ROOT)\n", "if str(ROOT) not in sys.path:\n", " sys.path.insert(0, str(ROOT))\n", "\n", "# Prepare only this notebook's declared inputs, without executing its equations.\n", "from scripts.notebook_reproduction import notebook_workspace\n", "\n", "ROOT = notebook_workspace(\"flow/oseen_variational.ipynb\", directory=ROOT)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dc077ef0", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "import os\n", "import sys\n", "ROOT = Path(os.environ.get(\"PYMHM_WORKSPACE\", Path.cwd())).resolve()\n", "if not (ROOT / \"examples\").is_dir():\n", " ROOT = next((p for p in ROOT.parents if (p / \"examples\").is_dir()), ROOT)\n", "if str(ROOT) not in sys.path: sys.path.insert(0, str(ROOT))\n", "OUTPUT = ROOT / \"build\" / \"native-teaching\"\n", "OUTPUT.mkdir(parents=True, exist_ok=True)\n", "print(\"Workspace:\", ROOT)\n" ] }, { "cell_type": "markdown", "id": "f7421319", "metadata": {}, "source": [ "## 1. Prescribe the physical solution before defining a matrix\n", "\n", "Take viscosity one, resistance one, and constant prescribed convection\n", "$\\beta=(2,3)$. The nonconservative Oseen equation is\n", "\n", "$$\n", "-\\Delta u+(\\beta\\cdot\\nabla)u+u+\\nabla p=f,\\qquad \\nabla\\cdot u=0.\n", "$$\n", "\n", "Choose\n", "\n", "$$\n", "u=(1-x+2y,-2-x+y),\\quad p=x+2y-0.8,\\quad f=u+(5,3).\n", "$$\n", "\n", "The gradient of velocity is constant, its divergence vanishes, and its\n", "Laplacian is zero. Hence $(\\beta\\cdot\\nabla)u=(4,1)$ and $\\nabla p=(1,2)$.\n", "The physical pressure mean is 0.7. Prescribe the exact velocity on every\n", "exterior macroface; pressure is fixed by its global mean, not by one mean per\n", "local problem.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "aadd7c4b", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import ufl\n", "import basix\n", "import basix.ufl\n", "from dolfinx import fem\n", "from pymhm import (\n", " TriangleMesh, FaceSpace, SkeletonSpace, Equation, MeshHierarchy,\n", " LocalContext, columns, compile_form, bind_interface, bind_problem, assemble, solve,\n", ")\n", "\n", "def exact_velocity(x):\n", " \"\"\"Return the solenoidal affine physical velocity.\"\"\"\n", " return np.column_stack((1-x[:,0]+2*x[:,1], -2-x[:,0]+x[:,1]))\n", "\n", "def exact_pressure(x):\n", " \"\"\"Return pressure with declared physical unit-square mean 0.7.\"\"\"\n", " return x[:,0]+2*x[:,1]-0.8\n", "\n", "macro=TriangleMesh.unit_square(2)\n", "skeleton=SkeletonSpace(macro,tuple(FaceSpace.uniform(1) for _ in macro.faces),2)\n", "hierarchy=MeshHierarchy(macro,tuple(macro.submesh(i,2) for i in range(len(macro.cells))))\n", "interface=bind_interface(skeleton,convention=\"normal\")\n", "\n", "\n", "taylor_hood_element=basix.ufl.mixed_element([\n", " basix.ufl.element(\"Lagrange\",\"triangle\",2,shape=(2,),lagrange_variant=basix.LagrangeVariant.equispaced),\n", " basix.ufl.element(\"Lagrange\",\"triangle\",1,lagrange_variant=basix.LagrangeVariant.equispaced),\n", "])\n" ] }, { "cell_type": "markdown", "id": "4d07e24c", "metadata": {}, "source": [ "## 2. Write the local mixed energy and both interface pairings\n", "\n", "Use the skew convection form\n", "\n", "$$\n", "\\begin{aligned}\n", "a_K((u,p),(v,q))={}&(\\nabla u,\\nabla v)_K+(u,v)_K\\\\\n", "&+\\tfrac12[(\\beta\\cdot\\nabla u,v)_K-(u,\\beta\\cdot\\nabla v)_K]\\\\\n", "&-(p,\\nabla\\cdot v)_K-(q,\\nabla\\cdot u)_K.\n", "\\end{aligned}\n", "$$\n", "\n", "Since $\\nabla\\cdot\\beta=0$, the effective reaction is one. The normal interface\n", "quantity is the **half-advection pseudotraction**\n", "\n", "$$\n", "\\lambda_K=(-\\nabla u+pI+\\tfrac12u\\otimes\\beta)n_K.\n", "$$\n", "\n", "The interface action is $+\\langle\\lambda_K,v\\rangle$; the global velocity pairing\n", "is $-\\langle\\mu_K,u\\rangle$. `trace_pairings` supplies the physical supports and\n", "`local.equations` applies the normal orientation once. No face permutation or\n", "orientation code is written here.\n", "\n", "Two translations are retained through physical velocity-integral moments.\n", "They are **coarse modes**, not null vectors when drag is positive. The local\n", "mixed pressure has no independently imposed mean. Named fields retain the\n", "native executed coefficient basis for later evaluation.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "855d20cf", "metadata": {}, "outputs": [], "source": [ "def oseen_volume_forms(W, domain):\n", " \"\"\"Write the common physical UFL operator, source and pressure integral.\"\"\"\n", " u,p=ufl.TrialFunctions(W)\n", " v,q=ufl.TestFunctions(W)\n", " x=ufl.SpatialCoordinate(domain)\n", " beta=ufl.as_vector((2.0,3.0))\n", " exact=ufl.as_vector((1-x[0]+2*x[1],-2-x[0]+x[1]))\n", " force=exact+ufl.as_vector((5.0,3.0))\n", " dx=ufl.Measure(\"dx\",domain=domain,metadata={\"quadrature_degree\":8})\n", " a=(ufl.inner(ufl.grad(u),ufl.grad(v))+ufl.inner(u,v)\n", " +0.5*(ufl.inner(ufl.grad(u)*beta,v)-ufl.inner(u,ufl.grad(v)*beta))\n", " -p*ufl.div(v)-q*ufl.div(u))*dx\n", " L=ufl.inner(force,v)*dx\n", " return a,L,q*dx\n", "\n", "def local_oseen(local: LocalContext):\n", " \"\"\"Declare one Taylor–Hood Oseen problem in its executed native basis.\"\"\"\n", " binding=local.native_space(taylor_hood_element)\n", " W=binding.space\n", " u,p=ufl.TrialFunctions(W)\n", " v,q=ufl.TestFunctions(W)\n", " a,L,mean_form=oseen_volume_forms(W,binding.mesh)\n", " dx=ufl.Measure(\"dx\",domain=binding.mesh,metadata={\"quadrature_degree\":8})\n", " b=local.trace_pairings(lambda phi,ds:ufl.inner(phi,v)*ds)\n", " c=local.trace_pairings(lambda phi,ds:-ufl.inner(phi,u)*ds,axis=\"rows\")\n", " _,velocity_dofs=W.sub(0).collapse()\n", " translations=np.zeros((binding.size,2))\n", " translations[velocity_dofs]=np.tile(np.eye(2),(len(velocity_dofs)//2,1))\n", " local.field(\"velocity\",binding,component=0)\n", " local.field(\"pressure\",binding,component=1)\n", " return local.equations(\n", " a=a,L=L,b=b,c=c,\n", " coarse_basis=translations,moments=columns(v[0]*dx,v[1]*dx),\n", " metadata={\"pressure_weights\":compile_form(q*dx)},\n", " )\n" ] }, { "cell_type": "markdown", "id": "b185fde8", "metadata": {}, "source": [ "## 3. Declare the global velocity boundary equation and pressure gauge\n", "\n", "Local `c` equations sum signed velocity continuity and exterior velocity\n", "moments. The additional global functional has the matching minus sign.\n", "The global physical pressure constraint is imposed **after assembly**, because\n", "its row depends on the reconstructed local pressure responses.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e69a12ac", "metadata": {}, "outputs": [], "source": [ "def global_equation(global_context):\n", " \"\"\"Insert the prescribed velocity moments in the declared global convention.\"\"\"\n", " boundary,_=global_context.boundary_data(exact_velocity,order=6)\n", " return Equation(0,-global_context.trace_load(boundary))\n", "\n", "problem=bind_problem(hierarchy,interface,local_oseen,\n", " global_equation=global_equation,retained=2)\n", "system=assemble(problem)\n", "pressure_weights=[record[\"pressure_weights\"] for record in system.local_metadata]\n", "gauge=system.mean_constraint(pressure_weights,0.7)\n", "solution=solve(system,constraints=(gauge,))\n", "velocity=solution.field(\"velocity\")\n", "pressure=solution.field(\"pressure\")\n" ] }, { "cell_type": "markdown", "id": "36aad409", "metadata": {}, "source": [ "## 4. Measure independent physical fields and original equations\n", "\n", "The exact affine functions belong to the declared local spaces. Check both\n", "fields at independent random physical points, and integrate the pressure mean\n", "through its assembled physical weights. Zero errors in a polynomial patch\n", "are roundoff controls and are not fitted as convergence slopes.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "af7dd598", "metadata": {}, "outputs": [], "source": [ "maximum={\"velocity\":0.0,\"pressure\":0.0}\n", "for vf,pf in zip(velocity,pressure,strict=True):\n", " bary=np.random.default_rng(127).dirichlet([1,1,1],size=40)\n", " vertices=vf.mesh.points[vf.mesh.cells[0]]\n", " points=bary@vertices\n", " maximum[\"velocity\"]=max(maximum[\"velocity\"],float(np.max(abs(vf.evaluate(points)-exact_velocity(points)))))\n", " maximum[\"pressure\"]=max(maximum[\"pressure\"],float(np.max(abs(pf.evaluate(points)-exact_pressure(points)))))\n", "mean=sum(float(w@values) for w,values in zip(pressure_weights,solution.fields,strict=True))\n", "assert max(maximum.values())<1e-10\n", "np.testing.assert_allclose(mean,0.7,rtol=1e-10,atol=1e-12)\n", "print({\"maximum_physical_errors\":maximum,\"pressure_mean\":mean,\n", " \"global_compatibility_relative\":solution.raw_residual})\n" ] }, { "cell_type": "markdown", "id": "64e23ef5", "metadata": {}, "source": [ "## 5. Independently assemble a fine classical Taylor–Hood reference\n", "\n", "The same written UFL volume operator is assembled on one globally conforming\n", "P2/P1 space, without local condensation. Strong exterior velocity and the same\n", "physical pressure mean make the reference a matching physical problem. Two\n", "fine meshes independently reproduce this affine solution; errors near roundoff\n", "are exactness controls and do not define a fitted rate.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "71fdfeb7", "metadata": {}, "outputs": [], "source": [ "from pymhm import MultiscaleProblem\n", "from pymhm.backends.spaces import bind_space\n", "from pymhm.postprocessing.fields import DiscreteField,FieldDefinition\n", "\n", "def classical_oseen(n):\n", " \"\"\"Assemble the same mixed forms globally with strong velocity and one pressure mean.\"\"\"\n", " fine=TriangleMesh.unit_square(n)\n", " binding=bind_space(fine,taylor_hood_element)\n", " try:\n", " W=binding.space\n", " a,L,pressure_mean=oseen_volume_forms(W,binding.mesh)\n", " A,F,m=compile_form(a),compile_form(L),compile_form(pressure_mean)\n", " V,vmap=W.sub(0).collapse()\n", " coordinates=V.tabulate_dof_coordinates()[:,:2]\n", " boundary=np.flatnonzero(np.any(np.isclose(coordinates,0)|np.isclose(coordinates,1),axis=1))\n", " fixed_dofs=np.asarray(vmap)[(2*boundary[:,None]+np.arange(2)).ravel()]\n", " values=exact_velocity(coordinates[boundary]).ravel()\n", " definition=MultiscaleProblem.from_global(Equation(A,F),binding.size,\n", " fixed=dict(zip(fixed_dofs.tolist(),values.tolist(),strict=True)),constraints=((m,0.7),))\n", " reference=solve(assemble(definition))\n", " velocity=DiscreteField(FieldDefinition(\"velocity\",descriptor=binding.descriptor(component=0)),reference.trace)\n", " pressure=DiscreteField(FieldDefinition(\"pressure\",descriptor=binding.descriptor(component=1)),reference.trace)\n", " points=np.random.default_rng(137).uniform(0,1,size=(100,2))\n", " checks={\"n\":n,\"velocity_max\":float(np.max(abs(velocity.evaluate(points)-exact_velocity(points)))),\n", " \"pressure_max\":float(np.max(abs(pressure.evaluate(points)-exact_pressure(points)))),\n", " \"pressure_integral\":float(m@reference.trace)}\n", " assert max(checks[\"velocity_max\"],checks[\"pressure_max\"])<1e-9\n", " np.testing.assert_allclose(checks[\"pressure_integral\"],0.7,rtol=1e-10,atol=1e-12)\n", " return velocity,pressure,checks\n", " finally:\n", " binding.close()\n", "\n", "reference_rows=[]\n", "for n in (8,16):\n", " reference_velocity,reference_pressure,row=classical_oseen(n)\n", " reference_rows.append(row)\n", "print(\"Independent conforming Taylor–Hood controls:\",reference_rows)\n" ] }, { "cell_type": "markdown", "id": "543e4107", "metadata": {}, "source": [ "## 5. Plot velocity and pressure\n", "\n", "Every panel displays the actual macro mesh. Macrocells are plotted independently,\n", "retaining their own incident field values. The analytical solution is the exact\n", "accuracy reference for this polynomial control.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a424812f", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from matplotlib.collections import LineCollection\n", "\n", "def scalar_panel(ax, fields, exact, title, *, error=False):\n", " \"\"\"Plot each macrocell independently, with original macro topology visible.\"\"\"\n", " artists=[]\n", " for field in fields:\n", " fine=field.mesh\n", " values=field.evaluate(fine.points)\n", " if values.ndim==2: values=np.linalg.norm(values,axis=1)\n", " if exact is not None:\n", " truth=exact(fine.points)\n", " if truth.ndim==2: truth=np.linalg.norm(truth,axis=1)\n", " values=abs(values-truth) if error else truth\n", " artist=ax.tripcolor(fine.points[:,0],fine.points[:,1],fine.cells,values,shading=\"gouraud\")\n", " artists.append(artist)\n", " lower=min(a.get_clim()[0] for a in artists);upper=max(a.get_clim()[1] for a in artists)\n", " for a in artists:a.set_clim(lower,upper)\n", " ax.add_collection(LineCollection(macro.points[macro.faces],colors=\"black\",linewidths=0.65,alpha=0.55))\n", " ax.set(title=title,xlabel=\"x\",ylabel=\"y\",aspect=\"equal\")\n", " return artists[-1]\n", "figure,axes=plt.subplots(2,4,figsize=(15,7),layout=\"constrained\")\n", "for row,fields,exact,name in [(0,velocity,exact_velocity,\"Velocity magnitude\"),(1,pressure,exact_pressure,\"Pressure\")]:\n", " for column,title,truth,error in [(0,\"Analytical\",exact,False),(1,\"MHM\",None,False),(2,\"Classical Taylor–Hood\",None,False),(3,\"Absolute error\",exact,True)]:\n", " chosen=(reference_velocity,) if row==0 and column==2 else (reference_pressure,) if row==1 and column==2 else fields\n", " artist=scalar_panel(axes[row,column],chosen,truth,f\"{title}\\n{name.lower()}\",error=error)\n", " figure.colorbar(artist,ax=axes[row,column],pad=0.025,shrink=0.85)\n", "figure.savefig(OUTPUT/\"oseen-affine-fields.png\",dpi=180)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Define the nonconstant smooth convergence problem\n", "\n", "Now write the same smooth problem used in the refinement curves. The stream\n", "function guarantees divergence-free velocity, and pressure has integral zero:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\psi&=-128x^2(1-x)^2y^2(1-y)^2,\\\\\n", "u_*&=(\\partial_y\\psi,-\\partial_x\\psi),\n", "&p_*&=(x-y)^6-1/28,\\\\\n", "f&=-\\Delta u_*+(\\nabla u_*)\\beta+u_*+\\nabla p_*,\n", "&\\beta&=(1,1)/\\sqrt2.\n", "\\end{aligned}\n", "$$\n", "\n", "Velocity is zero on the exterior. The local equal-order P3/P3 fields require\n", "the consistent Oseen stabilization; Taylor–Hood is a different stable local\n", "choice. Write both the trial and adjoint-test residual before assembling:\n", "\n", "$$\n", "\\begin{aligned}\n", "R(u,p)&=-\\Delta u+(\\nabla u)\\beta+u+\\nabla p,\\\\\n", "R^*(v,q)&=-\\Delta v-(\\nabla v)\\beta+v+\\nabla q,\\\\\n", "a_K^{\\rm US}&=a_K+\\kappa_K(\\nabla\\cdot u,\\nabla\\cdot v)_K\n", " -\\delta_K(R(u,p),R^*(v,q))_K,\\\\\n", "L_K^{\\rm US}&=(f,v)_K-\\delta_K(f,R^*(v,q))_K.\n", "\\end{aligned}\n", "$$\n", "\n", "The minus sign and stabilized source are part of the formulation. For the\n", "computed polynomial inverse constant $m_K$, set $d_K=4/m_K$, $b_K=h_K$,\n", "$\\delta_K=h_K^2/[\\max(h_K^2,d_K)+\\max(d_K,b_K)]$ and\n", "$\\kappa_K=b_K\\min(1,b_K/d_K)$. The generic inverse-bound routine uses Basix\n", "polynomial derivatives; it does not choose the physical equation for us.\n" ] }, { "cell_type": "code", "metadata": { "tags": [ "pymhm-oseen-smooth-definition" ] }, "source": [ "from dolfinx import default_scalar_type\n", "from pymhm.fem.inequalities import laplacian_inverse_bound\n", "from pymhm.fem.scalar.operators import triangle_quadrature\n", "from pymhm.fem.scalar.triangle import tabulate\n", "from examples.solve_oseen import crisscross\n", "\n", "\n", "def smooth_velocity(points):\n", " \"\"\"Differentiate the analytical stream function independently of assembly.\"\"\"\n", " x,y=points.T\n", " sx,sy=x*x*(1-x)**2,y*y*(1-y)**2\n", " dsx,dsy=2*x*(1-x)*(1-2*x),2*y*(1-y)*(1-2*y)\n", " return np.column_stack((-128*sx*dsy,128*dsx*sy))\n", "\n", "\n", "def smooth_pressure(points):\n", " \"\"\"Return the independently stated pressure with zero physical mean.\"\"\"\n", " return (points[:,0]-points[:,1])**6-1/28\n", "\n", "\n", "def smooth_forms(W,domain,*,delta=0.0,kappa=0.0):\n", " \"\"\"Declare the same physical volume equation and optional consistent residuals.\"\"\"\n", " u,p=ufl.TrialFunctions(W)\n", " v,q=ufl.TestFunctions(W)\n", " x=ufl.SpatialCoordinate(domain)\n", " beta=ufl.as_vector((1/np.sqrt(2),1/np.sqrt(2)))\n", " psi=-128*x[0]**2*(1-x[0])**2*x[1]**2*(1-x[1])**2\n", " truth=ufl.as_vector((psi.dx(1),-psi.dx(0)))\n", " truth_p=(x[0]-x[1])**6-1/28\n", " force=-ufl.div(ufl.grad(truth))+ufl.grad(truth)*beta+truth+ufl.grad(truth_p)\n", " dx=ufl.Measure('dx',domain=domain,metadata={'quadrature_degree':20})\n", " physical=(ufl.inner(ufl.grad(u),ufl.grad(v))+ufl.inner(u,v)\n", " +0.5*(ufl.inner(ufl.grad(u)*beta,v)-ufl.inner(u,ufl.grad(v)*beta))\n", " -p*ufl.div(v)-q*ufl.div(u))*dx\n", " residual=-ufl.div(ufl.grad(u))+ufl.grad(u)*beta+u+ufl.grad(p)\n", " adjoint=-ufl.div(ufl.grad(v))-ufl.grad(v)*beta+v+ufl.grad(q)\n", " a=physical+(kappa*ufl.div(u)*ufl.div(v)-delta*ufl.inner(residual,adjoint))*dx\n", " L=(ufl.inner(force,v)-delta*ufl.inner(force,adjoint))*dx\n", " return a,L,q*dx\n", "\n", "\n", "equal_order=basix.ufl.mixed_element([\n", " basix.ufl.element('Lagrange','triangle',3,shape=(2,),\n", " lagrange_variant=basix.LagrangeVariant.equispaced),\n", " basix.ufl.element('Lagrange','triangle',3,\n", " lagrange_variant=basix.LagrangeVariant.equispaced),\n", "])\n", "\n", "\n", "def smooth_local_equations(local):\n", " \"\"\"Write enriched local momentum, incompressibility and signed face pairings.\"\"\"\n", " fine=local.mesh\n", " bary,weights=triangle_quadrature(10)\n", " _,_,_,gradient,hessian=tabulate(fine,3,bary)\n", " h=np.max(fine.lengths[fine.cell_faces],axis=1)\n", " inverse=laplacian_inverse_bound(gradient,hessian,weights,h)\n", " d=4/inverse\n", " delta=float((h*h/(np.maximum(h*h,d)+np.maximum(d,h)))[0])\n", " kappa=float((h*np.minimum(1,h/d))[0])\n", " binding=local.native_space(equal_order)\n", " W=binding.space\n", " u,p=ufl.TrialFunctions(W)\n", " v,q=ufl.TestFunctions(W)\n", " delta_coefficient=fem.Constant(binding.mesh,default_scalar_type(delta))\n", " kappa_coefficient=fem.Constant(binding.mesh,default_scalar_type(kappa))\n", " a,L,pressure_mean=smooth_forms(W,binding.mesh,delta=delta_coefficient,kappa=kappa_coefficient)\n", " dx=ufl.Measure('dx',domain=binding.mesh,metadata={'quadrature_degree':20})\n", " b=local.trace_pairings(lambda phi,ds:ufl.inner(phi,v)*ds)\n", " c=local.trace_pairings(lambda phi,ds:-ufl.inner(phi,u)*ds,axis='rows')\n", " _,velocity_rows=W.sub(0).collapse()\n", " translations=np.zeros((binding.size,2))\n", " translations[velocity_rows]=np.tile(np.eye(2),(len(velocity_rows)//2,1))\n", " local.field('velocity',binding,component=0)\n", " local.field('pressure',binding,component=1)\n", " return local.equations(\n", " a=a,L=L,b=b,c=c,coarse_basis=translations,\n", " moments=columns(v[0]*dx,v[1]*dx),\n", " metadata={'pressure_weights':compile_form(pressure_mean)},\n", " )\n", "\n", "\n", "smooth_macro=crisscross(8)\n", "smooth_trace=SkeletonSpace(\n", " smooth_macro,tuple(FaceSpace.uniform(1) for _ in smooth_macro.faces),components=2,\n", ")\n", "smooth_hierarchy=MeshHierarchy(\n", " smooth_macro,tuple(smooth_macro.submesh(i,1) for i in range(len(smooth_macro.cells))),\n", ")\n", "smooth_problem=bind_problem(\n", " smooth_hierarchy,bind_interface(smooth_trace,convention='normal'),\n", " smooth_local_equations,global_equation=Equation(0,0),retained=2,\n", ")\n", "smooth_system=assemble(smooth_problem)\n", "smooth_weights=[item['pressure_weights'] for item in smooth_system.local_metadata]\n", "smooth_gauge=smooth_system.mean_constraint(smooth_weights,0.0)\n", "smooth_solution=solve(smooth_system,constraints=(smooth_gauge,))\n", "smooth_u,smooth_p=smooth_solution.field('velocity'),smooth_solution.field('pressure')\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Check the physical norms and independently refine Taylor–Hood\n", "\n", "The product norm used by Araya et al. (2021), Section 2.2, is\n", "\n", "$$\n", "\\lVert(e_u,e_p)\\rVert_{V\\times Q}^2\n", "=\\tfrac12\\lVert e_u\\rVert_{L^2}^2\n", " +\\lVert\\nabla e_u\\rVert_{L^2(\\mathcal T_H)}^2\n", " +\\lVert e_p\\rVert_{L^2}^2.\n", "$$\n", "\n", "The smooth P1-trace product-norm order two and additional velocity L2 order\n", "three are numerical observations in Section 5.1 of Araya et al. (2021).\n", "Theorem 7 gives a posteriori reliability and efficiency; it does not prove\n", "these a priori orders. Integrate the physical fields at independent\n", "Gaussian points, keeping the incident local cell identifiers explicit. The\n", "original-row observer checks the complete executed mixed blocks and pressure\n", "gauge separately; an algebraic residual is not this physical product norm.\n", "\n", "For the classical baseline, use the same `smooth_forms` with delta and kappa\n", "zero in the stable global P2/P1 Taylor–Hood space. Strong zero velocity and\n", "one zero pressure integral match the physical problem. Its refinement errors\n", "must decrease before using the finest numerical field as a baseline.\n" ] }, { "cell_type": "code", "metadata": { "tags": [ "pymhm-oseen-smooth-validation" ] }, "source": [ "import json\n", "from types import SimpleNamespace\n", "from examples.core_elasticity_field_archive import ProductionObservation\n", "from examples.minimal_flow_originals import original_diagnostics\n", "from pymhm.fem.vector.curl import quadrature,physical_points\n", "\n", "\n", "def smooth_errors(velocity_fields,pressure_fields,order):\n", " \"\"\"Integrate physical velocity, pressure and broken-gradient errors independently.\"\"\"\n", " sums=np.zeros(3,dtype=np.longdouble)\n", " for vf,pf in zip(velocity_fields,pressure_fields,strict=True):\n", " fine=vf.mesh\n", " bary,weights,_=quadrature(fine,1,order)\n", " points=physical_points(fine,bary)\n", " cells=np.repeat(np.arange(len(fine.cells)),points.shape[1])\n", " flat=points.reshape(-1,2)\n", " numerical_u,numerical_gradient=vf.values_and_gradient(flat,cells=cells)\n", " numerical_p=pf.evaluate(flat,cells=cells)\n", " x,y=flat.T\n", " sx,sy=x*x*(1-x)**2,y*y*(1-y)**2\n", " dx,dy=2*x*(1-x)*(1-2*x),2*y*(1-y)*(1-2*y)\n", " ddx,ddy=2-12*x+12*x*x,2-12*y+12*y*y\n", " exact_gradient=-128*np.array([[dx*dy,sx*ddy],[-ddx*sy,-dx*dy]]).transpose(2,0,1)\n", " eu=numerical_u-smooth_velocity(flat)\n", " ep=numerical_p-smooth_pressure(flat)\n", " eg=numerical_gradient-exact_gradient\n", " w=weights.ravel()\n", " sums+=np.array([np.sum(w*np.sum(eu*eu,axis=1),dtype=np.longdouble),\n", " np.sum(w*ep*ep,dtype=np.longdouble),\n", " np.sum(w*np.sum(eg*eg,axis=(1,2)),dtype=np.longdouble)])\n", " return {'velocity_l2':float(np.sqrt(sums[0])),\n", " 'pressure_l2':float(np.sqrt(sums[1])),\n", " 'mixed_error':float(np.sqrt(sums[0]/2+sums[1]+sums[2]))}\n", "\n", "\n", "observation=ProductionObservation(\n", " system=smooth_system,applied_boundary=np.zeros(smooth_system.trace_size),\n", " mean_weights=[tuple(smooth_weights)],mean_values=[0.0],constraints=[smooth_gauge],\n", ")\n", "smooth_original=original_diagnostics(\n", " SimpleNamespace(hybrid=smooth_solution),observation,\n", " [{'mixed_velocity_pressure':len(values)} for values in smooth_solution.fields],\n", ")\n", "assert smooth_original['accepted']\n", "norms={order:smooth_errors(smooth_u,smooth_p,order) for order in (10,14)}\n", "for name in norms[14]:\n", " np.testing.assert_allclose(norms[10][name],norms[14][name],rtol=1e-8,atol=1e-12)\n", "record_path=ROOT/'examples/results/oseen/smooth-nu1-l1-uniform.json'\n", "record=json.loads(record_path.read_text())\n", "archived=next(row for row in record['rows'] if row['resolution']==8)\n", "for name in norms[14]:\n", " np.testing.assert_allclose(norms[14][name],archived[name],rtol=1e-8,atol=1e-12)\n", "print('Current by-hand UFL and archived smooth norms:',norms[14])\n", "print('Original mixed rows:',{k:v for k,v in smooth_original.items() if k!='blocks'})\n", "\n", "\n", "def classical_smooth_oseen(n):\n", " \"\"\"Assemble the declared physical UFL operator globally in stable Taylor–Hood.\"\"\"\n", " fine=TriangleMesh.unit_square(n)\n", " binding=bind_space(fine,taylor_hood_element)\n", " try:\n", " W=binding.space\n", " a,L,mean=smooth_forms(W,binding.mesh)\n", " A,F,m=compile_form(a),compile_form(L),compile_form(mean)\n", " V,parent=W.sub(0).collapse()\n", " points=V.tabulate_dof_coordinates()[:,:2]\n", " exterior=np.flatnonzero(np.any(np.isclose(points,0)|np.isclose(points,1),axis=1))\n", " fixed=np.asarray(parent)[(2*exterior[:,None]+np.arange(2)).ravel()]\n", " definition=MultiscaleProblem.from_global(\n", " Equation(A,F),binding.size,fixed={int(i):0.0 for i in fixed},constraints=((m,0.0),),\n", " )\n", " state=solve(assemble(definition))\n", " velocity=DiscreteField(FieldDefinition('velocity',descriptor=binding.descriptor(component=0)),state.trace)\n", " pressure=DiscreteField(FieldDefinition('pressure',descriptor=binding.descriptor(component=1)),state.trace)\n", " row={'n':n,**smooth_errors((velocity,),(pressure,),14),\n", " 'global_compatibility_relative':state.raw_residual,'pressure_integral':float(m@state.trace)}\n", " np.testing.assert_allclose(row['pressure_integral'],0.0,atol=1e-10)\n", " return velocity,pressure,row\n", " finally:\n", " binding.close()\n", "\n", "\n", "smooth_classical=[]\n", "for n in (16,32,64):\n", " smooth_reference_u,smooth_reference_p,row=classical_smooth_oseen(n)\n", " smooth_classical.append(row)\n", "assert all(smooth_classical[i+1][name]