From ac76909406c69e75ddd7198587842d67bd0f1ef2 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Fri, 8 Jul 2022 15:43:10 +0200 Subject: [PATCH 1/2] improve export.triplot This patch generalizes export.triplot by allowing a matplotlib axes object, rather than a file name, as the first argument. This makes it possible to stack additional items on top of the base figure without having to reimplement triplot locally. --- nutils/export.py | 77 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 51 insertions(+), 26 deletions(-) diff --git a/nutils/export.py b/nutils/export.py index c1494e958..d6d1b795e 100644 --- a/nutils/export.py +++ b/nutils/export.py @@ -47,34 +47,59 @@ def plotlines_(ax, xy, lines, **kwargs): return lc -def triplot(name, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None): - if (tri is None) != (values is None): +def _triplot_1d(ax, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None): + if plabel: + ax.set_xlabel(plabel) + if vlabel: + ax.set_ylabel(vlabel) + if hull is not None: + for x in points[hull[:, 0], 0]: + ax.axvline(x, color=linecolor, linewidth=linewidth) + if tri is not None: + plotlines_(ax, [points[:, 0], values], tri) + ax.autoscale(enable=True, axis='x', tight=True) + if clim is None: + ax.autoscale(enable=True, axis='y', tight=False) + else: + ax.set_ylim(clim) + + +def _triplot_2d(ax, points, values=None, *, tri=None, hull=None, cmap=None, clim=None, linewidth=.1, linecolor='k', plabel=None, vlabel=None): + ax.set_aspect('equal') + if plabel: + ax.set_xlabel(plabel) + ax.set_ylabel(plabel) + if tri is not None: + im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True) + if clim is not None: + im.set_clim(clim) + else: + im = None + if hull is not None: + plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5) + ax.autoscale(enable=True, axis='both', tight=True) + return im + + +def triplot(name, points, values=None, **kwargs): + + if points.shape[1] == 1: + _triplot = _triplot_1d + elif points.shape[1] == 2: + _triplot = _triplot_2d + else: + raise Exception(f'invalid spatial dimension: {nd}') + + if (kwargs.get('tri') is None) != (values is None): raise Exception('tri and values can only be specified jointly') + + if not isinstance(name, str): + return _triplot(name, points, values, **kwargs) + with mplfigure(name) as fig: - if points.shape[1] == 1: - ax = fig.add_subplot(111, xlabel=plabel, ylabel=vlabel) - if tri is not None: - plotlines_(ax, [points[:, 0], values], tri) - if hull is not None: - for x in points[hull[:, 0], 0]: - ax.axvline(x, color=linecolor, linewidth=linewidth) - ax.autoscale(enable=True, axis='x', tight=True) - if clim is None: - ax.autoscale(enable=True, axis='y', tight=False) - else: - ax.set_ylim(clim) - elif points.shape[1] == 2: - ax = fig.add_subplot(111, xlabel=plabel, ylabel=plabel, aspect='equal') - if tri is not None: - im = ax.tripcolor(*points.T, tri, values, shading='gouraud', cmap=cmap, rasterized=True) - if clim is not None: - im.set_clim(clim) - fig.colorbar(im, label=vlabel) - if hull is not None: - plotlines_(ax, points.T, hull, colors=linecolor, linewidths=linewidth, alpha=1 if tri is None else .5) - ax.autoscale(enable=True, axis='both', tight=True) - else: - raise Exception('invalid spatial dimension: {}'.format(points.shape[1])) + im = _triplot(fig.add_subplot(111), points, values, **kwargs) + if im: + fig.colorbar(im, label=kwargs.get('vlabel')) @util.positional_only From e1895d7fb04e4e17474694b24bfa875f867a0d76 Mon Sep 17 00:00:00 2001 From: Gertjan van Zwieten Date: Mon, 1 Aug 2022 14:31:43 +0200 Subject: [PATCH 2/2] unify drivencavity examples This patch merges drivencavity-compatible into drivencavity. The new example script contains boolean flags for compatible spaces and strong imposition of tangential velocity. It also improves the vizualisation and adds two cross-section figures for comparison with benchmark results. --- examples/drivencavity-compatible.py | 157 ----------------- examples/drivencavity.py | 257 ++++++++++++++++++++-------- 2 files changed, 190 insertions(+), 224 deletions(-) delete mode 100644 examples/drivencavity-compatible.py diff --git a/examples/drivencavity-compatible.py b/examples/drivencavity-compatible.py deleted file mode 100644 index 1a2e9c940..000000000 --- a/examples/drivencavity-compatible.py +++ /dev/null @@ -1,157 +0,0 @@ -#! /usr/bin/env python3 -# -# In this script we solve the lid driven cavity problem for stationary Stokes -# and Navier-Stokes flow. That is, a unit square domain, with no-slip left, -# bottom and right boundaries and a top boundary that is moving at unit -# velocity in positive x-direction. -# -# The script is identical to :ref:`examples/drivencavity.py` except that it -# uses the Raviart-Thomas discretization providing compatible velocity and -# pressure spaces resulting in a pointwise divergence-free velocity field. - -from nutils import mesh, function, solver, export, cli, testing -from nutils.expression_v2 import Namespace -import numpy -import treelog - - -def main(nelems: int, degree: int, reynolds: float): - ''' - Driven cavity benchmark problem using compatible spaces. - - .. arguments:: - - nelems [12] - Number of elements along edge. - degree [2] - Polynomial degree for velocity; the pressure space is one degree less. - reynolds [1000] - Reynolds number, taking the domain size as characteristic length. - ''' - - verts = numpy.linspace(0, 1, nelems+1) - domain, geom = mesh.rectilinear([verts, verts]) - - ns = Namespace() - ns.δ = function.eye(domain.ndims) - ns.Σ = function.ones([domain.ndims]) - ns.x = geom - ns.Re = reynolds - - ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('u', 'v'), function.vectorize([ - domain.basis('spline', degree=(degree, degree-1), removedofs=((0, -1), None)), - domain.basis('spline', degree=(degree-1, degree), removedofs=(None, (0, -1)))])) - ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1)) - ns.add_field(('λ', 'μ')) - - ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij' - ns.uwall = numpy.stack([domain.boundary.indicator('top'), 0]) - ns.N = 5 * degree * nelems # nitsche constant based on element size = 1/nelems - ns.nitsche_i = '(N v_i - (∇_j(v_i) + ∇_i(v_j)) n_j) / Re' - - res = domain.integral('∇_j(v_i) σ_ij dV' @ ns, degree=2*degree) - res += domain.boundary.integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dS' @ ns, degree=2*degree) - res += domain.integral('q (∇_k(u_k) + λ) dV' @ ns, degree=2*degree) - res += domain.integral('μ p dV' @ ns, degree=2*degree) - - with treelog.context('stokes'): - args0 = solver.solve_linear('u:v,p:q,λ:μ', res) - postprocess(domain, ns, **args0) - - res += domain.integral('v_i ∇_j(u_i) u_j dV' @ ns, degree=3*degree) - with treelog.context('navierstokes'): - args1 = solver.newton('u:v,p:q,λ:μ', res, arguments=args0).solve(tol=1e-10) - postprocess(domain, ns, **args1) - - return args0, args1 - -# Postprocessing in this script is separated so that it can be reused for the -# results of Stokes and Navier-Stokes, and because of the extra steps required -# for establishing streamlines. - - -def postprocess(domain, ns, every=.05, spacing=.01, **arguments): - - div = domain.integral('∇_k(u_k)^2 dV' @ ns, degree=1).eval(**arguments)**.5 - treelog.info('velocity divergence: {:.2e}'.format(div)) # confirm that velocity is pointwise divergence-free - - ns = ns.copy_() # copy namespace so that we don't modify the calling argument - ns.add_field('stream', domain.basis('std', degree=2)[1:]) # stream function (sans first dof to obtain non-singular system) - ns.ε = function.levicivita(2) - sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(stream))^2 dV' @ ns, degree=4) - arguments = solver.optimize('stream,', sqr, arguments=arguments) # compute streamlines - - bezier = domain.sample('bezier', 9) - x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'stream'] @ ns, **arguments) - with export.mplfigure('flow.png') as fig: # plot velocity as field, pressure as contours, streamlines as dashed - ax = fig.add_axes([.1, .1, .8, .8], yticks=[], aspect='equal') - import matplotlib.collections - ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2)) - ax.tricontour(x[:, 0], x[:, 1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9) - caxu = fig.add_axes([.1, .1, .03, .8], title='velocity') - imu = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, u, shading='gouraud', cmap='jet') - fig.colorbar(imu, cax=caxu) - caxu.yaxis.set_ticks_position('left') - caxp = fig.add_axes([.87, .1, .03, .8], title='pressure') - imp = ax.tricontour(x[:, 0], x[:, 1], bezier.tri, p, 16, cmap='gray', linestyles='solid') - fig.colorbar(imp, cax=caxp) - -# If the script is executed (as opposed to imported), :func:`nutils.cli.run` -# calls the main function with arguments provided from the command line. To -# keep with the default arguments simply run :sh:`python3 -# drivencavity-compatible.py`. - - -if __name__ == '__main__': - cli.run(main) - -# Once a simulation is developed and tested, it is good practice to save a few -# strategic return values for regression testing. The :mod:`nutils.testing` -# module, which builds on the standard :mod:`unittest` framework, facilitates -# this by providing :func:`nutils.testing.TestCase.assertAlmostEqual64` for the -# embedding of desired results as compressed base64 data. - - -class test(testing.TestCase): - - def test_p2(self): - args0, args1 = main(nelems=3, reynolds=100, degree=2) - with self.subTest('stokes-velocity'): - self.assertAlmostEqual64(args0['u'], ''' - eNrzu9Bt8OuUndkD/eTTSqezzP2g/E3698/ZmZlf2GjSaHJS3/90/Wm/C4qGh04CAErOF+g=''') - with self.subTest('stokes-pressure'): - self.assertAlmostEqual64(args0['p'], ''' - eNoBIADf/0fSMC4P0ZLKjCk4JC7Pwy901sjb0jA90Lkt0NHxLm41+KgP+Q==''') - with self.subTest('stokes-multiplier'): - self.assertAlmostEqual(args0['λ'], 0) - with self.subTest('navier-stokes-velocity'): - self.assertAlmostEqual64(args1['u'], ''' - eNoBMADP/2XOWjJSy5k1jS+yyzvLODfgL1rO0MrINpsxHM2ZNSrPqDTANCPVQsxCzeAvcc04yT3AF9c=''') - with self.subTest('navier-stokes-pressure'): - self.assertAlmostEqual64(args1['p'], ''' - eNpbqpasrWa46TSTfoT+krO/de/pntA3Pnf2zFNtn2u2hglmAOKVDlE=''') - with self.subTest('navier-stokes-multiplier'): - self.assertAlmostEqual(args1['λ'], 0) - - def test_p3(self): - args0, args1 = main(nelems=3, reynolds=100, degree=3) - with self.subTest('stokes-velocity'): - self.assertAlmostEqual64(args0['u'], ''' - eNqbo3f7/AeDb2dmGGtd7r+odk7icoQJgjUHLpty8b/BvrMzjF/rll9qMD5kwAAFopc6dRvO2J2fo8d4 - 3sko4wwAjL4lyw==''', atol=1e-14) - with self.subTest('stokes-pressure'): - self.assertAlmostEqual64(args0['p'], ''' - eNpbrN+us+Z8qWHMyfcXrly013l1ZocRAxwI6uvoHbwsZuxxNvZC5eUQg+5zS8wAElAT9w==''') - with self.subTest('stokes-multiplier'): - self.assertAlmostEqual(args0['λ'], 0) - with self.subTest('navier-stokes-velocity'): - self.assertAlmostEqual64(args1['u'], ''' - eNoBUACv/14yGcxyNPbJYTahLj/LSDE7yy43SM9WMsXJoDR+N3Iw8s1hM5zJODeizcE0X8phNrQwUDOO - NbMzJi+ty4s1oDFqzxIzysjWzXIwFM3tNMjIpDMlJw==''') - with self.subTest('navier-stokes-pressure'): - self.assertAlmostEqual64(args1['p'], ''' - eNoBMgDN/yoxvDHizaEy88kDLgoviCt4znIzMy7jL7nTRMzqzhAwBTA/LE0w6czY2GQwoNEYMHE3NDUW - DQ==''') - with self.subTest('navier-stokes-multiplier'): - self.assertAlmostEqual(args1['λ'], 0) diff --git a/examples/drivencavity.py b/examples/drivencavity.py index fc2731b29..ea7ac718a 100644 --- a/examples/drivencavity.py +++ b/examples/drivencavity.py @@ -1,23 +1,73 @@ #! /usr/bin/env python3 # # In this script we solve the lid driven cavity problem for stationary Stokes -# and Navier-Stokes flow. That is, a unit square domain, with no-slip left, -# bottom and right boundaries and a top boundary that is moving at unit -# velocity in positive x-direction. +# and Navier-Stokes flow. This benchmark problem consists of a square domain +# with fixed left, bottom and right boundaries and a top boundary that is +# moving at unit velocity in positive x-direction. Reference results can be +# found for instance at https://www.acenumerics.com/the-benchmarks.html. +# +# The general conservation laws are: +# +# Dρ/Dt = 0 (mass) +# ρ Du_i/Dt = ∇_j(σ_ij) (momentum) +# +# where we used the material derivative D·/Dt := ∂·/∂t + ∇_j(· u_j). The stress +# tensor is σ_ij := μ (∇_j(u_i) + ∇_i(u_j) - 2 ∇_k(u_k) δ_ij / δ_nn) - p δ_ij, +# with pressure p and dynamic viscosity μ, and ρ is the density. +# +# Assuming steady, incompressible flow, we take density to be constant. Further +# introducing a reference length L and reference velocity U, we make the +# equations dimensionless by taking spatial coordinates relative to L, velocity +# relative to U, and pressure relative to ρ U^2. This reduces the conservation +# laws to: +# +# ∇_k(u_k) = 0 (mass) +# Du_i/Dt = ∇_j(σ_ij) (momentum) +# +# where the material derivative simplifies to D·/Dt := ∇_j(·) u_j, and the +# stress tensor becomes σ_ij := (∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij, with +# Reynolds number Re = ρ U L / μ. +# +# The weak form is obtained by multiplication of a test function and partial +# integration of the right hand side of the momentum balance. +# +# ∀ q: ∫_Ω q ∇_k(u_k) = 0 (mass) +# ∀ v: ∫_Ω (v_i Du_i/Dt + ∇_j(v_i) σ_ij) = ∫_Γ v_i σ_ij n_j (momentum) +# +# A remaining issue with this system is that its solution is not unique due to +# its strictly kinematic boundary conditions: since ∫_Ω ∇_k(u_k) = 0 for any u +# that satisfies the non-penetrating boundary conditions, any pressure space +# that contains unity results in linear dependence. Furthermore, the strong +# form shows that constant pressure shifts do not affect the momentum balance. +# Both issues are solved by arbitrarily removing one of the basis functions. +# +# The normal velocity components at the wall are strongly constrained. The +# tangential components are either strongly or weakly constrained depending on +# the `strongbc` parameter. Since the tangential components at the top are +# incompatible with the normal components at the left and right boundary, the +# constraints are constructed in two steps to avoid Gibbs type oscillations, +# and to make sure that the non-penetrating condition takes precedence. +# +# Depending on the `compatible` parameter, the script uses either a Taylor-Hood +# (False) or a Raviart-Thomas (True) discretization. In case of TH, the system +# is consistently modified by adding ∫_Ω .5 u_i v_i ∇_j(u_j) to yield the skew- +# symmetric advective term ∫_Ω .5 (v_i ∇_j(u_i) - u_i ∇_j(v_i)) u_j. In case of +# RT, the discretization guarantees a pointwise divergence-free velocity field, +# and skew symmetry is implied. from nutils import mesh, function, solver, export, cli, testing from nutils.expression_v2 import Namespace +import treelog as log import numpy -import treelog -def main(nelems: int, etype: str, degree: int, reynolds: float): +def main(nelems: int, etype: str, degree: int, reynolds: float, compatible: bool, strongbc: bool): ''' Driven cavity benchmark problem. .. arguments:: - nelems [12] + nelems [32] Number of elements along edge. etype [square] Element type (square/triangle/mixed). @@ -25,36 +75,69 @@ def main(nelems: int, etype: str, degree: int, reynolds: float): Polynomial degree for velocity; the pressure space is one degree less. reynolds [1000] Reynolds number, taking the domain size as characteristic length. + strongbc [no] + Use strong boundary constraints + compatible [no] + Use compatible spaces and weakly imposed boundary conditions; requires + etype=square and strongbc=no ''' + if compatible and (strongbc or etype != 'square'): + raise Exception(f'compatible mode requires square elements and weak boundary conditions') + domain, geom = mesh.unitsquare(nelems, etype) ns = Namespace() ns.δ = function.eye(domain.ndims) + ns.ε = function.levicivita(2) ns.Σ = function.ones([domain.ndims]) ns.Re = reynolds + ns.uwall = numpy.stack([domain.boundary.indicator('top'), 0]) ns.x = geom ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS')) - ns.add_field(('u', 'v'), domain.basis('std', degree=degree), shape=(domain.ndims,)) - ns.add_field(('p', 'q'), domain.basis('std', degree=degree-1)) + if not compatible: + ns.add_field(('u', 'v'), domain.basis('std', degree=degree), shape=(domain.ndims,)) + ns.add_field(('p', 'q'), domain.basis('std', degree=degree-1)[1:]) + ns.add_field('ψ', domain.basis('std', degree=2)[1:]) + else: + ns.add_field(('u', 'v'), function.vectorize([ + domain.basis('spline', degree=(degree, degree-1)), + domain.basis('spline', degree=(degree-1, degree))])) + ns.add_field(('p', 'q'), domain.basis('spline', degree=degree-1)[1:]) + ns.add_field('ψ', domain.basis('spline', degree=degree)[1:]) ns.σ_ij = '(∇_j(u_i) + ∇_i(u_j)) / Re - p δ_ij' - sqr = domain.boundary.integral('u_k u_k dS' @ ns, degree=degree*2) - sqr += domain.locate(ns.x, [(1,1)], weights=numpy.array([1.]), tol=1e-10).integral(ns.p**2) - cons = solver.optimize('u,p', sqr, droptol=1e-15) + # weak formulation for Stokes flow, over-integrating for improved + # efficiency when we turn this to Navier-Stokes later on + res = domain.integral('∇_j(v_i) σ_ij dV' @ ns, degree=degree*3) + res += domain.integral('q ∇_k(u_k) dV' @ ns, degree=degree*3) - # project lid velocity separately to avoid Gibbs wiggles around boundary incompatibilities - sqr = domain.boundary['top'].integral('(u_0 - 1)^2 dS' @ ns, degree=degree*2) - lidcons = solver.optimize('u,', sqr, droptol=1e-15) - cons['u'] = numpy.choose(numpy.isnan(lidcons['u']), [lidcons['u'], cons['u']]) + # strong enforcement of non-penetrating boundary conditions + sqr = domain.boundary.integral('(u_k n_k)^2 dS' @ ns, degree=degree*2) + cons = solver.optimize('u,', sqr, droptol=1e-15) - res = domain.integral('(∇_j(v_i) σ_ij + q ∇_k(u_k)) dV' @ ns, degree=degree*2) - with treelog.context('stokes'): + if strongbc: + # strong enforcement of tangential boundary conditions + sqr = domain.boundary.integral('(ε_ij n_i (u_j - uwall_j))^2 dS' @ ns, degree=degree*2) + tcons = solver.optimize('u,', sqr, droptol=1e-15) + cons['u'] = numpy.choose(numpy.isnan(cons['u']), [cons['u'], tcons['u']]) + else: + # weak enforcement of tangential boundary conditions via Nitsche's method + ns.N = 5 * degree * nelems # Nitsche constant based on element size = 1/nelems + ns.nitsche_i = '(N v_i - (∇_j(v_i) + ∇_i(v_j)) n_j) / Re' + res += domain.boundary.integral('(nitsche_i (u_i - uwall_i) - v_i σ_ij n_j) dS' @ ns, degree=2*degree) + + with log.context('stokes'): args0 = solver.solve_linear('u:v,p:q', res, constrain=cons) postprocess(domain, ns, **args0) - res += domain.integral('.5 (v_i ∇_j(u_i) - u_i ∇_j(v_i)) u_j dV' @ ns, degree=degree*3) - with treelog.context('navierstokes'): + # change to Navier-Stokes by adding convection + res += domain.integral('v_i ∇_j(u_i) u_j dV' @ ns, degree=degree*3) + if not compatible: + # add consistent term for skew-symmetry + res += domain.integral('.5 u_i v_i ∇_j(u_j) dV' @ ns, degree=degree*3) + + with log.context('navier-stokes'): args1 = solver.newton('u:v,p:q', res, arguments=args0, constrain=cons).solve(tol=1e-10) postprocess(domain, ns, **args1) @@ -64,35 +147,43 @@ def main(nelems: int, etype: str, degree: int, reynolds: float): # results of Stokes and Navier-Stokes, and because of the extra steps required # for establishing streamlines. +def postprocess(domain, ns, **arguments): -def postprocess(domain, ns, every=.05, spacing=.01, **arguments): - - ns = ns.copy_() # copy namespace so that we don't modify the calling argument - ns.add_field('stream', domain.basis('std', degree=2)[1:]) # stream function (sans first dof to obtain non-singular system) - ns.ε = function.levicivita(2) - sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(stream))^2 dV' @ ns, degree=4) - arguments = solver.optimize('stream,', sqr, arguments=arguments) # compute streamlines + # reconstruct velocity streamlines + sqr = domain.integral('Σ_i (u_i - ε_ij ∇_j(ψ))^2 dV' @ ns, degree=4) + arguments = solver.optimize('ψ,', sqr, arguments=arguments) bezier = domain.sample('bezier', 9) - x, u, p, stream = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'stream'] @ ns, **arguments) - with export.mplfigure('flow.png') as fig: # plot velocity as field, pressure as contours, streamlines as dashed - ax = fig.add_axes([.1, .1, .8, .8], yticks=[], aspect='equal') - import matplotlib.collections - ax.add_collection(matplotlib.collections.LineCollection(x[bezier.hull], colors='w', linewidths=.5, alpha=.2)) - ax.tricontour(x[:, 0], x[:, 1], bezier.tri, stream, 16, colors='k', linestyles='dotted', linewidths=.5, zorder=9) - caxu = fig.add_axes([.1, .1, .03, .8], title='velocity') - imu = ax.tripcolor(x[:, 0], x[:, 1], bezier.tri, u, shading='gouraud', cmap='jet') - fig.colorbar(imu, cax=caxu) - caxu.yaxis.set_ticks_position('left') - caxp = fig.add_axes([.87, .1, .03, .8], title='pressure') - imp = ax.tricontour(x[:, 0], x[:, 1], bezier.tri, p, 16, cmap='gray', linestyles='solid') - fig.colorbar(imp, cax=caxp) + x, u, p, ψ = bezier.eval(['x_i', 'sqrt(u_i u_i)', 'p', 'ψ'] @ ns, **arguments) + + with export.mplfigure('flow.png', dpi=150) as fig: # plot velocity as field, pressure as contours, streamlines as dashed + ax = fig.add_subplot(111, aspect='equal') + im = export.triplot(ax, x, u, tri=bezier.tri, hull=bezier.hull, cmap='hot_r', clim=(0,1)) + fig.colorbar(im, label='velocity') + ax.tricontour(x[:, 0], x[:, 1], bezier.tri, ψ, levels=numpy.percentile(ψ, numpy.arange(2,100,3)), colors='k', linestyles='solid', linewidths=.5, zorder=9) + + x = numpy.linspace(0, 1, 1001) + v = domain.locate(ns.x, numpy.stack([x, numpy.repeat(.5, len(x))], 1), tol=1e-10).eval(ns.u[1], **arguments) + imin = numpy.argmin(v) + imax = numpy.argmax(v) + with export.mplfigure('cross-hor.png', dpi=150) as fig: + ax = fig.add_subplot(111, xlim=(0,1), title='horizontal cross section at y=0.5', xlabel='x-coordinate', ylabel='vertical velocity') + ax.plot(x, v) + ax.annotate(f'x={x[imax]:.3f}\nv={v[imax]:.3f}', xy=(x[imax], v[imax]), xytext=(x[imax], 0), arrowprops=dict(arrowstyle="->"), ha='center', va='center') + ax.annotate(f'x={x[imin]:.3f}\nv={v[imin]:.3f}', xy=(x[imin], v[imin]), xytext=(x[imin], 0), arrowprops=dict(arrowstyle="->"), ha='center', va='center') + + y = numpy.linspace(0, 1, 1001) + u = domain.locate(ns.x, numpy.stack([numpy.repeat(.5, len(y)), y], 1), tol=1e-10).eval(ns.u[0], **arguments) + imin = numpy.argmin(u) + with export.mplfigure('cross-ver.png', dpi=150) as fig: + ax = fig.add_subplot(111, ylim=(0,1), title='vertical cross section at x=0.5', ylabel='y-coordinate', xlabel='horizontal velocity') + ax.plot(u, y) + ax.annotate(f'y={y[imin]:.3f}\nu={u[imin]:.3f}', xy=(u[imin], y[imin]), xytext=(0, y[imin]), arrowprops=dict(arrowstyle="->"), ha='left', va='center') # If the script is executed (as opposed to imported), :func:`nutils.cli.run` # calls the main function with arguments provided from the command line. To # keep with the default arguments simply run :sh:`python3 drivencavity.py`. - if __name__ == '__main__': cli.run(main) @@ -102,50 +193,82 @@ def postprocess(domain, ns, every=.05, spacing=.01, **arguments): # this by providing :func:`nutils.testing.TestCase.assertAlmostEqual64` for the # embedding of desired results as compressed base64 data. - class test(testing.TestCase): - def test_square(self): - args0, args1 = main(nelems=3, etype='square', reynolds=100, degree=3) + def test_baseline(self): + args0, args1 = main(nelems=3, etype='square', degree=2, reynolds=100, compatible=False, strongbc=False) with self.subTest('stokes-velocity'): self.assertAlmostEqual64(args0['u'], ''' - eNpjYCAMgswhdJ1O+uWtl7/rp13hMyq4rmx8/cI/44qLW0wMLliZMhrxGMHUPT//WmfruW8GWudSjJ6d - FTeWP/vf+NH5TuNVhurGPqZvL8LUbTi3X+P1WV2D2rPthjZnw4yenflpdODSFaO9RpuNZple14Wp+362 - VzP87Ce9b2f0DHaduW+w7EyEIYPeH4MW4z6Dh6apSOqKr4Wf5bv47cyl87vOKJ5fdmbFOQY9lvMtxkXn - H5rOvoSw1/H667OXz9eerTxnc3bV2Wdn2M8euKRzdq+R79lZppqXEP4Qvbz1HNd5rXNzzj47+/KM/FnG - M4/Ol59ZZXjzjI+psB4iXGbqbL3MeSHtyqezBdfvnrl+gelMxUWf0wYXjp1iNPpyFqaOmHAGAFFhkvE=''') + eNo9jSsLwlAcxQ82gwOjoFZ1ONA9NK6vC2Ky2cU2QeziFxDBJIKfZLv/O92DYRqo6waDGBR3GXjKgfPg + Bwh9A6AYA80TMOwCU0OkcwVonGMl8Uaa5C90i9+MipP2tn1gEgFXvgpDGqgRK+smPTWbaqqcdWMfeBPw + YcDGBdLMJR3Y/X+zdkhHHrEHM6lENt25+OU8OUi8PUn+klm87lacqiN4uQrZ4tUCLh3g4AFrV6Q/uctG + gQ==''') with self.subTest('stokes-pressure'): self.assertAlmostEqual64(args0['p'], ''' - eNp7dOb9mRdnHp2pPbPw9MzTN848OXPpzJ4z1Wcizqw/c//Ma6DM9TMHzsw/s+PMFxTIdfbvGfazwmf1 - zkaftTgrdJblrORZ47NTz94463qW/ezPM4xAcsLZjZcZGAAX7kL6''') + eNoBHgDh/+vVsNEXy6jTbiq1z7Av9C0mLJkw1NDTLEEtEC/xNAwED0s=''') with self.subTest('navier-stokes-velocity'): self.assertAlmostEqual64(args1['u'], ''' - eNpjYCAMgswh9ErtA5c9ruTp/7xiZbhX28jo8MVbRn1XLIxVLhcbBxuIGsDUHbhgoxNx3tnA9vwcQ4fz - PkbC59mNP2rONOIxnGiUaOx9GaYu5PxOjfJzB/Xdz5kbWp9jNYo9t81ont4so1OGXUbNJtX6MHVd515p - XT4rqT//7EWDprPzDR+eTTY6pBdltNhoq9ELE+3zMHWe58rVl53lv2R1Vvbq+zOTdCLPmhpONkgwlDMO - NawyvWAIU/f7nMRN03PyF3rOSp6XOqt3bv+ZA+cirnSclzf2vxhvGgo3T/pC5xXW80Ln751ddzb1rOpZ - wzOXTp89l3iu1vic0WrTImOYugCd8uvrriy/aHf1w9nkizPPvDl/7rTe+cRTBmf3nVQx0T4DU0dMOAMA - p1CDeg==''') + eNpjYAABz4sMDP81GRh69BgYTgBxkQlIdCVQ5Nb5R/oSV+KNW8/+NXlidMOE4WS+mZM5A8M1JQaGuLMt + ejPOpBgtOs1vonRe0IT59ErjfWZAuasMDB1n9M9vOG2kL3B6l1H0uUXGRQZXjVYB9a2/xMCQfW7FueVn + nM/5nf5xevqZxDOq5w4bCwLlOoD6XDV/n1t//s5ZvjPzTjmdDjx55+Slky/MGaDgHFB3vz4DgynQfS9O + A3WcBIkCAB7aSkk=''') with self.subTest('navier-stokes-pressure'): self.assertAlmostEqual64(args1['p'], ''' - eNoNiT0OQEAYBdd5xD1cSOIMiIJGVJuIiliS9VNYKolkNG6j9BUvmZlnmJnoSAnx6RmFNSUVjgErRVOQ - iFtWKTUZMQ2n/BuGnJaFi52bl48D73Fis+wjCpT6AWgsRHE=''') + eNpz01W9oHVmuU7SJYtzgherdcr0n59dfiZT11yP97yCGQDN0Azu''') def test_mixed(self): - args0, args1 = main(nelems=3, etype='mixed', reynolds=100, degree=2) + args0, args1 = main(nelems=3, etype='mixed', degree=2, reynolds=100, compatible=False, strongbc=False) + with self.subTest('stokes-velocity'): + self.assertAlmostEqual64(args0['u'], ''' + eNpjYAABHx0Ghrbz+lcYGJZpR2hrnDtm/EObgWHvWSGD1WeuGTIwmJ9jYLAwnH32p8nKMzFmFqfVTJTP + aBszMOw0AenebmJh9tuMgWHGWX9jQ3MGhr9nLA35zxRcWHOm4mzOBQaG55cZGCTPGV6fdUrhwtEzvhe2 + n+Y8k3RG+Mwio99nuoHqg4G48WzCmTignYUXDfXNzoedATuL4bMeA0Op9qczWqfXnTl2ioHhINAdHufv + ntx18qoZSH7FSRAJAB13Sc0=''') + with self.subTest('stokes-pressure'): + self.assertAlmostEqual64(args0['p'], ''' + eNp7pKl+nf1KznmxS62ns/W+az/TNTL4ondU1/46t6GKKQDiJg1H''') + with self.subTest('navier-stokes-velocity'): + self.assertAlmostEqual64(args1['u'], ''' + eNpjYACBVVcYGH5fZAKSChcLLv05p2jsp8XA4Hhe4YrV2QmGDAxHzzMwuBpknHcz4T3nYGp7JslY+ewy + YwaGamOQ7jyjzaZJZgwMYaeFTR4D6TfnvhgUnanXTzuz4WzmBQaG6MsMDEcueOpxnF5iwHZ+kfHaUypn + n5xefpbrzEYjZ3MGhiogDr/YYbxbjYHhrH6lYcY55zNgZzGcNWBgUL0Uctr3zLzTt08xMOScZmCYdnbl + qQMnpcxB8konQSQACVZG3A==''') + with self.subTest('navier-stokes-pressure'): + self.assertAlmostEqual64(args1['p'], ''' + eNrbqjVZs1/ry/n48z1nSrW9L83RkTmneNZMO/TCOUNbMwDktQ3z''') + + def test_compatible(self): + args0, args1 = main(nelems=3, etype='square', degree=2, reynolds=100, compatible=True, strongbc=False) + with self.subTest('stokes-velocity'): + self.assertAlmostEqual64(args0['u'], ''' + eNpjYIAAvwvdBr9O2Zk90E8+rXQ6yxzGZ4CDTfr3z0H45hc2mjSagFgn9f1P15+G6Fc0PHSSgQEAx7kX + 6A==''') + with self.subTest('stokes-pressure'): + self.assertAlmostEqual64(args0['p'], ''' + eNoL1u+7NOfUR929ugvORxlU6W7V1TcUuyiif/PKCf1yUwDfRw2t''') + with self.subTest('navier-stokes-velocity'): + self.assertAlmostEqual64(args1['u'], ''' + eNpjYICA1HNRRkGnZ5r26m86bX3awvyBftS5C6dOmDHAwWxDmbMzTUEsrfMrTA6YgFjKV53OOJ0FsR7o + F561OMnAAAC5tRfX''') + with self.subTest('navier-stokes-pressure'): + self.assertAlmostEqual64(args1['p'], ''' + eNoz1VYx3HT6t16w/uKz73Uv6R7RNzx35swh7XdXrQ0TzADbMQ6l''') + + def test_strong(self): + args0, args1 = main(nelems=3, etype='square', degree=2, reynolds=100, compatible=False, strongbc=True) with self.subTest('stokes-velocity'): self.assertAlmostEqual64(args0['u'], ''' - eNpjYEAFEy++uHzs/EYjEFv73Hm9T2eVDGFyMWdfGJ/TVTGxMvQyqjxbAlYTZM7AsNWIxwhEG5hs0wTR - 1Wd5DDTOHrx05OzmczC9Cud8riSccbrqYJR2dfMZ07M+hkznLpuongepB+Eyk/TzIHVbL4QbrLqw9Qyy - m+aee31awWALXEzP+NIZkB6Y/QAD2Dbr''') + eNpjYMAPDl2wNEg9p2D8QeOQcafBJ9OTJ6abB5lD5E6eVb348oyVkfSZf8YFZ6RNZp6+ZwiTO3KGgUEP + iNedYmDoPc3AsMGEgQGh77beyzPHzkqfYTpTcObp6Zmnlc7B5A5dOH4+9dyDMx807M50GvCdOnki8wRM + DhcAAEYiNtQ=''') with self.subTest('stokes-pressure'): self.assertAlmostEqual64(args0['p'], ''' - eNoBIADf/1zNa81jzWHNs8w3zV3MUs2HzZrNzc26zanNd82qzgAAR/MTmQ==''') + eNoBHgDh/3fRlNYxy7PR0NVKz1ktVi1E1HowsdGJ07Qt/9PINA6QEUk=''') with self.subTest('navier-stokes-velocity'): self.assertAlmostEqual64(args1['u'], ''' - eNpjYEAFvRcfXLa8eNsQxP5xLkBv17ktBjC5iHOfjbr1XxotNbAy0j4nbQQSCzJnYLA0lNIH0XLGi3VB - dNK5HoNjZ2frXj+77BxM775zfbq6Z+cYxBpd1vc543jursGSC0omz86D1IPwBJOzYDszL5oaLbgQcQbZ - TTZnec9svix4FsZXNbl9GqQHZj8AAcY1/g==''') + eNpjYMAPAs4H6M8zaDJOO91vKma628TihKx5kDlEbv2Z5fqFZ24aKZ/mNplmmGJy5eRbY5jcpdOGF+tP + /9BPPW1gXH6W2eSpkds5mBz72fvnUs/EnHt+etXpCWdZzgSd3W8Ck1M9L3zGGyi/4Pz80+ZnNp44c9L8 + BEwOFwAA4RM3QA==''') with self.subTest('navier-stokes-pressure'): self.assertAlmostEqual64(args1['p'], ''' - eNoBIADf/wXMDswMzAfMzMvry73LAcwIzCDM/ssJzDTM9MvRzAAAHqQRsw==''') + eNoBHgDh//ot489SzMEsntHezWTPAC+jL+XN/8wF02UxTc1JNhf0ELc=''')