diff --git a/docs/_static/tutorial_11.png b/docs/_static/tutorial_11.png new file mode 100644 index 0000000..13e17ca Binary files /dev/null and b/docs/_static/tutorial_11.png differ diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 9ee917f..32200d0 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -23,3 +23,4 @@ Before attempting the tutorials, make sure to review :ref:`getting_started` and, tutorial_8 tutorial_9 tutorial_10 + tutorial_11 diff --git a/docs/tutorials/tutorial_11.rst b/docs/tutorials/tutorial_11.rst new file mode 100644 index 0000000..e1e37ec --- /dev/null +++ b/docs/tutorials/tutorial_11.rst @@ -0,0 +1,129 @@ +.. Contains the eleventh tutorial. +.. _tutorial_11: + +Tutorial 11 - Solving the k-epsilon RANS Turbulence Model (2D) +================================================================ + +The files for this tutorial can be found in "examples/backward_facing_step". + +Governing Equations +-------------------- + +This tutorial demonstrates OpenCMP's high-Reynolds-number k-epsilon RANS turbulence model. The mean-flow momentum equation gains a turbulent (eddy) viscosity :math:`\nu_t` on top of the molecular kinematic viscosity :math:`\nu`, and two transport equations close the model -- turbulent kinetic energy :math:`k` and its dissipation rate :math:`\epsilon`: + +.. math:: + \frac{\partial \bm{u}}{\partial t} + \bm{\nabla} \cdot \left( \bm{u} \bm{w} \right) - \bm{\nabla} \cdot \left[ \left( \nu + \nu_t \right) \bm{\nabla} \bm{u} \right] + \bm{\nabla} p &= \bm{f} \mbox{ in } \Omega \\ + \bm{\nabla} \cdot \bm{u} &= 0 \mbox{ in } \Omega \\ + \frac{\partial k}{\partial t} + \bm{\nabla} \cdot \left( \bm{u} k \right) - \bm{\nabla} \cdot \left[ \left( \nu + \frac{\nu_t}{\sigma_k} \right) \bm{\nabla} k \right] &= P_k - \epsilon \\ + \frac{\partial \epsilon}{\partial t} + \bm{\nabla} \cdot \left( \bm{u} \epsilon \right) - \bm{\nabla} \cdot \left[ \left( \nu + \frac{\nu_t}{\sigma_\epsilon} \right) \bm{\nabla} \epsilon \right] &= C_1 \frac{\epsilon}{k} P_k - C_2 \frac{\epsilon^2}{k} + +with production :math:`P_k = 2 \nu_t \, \bm{S} : \bm{S}`, :math:`\bm{S} = \frac{1}{2} \left( \bm{\nabla} \bm{u} + \bm{\nabla} \bm{u}^T \right)`, and the algebraic closure :math:`\nu_t = C_\mu k^2 / \epsilon`. + +The standard k-epsilon model is only valid in the fully turbulent log layer, so instead of resolving the viscous sublayer with mesh refinement, OpenCMP applies a wall function on the boundary marked ``wall``: near-wall cells get an algebraic wall-law eddy viscosity and dissipation, while velocity keeps its ordinary no-slip condition. Controlled by the ``[OTHER]`` switches ``wall_function`` and ``wall_boundary``. + +The example is a 2D backward-facing step: a duct of height :math:`H = 1` that abruptly expands to height :math:`2H`, generated by "backward_facing_step.geo". ``wall`` covers the upstream floor, the step face, and the downstream floor/ceiling; ``inlet`` and ``outlet`` are the two open ends. Kinematic viscosity is 0.00002, inlet velocity is a uniform 1.0 in x. + +The Main Configuration Files +------------------------------ + +"config_IC" runs a Stokes solve for a divergence-free initial guess, saving velocity and pressure to separate ".sol" files (``split_components = True``) so they can be reloaded individually. "config" runs the main k-epsilon solve. + +"config" adds the ``k`` and ``epsilon`` finite elements (discontinuous ``L2``, required for the wall-law dissipation) and switches the model to ``KEpsilonINS``:: + + [FINITE ELEMENT SPACE] + elements = u -> HDiv + p -> L2 + k -> L2 + epsilon -> L2 + interpolant_order = 3 + + [SOLVER] + linearization_method = Oseen + nonlinear_solver = NoMixing + nonlinear_tolerance = relative -> 1e-4 + absolute -> 1e-4 + nonlinear_max_iterations = 200 + relaxation_factors = 0.5, 0.5, 0.3, 0.3 + + [OTHER] + model = KEpsilonINS + wall_function = True + wall_boundary = wall + production_limiter = True + auto_turbulence_inlet = inlet + +Four relaxation factors are given, one per component (:math:`\bm{u}`, :math:`p`, :math:`k`, :math:`\epsilon`); under-relaxing keeps the segregated k-epsilon iteration from diverging while :math:`\nu_t` is still adjusting. ``auto_turbulence_inlet`` is new -- see below. + +The Boundary and Initial Condition Configuration Files +----------------------------------------------------------- + +Velocity is no-slip on ``wall`` and uniform at ``inlet``, with a "do-nothing" stress condition at ``outlet``. No values are given for :math:`k` or :math:`\epsilon`, only a zero-flux (Neumann) fallback:: + + [DIRICHLET] + u = inlet -> [1.0, 0.0] + wall -> [0.0, 0.0] + + [STRESS] + u = outlet -> [0.0, 0.0] + + [NEUMANN] + k = outlet -> 0.0 + wall -> 0.0 + epsilon = outlet -> 0.0 + wall -> 0.0 + +Similarly, "ic_dir/ic_config" only supplies :math:`\bm{u}` and :math:`p` (reloaded from the Stokes solve); :math:`k` and :math:`\epsilon` are left unset:: + + [KEpsilonINS] + u = all -> output/components_sol/u.sol + p = all -> output/components_sol/p.sol + +Setting ``auto_turbulence_inlet = inlet`` in "config" tells ``KEpsilonINS`` to fill in the missing inlet Dirichlet values (and initial condition, since none were explicitly given) automatically from the prescribed inlet velocity: it computes the bulk velocity through that boundary, then the standard internal-flow estimate + +.. math:: + Re = \frac{U_{bulk} D_H}{\nu}, \quad I = 0.16\, Re^{-1/8}, \quad k = \frac{3}{2}(U_{bulk} I)^2, \quad \epsilon = \frac{C_\mu^{3/4} k^{3/2}}{r \, D_H} + +where the hydraulic diameter :math:`D_H` and length-scale ratio :math:`r` come from the model configuration file. Any value the user *does* specify explicitly in ``bc_config``/``ic_config`` still takes precedence. + +The Model Configuration File +-------------------------------- + +The k-epsilon closure constants and numerical safeguards are all shown at their default values (omitting them changes nothing). ``turbulence_hydraulic_diameter`` and ``turbulence_length_scale_ratio`` feed the automatic inlet formula above; for this 2D channel the hydraulic diameter is :math:`2H = 2`:: + + [PARAMETERS] + kinematic_viscosity = all -> 0.00002 + c_mu = all -> 0.09 + c_1 = all -> 1.44 + c_2 = all -> 1.92 + sigma_k = all -> 1.0 + sigma_epsilon = all -> 1.3 + kappa = all -> 0.4187 + e_log = all -> 9.793 + k_floor = all -> 1e-8 + epsilon_floor = all -> 1e-6 + max_viscosity_ratio = all -> 200000 + production_limit_coefficient = all -> 10 + max_epsilon_k_ratio = all -> 10 + turbulence_hydraulic_diameter = all -> 2.0 + turbulence_length_scale_ratio = all -> 0.07 + + [FUNCTIONS] + source = u -> [0.0, 0.0] + k -> 0.0 + epsilon -> 0.0 + + +Running the Simulation +-------------------------- + +From "examples/backward_facing_step": + +1) Run the Stokes solve: :code:`python3 -m opencmp config_IC` +2) Run the k-epsilon solve: :code:`python3 -m opencmp config` + +Since ``transient = False``, this is a steady-state solve driven by Picard (Oseen) iteration rather than time steps. The result shows the expected recirculation bubble just downstream of the step: + +.. image:: ../_static/tutorial_11.png + :width: 700 + :align: center + :alt: Steady-state velocity magnitude for the backward-facing-step k-epsilon solve. diff --git a/examples/backward_facing_step/backward_facing_step.geo b/examples/backward_facing_step/backward_facing_step.geo new file mode 100644 index 0000000..f296bef --- /dev/null +++ b/examples/backward_facing_step/backward_facing_step.geo @@ -0,0 +1,51 @@ +// Adjustable coarse/fine proxy for the backward-facing-step tutorial. +// Generate with: gmsh -2 backward_facing_step_proxy.geo -format msh2 + +SetFactory("Built-in"); + +H = 1.0; +L_up = 2.0; +L_down = 6.0; + +// Primary mesh controls. Reduce these values for a finer mesh. +h_wall = 0.2; +h_core = 0.2; + +Point(1) = {-L_up, H, 0, h_wall}; +Point(2) = {0, H, 0, h_wall}; +Point(3) = {0, 0, 0, h_wall}; +Point(4) = {L_down, 0, 0, h_wall}; +Point(5) = {L_down, 2 * H, 0, h_wall}; +Point(6) = {-L_up, 2 * H, 0, h_wall}; + +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 5}; +Line(5) = {5, 6}; +Line(6) = {6, 1}; + +Curve Loop(1) = {1, 2, 3, 4, 5, 6}; +Plane Surface(1) = {1}; + +Physical Curve("inlet") = {6}; +Physical Curve("outlet") = {4}; +Physical Curve("wall") = {1, 2, 3, 5}; +Physical Surface("surface") = {1}; + +Field[1] = Distance; +Field[1].CurvesList = {1, 2, 3, 5}; +Field[1].Sampling = 250; + +Field[2] = Threshold; +Field[2].InField = 1; +Field[2].SizeMin = h_wall; +Field[2].SizeMax = h_core; +Field[2].DistMin = 0.10 * H; +Field[2].DistMax = 0.50 * H; +Background Field = 2; + +Mesh.CharacteristicLengthExtendFromBoundary = 0; +Mesh.Algorithm = 6; +Mesh.ElementOrder = 1; +Mesh.MshFileVersion = 2.2; diff --git a/examples/backward_facing_step/backward_facing_step.msh b/examples/backward_facing_step/backward_facing_step.msh new file mode 100644 index 0000000..d927dbb --- /dev/null +++ b/examples/backward_facing_step/backward_facing_step.msh @@ -0,0 +1,1433 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +4 +1 1 "inlet" +1 2 "outlet" +1 3 "wall" +2 4 "surface" +$EndPhysicalNames +$Nodes +473 +1 -2 1 0 +2 0 1 0 +3 0 0 0 +4 6 0 0 +5 6 2 0 +6 -2 2 0 +7 -1.799999999999167 1 0 +8 -1.6 1 0 +9 -1.400000000001387 1 0 +10 -1.200000000002774 1 0 +11 -1.000000000004117 1 0 +12 -0.8000000000033287 1 0 +13 -0.6000000000024965 1 0 +14 -0.4000000000016644 1 0 +15 -0.200000000000832 1 0 +16 0 0.7999999999999998 0 +17 0 0.6000000000013869 0 +18 0 0.4000000000016644 0 +19 0 0.2000000000008322 0 +20 0.1999999999991098 0 0 +21 0.3999999999982407 0 0 +22 0.5999999999975234 0 0 +23 0.7999999999967187 0 0 +24 0.9999999999956506 0 0 +25 1.199999999994909 0 0 +26 1.399999999994197 0 0 +27 1.599999999993484 0 0 +28 1.799999999992772 0 0 +29 1.999999999992059 0 0 +30 2.199999999991346 0 0 +31 2.399999999990634 0 0 +32 2.599999999989921 0 0 +33 2.799999999989208 0 0 +34 2.999999999988583 0 0 +35 3.199999999989264 0 0 +36 3.399999999990031 0 0 +37 3.599999999990798 0 0 +38 3.799999999991565 0 0 +39 3.999999999992331 0 0 +40 4.199999999993098 0 0 +41 4.399999999993865 0 0 +42 4.599999999994631 0 0 +43 4.799999999995399 0 0 +44 4.999999999996166 0 0 +45 5.199999999996932 0 0 +46 5.399999999997699 0 0 +47 5.599999999998466 0 0 +48 5.799999999999233 0 0 +49 6 0.1999999999996293 0 +50 6 0.3999999999991157 0 +51 6 0.5999999999985328 0 +52 6 0.7999999999979498 0 +53 6 0.9999999999973885 0 +54 6 1.199999999997894 0 +55 6 1.39999999999842 0 +56 6 1.599999999998947 0 +57 6 1.799999999999474 0 +58 5.799999999999445 2 0 +59 5.599999999998889 2 0 +60 5.399999999998333 2 0 +61 5.199999999997778 2 0 +62 4.999999999997396 2 0 +63 4.799999999998888 2 0 +64 4.600000000000552 2 0 +65 4.400000000002217 2 0 +66 4.200000000003882 2 0 +67 4.00000000000546 2 0 +68 3.800000000006102 2 0 +69 3.600000000006657 2 0 +70 3.400000000007211 2 0 +71 3.200000000007766 2 0 +72 3.000000000008321 2 0 +73 2.800000000008876 2 0 +74 2.60000000000943 2 0 +75 2.400000000009986 2 0 +76 2.200000000010541 2 0 +77 2.000000000011009 2 0 +78 1.800000000010541 2 0 +79 1.600000000009986 2 0 +80 1.400000000009431 2 0 +81 1.200000000008877 2 0 +82 1.000000000008322 2 0 +83 0.8000000000077669 2 0 +84 0.6000000000072117 2 0 +85 0.4000000000066581 2 0 +86 0.2000000000061029 2 0 +87 5.547562409446982e-12 2 0 +88 -0.199999999995006 2 0 +89 -0.3999999999955612 2 0 +90 -0.5999999999961165 2 0 +91 -0.7999999999966718 2 0 +92 -0.9999999999972262 2 0 +93 -1.199999999997781 2 0 +94 -1.399999999998336 2 0 +95 -1.59999999999889 2 0 +96 -1.799999999999445 2 0 +97 -2 1.8 0 +98 -2 1.6 0 +99 -2 1.4 0 +100 -2 1.2 0 +101 4.715889888994249 1.820945090216231 0 +102 4.126254593927985 1.80991484138901 0 +103 2.29999999999099 0.1732050807562706 0 +104 2.898220202048412 0.1770020610447061 0 +105 1.076568448514803 0.1869321149206973 0 +106 1.725488707240545 0.1980115266239818 0 +107 3.482343623972824 1.809397822998246 0 +108 2.49815492179354 1.851001200578305 0 +109 1.905982299339834 1.864155400698314 0 +110 1.296315300350435 1.837012656419556 0 +111 3.899455533390836 0.1692275483986728 0 +112 4.501353593709718 0.2020259114662162 0 +113 5.099999999996546 0.1732050807575517 0 +114 5.826794919243029 1.099999999997815 0 +115 0.1732050807562872 0.8999999999999999 0 +116 0.6965691307747465 1.82072802706088 0 +117 5.303637951712282 1.820333885191719 0 +118 0.1732050807567705 1.099999999999797 0 +119 0.3464101615134161 0.9999999999994565 0 +120 3.299999999989882 0.1732050807574185 0 +121 5.825826487855652 0.6975213828789675 0 +122 2.891905155833481 1.815944470027591 0 +123 0.6981596001939732 0.1777086041428091 0 +124 -1.489686081938773 1.163549923920361 0 +125 -1.100000000003446 1.173205080755705 0 +126 -0.2999835968974923 1.825730211136559 0 +127 -0.6999999999962899 1.826794919242821 0 +128 0.3219644455985776 1.809667120872487 0 +129 0.3464101615128699 0.7999999999994878 0 +130 0.5196152422700365 0.8999999999989881 0 +131 0.5196152422705889 1.099999999998985 0 +132 0.6928203230272075 0.9999999999984999 0 +133 -1.299999999997694 1.826794919242693 0 +134 -0.5043371746586531 1.16071004212909 0 +135 5.836293387951234 1.488665085739994 0 +136 5.513068339144769 0.1981959154845498 0 +137 0.6928203230277769 1.199999999998502 0 +138 0.8660254037843816 1.099999999998013 0 +139 0.8660254037838151 0.8999999999980097 0 +140 1.039230484540986 0.9999999999975215 0 +141 1.039230484541552 1.199999999997523 0 +142 1.212435565298157 1.099999999997034 0 +143 1.212435565297593 0.8999999999970348 0 +144 1.388163382683139 1.004369508010981 0 +145 1.386061102160056 1.200728251332285 0 +146 1.559336258934117 1.100849626554419 0 +147 1.561870674898172 0.9052393637726254 0 +148 1.732636720936703 1.001014831718056 0 +149 1.732652617810468 0.8010423659120773 0 +150 1.905453842260253 0.9003428662682451 0 +151 1.905386532876253 1.100226282994355 0 +152 2.078515735496588 1.000094858206663 0 +153 2.078491870910347 1.200053523529815 0 +154 2.251680327879773 1.100024730285317 0 +155 2.251677557248011 0.9000199314113152 0 +156 2.424875428170985 1.000007443611683 0 +157 2.42487376475952 0.8000045624994914 0 +158 2.598077366642508 0.9000020010137799 0 +159 2.598077120164044 1.100001574099576 0 +160 2.771281636126631 1.000000595847151 0 +161 1.905389182354179 0.7002308720266913 0 +162 2.77128150091502 1.200000361652794 0 +163 2.944486465003926 1.100000159577923 0 +164 2.944486445558745 0.9000001258988586 0 +165 3.117691481095167 1.000000047573735 0 +166 3.11769147355894 1.20000003451963 0 +167 3.290896542281853 1.100000013676174 0 +168 3.290896540275674 0.9000000102023461 0 +169 3.464101617436827 1.000000003973375 0 +170 3.458948662971378 0.8089251796816016 0 +171 3.636447870916344 0.9014875306024608 0 +172 3.637163558781539 1.100247922422684 0 +173 3.810344782969403 1.000289242163832 0 +174 3.810460088185946 1.20008952742414 0 +175 3.98368041038376 1.100063128257312 0 +176 3.9836829506234 0.9000587283962558 0 +177 4.156921938164967 0.9999999999887335 0 +178 4.156916287033673 0.8000097880565726 0 +179 4.330126077066453 0.9000016313328891 0 +180 4.330126861946278 1.100000271879021 0 +181 4.503331916540244 1.000000317193668 0 +182 4.503331912179247 0.8000003247461893 0 +183 4.67653718043535 0.8999999999872841 0 +184 4.676537149912823 1.100000052854934 0 +185 4.849742256105335 1.000000008798077 0 +186 4.849742255258027 1.200000010266623 0 +187 5.022947341471486 1.100000000814582 0 +188 5.022947341021563 0.9000000015929113 0 +189 4.676537149184971 0.7000000541136884 0 +190 1.548944206232109 0.7234994822369983 0 +191 5.188039063141041 1.000000000391674 0 +192 5.194800195956938 0.8000000003212473 0 +193 5.186686836701978 1.200000000191548 0 +194 1.558885249764171 1.300068455756421 0 +195 2.598076842927763 0.7000010939142114 0 +196 3.617104269811613 0.6837066347770926 0 +197 0.004465819873861432 1.207735026918892 0 +198 2.251673579818007 1.300013042298507 0 +199 0.5196152422711537 1.299999999998988 0 +200 4.676537174360107 1.300000010511861 0 +201 3.29089653902067 1.300000008026661 0 +202 3.286402000738564 0.7119429741441921 0 +203 0.5196152422694704 0.6999999999989922 0 +204 3.983702168160518 1.300025442606306 0 +205 4.849742259191726 1.400000003454201 0 +206 -1.836293387951482 1.488665085741304 0 +207 3.810500713699945 1.400019161664957 0 +208 5.372590891217914 0.9210870397500105 0 +209 5.369671030296724 0.7035145066687083 0 +210 3.117691457720441 1.400000007085566 0 +211 2.078467374384346 1.400011094301195 0 +212 4.503332063219504 0.6000000631352367 0 +213 4.673233096278839 0.5057228416170048 0 +214 2.424871674884895 0.6000009427314389 0 +215 0.6928203230283385 1.399999999998502 0 +216 2.598076262366138 0.5000000883529586 0 +217 5.195979306052785 0.6005857511556347 0 +218 1.385717309231959 1.400132784512563 0 +219 1.558865091167685 1.50003354004212 0 +220 0.1687392608823997 0.307735026920587 0 +221 0.1739493840693875 1.301289171152911 0 +222 0.006604759422929646 1.4002133391062 0 +223 3.983710264421383 0.7000114194016973 0 +224 4.156919897477751 0.6000035345687855 0 +225 2.598076397623053 1.30000032262089 0 +226 2.771281405873768 0.6000001970394528 0 +227 2.771281309251153 1.400000029680259 0 +228 2.251666654418516 1.500001047153748 0 +229 0.5196152422717312 1.49999999999899 0 +230 3.290896534760902 1.500000000647548 0 +231 5.829907236381057 0.3192180115220712 0 +232 -0.1810563553985545 1.312673029721057 0 +233 5.699152526142959 1.830645590988407 0 +234 -1.699152526142983 1.830645590988467 0 +235 5.371887606885953 0.5050251110450485 0 +236 0.1768500303165584 1.496058672223337 0 +237 0.004627693960895851 1.623460803044036 0 +238 3.985881955872225 1.496249939477035 0 +239 4.15728033970181 1.399379230339546 0 +240 3.80740425487539 1.628539324165006 0 +241 4.849191575290535 0.6009538159461567 0 +242 1.732060622120716 1.400016999296639 0 +243 3.464101615974917 1.400000001439321 0 +244 0.346410161512376 0.5999999999997115 0 +245 0.51961524226905 0.4999999999992684 0 +246 0.6928203230260978 0.5999999999985435 0 +247 4.673233106221126 1.494277180237408 0 +248 4.502781419631145 1.399046198450212 0 +249 4.4988302178202 1.597630227567447 0 +250 2.598396876698366 1.496448710809146 0 +251 4.830149691825156 0.3730921704639418 0 +252 3.983716857406679 0.4999999999892177 0 +253 4.17384675532625 0.3717704213845126 0 +254 3.120686757099165 1.594811984774291 0 +255 4.832908144439959 1.630585413706644 0 +256 5.020141655491436 1.505097569517528 0 +257 5.012752367952463 1.682102818243168 0 +258 5.534923163264589 0.6014232696087509 0 +259 2.085357086957448 1.598534757365323 0 +260 2.25829019712397 1.687746082124299 0 +261 1.731003674176607 1.598901204373453 0 +262 1.549628633468017 1.688497413455346 0 +263 5.198691634203925 1.628840110285881 0 +264 2.781218271986689 1.626678849037754 0 +265 2.789740053531362 0.3665442331342964 0 +266 2.951965145491068 0.5071607391810755 0 +267 2.601088089271256 0.3130596481368549 0 +268 3.295953330767563 1.695123058875572 0 +269 5.197739293537637 0.3717561248405304 0 +270 4.156921938166652 1.599999999988729 0 +271 5.532395030868232 0.4085446626195848 0 +272 5.708580577717816 0.5048245507830141 0 +273 4.502364804740366 0.4125866724945864 0 +274 0.6934451243197839 1.603454671175565 0 +275 0.8661295373339907 1.500575778527523 0 +276 0.5172046889956422 1.699992421214843 0 +277 0.6942940546944972 0.3949652625575458 0 +278 0.8646278532758763 0.497667928336997 0 +279 0.5196152422683276 0.2999999999989925 0 +280 3.110143347565837 1.805897137066593 0 +281 0.8791012364196535 0.2977909622138311 0 +282 1.046543270869362 0.4066990725784249 0 +283 1.237137711463851 0.3352183737150695 0 +284 1.200173951640256 0.5277586787804298 0 +285 1.384888354527269 0.4570447909359722 0 +286 1.405237467599955 0.3004054321434585 0 +287 1.580447146553319 0.3791446391038401 0 +288 1.37221104299095 0.6303266725864406 0 +289 0.8660254037860764 1.699999999998013 0 +290 1.037175245918001 1.596506124600052 0 +291 1.060946716250862 1.807557473227634 0 +292 5.299999999997024 0.1732050807573886 0 +293 0.3464101615151304 1.59999999999948 0 +294 -0.4999808630483291 1.826607997558866 0 +295 -0.5999804074068338 1.65354921453885 0 +296 -0.7999967345646268 1.65358306782786 0 +297 -0.7026574683113098 1.481913348189344 0 +298 -0.9004423671432424 1.48391331892005 0 +299 -0.5039918885103523 1.485905478777195 0 +300 -1.000073183615576 1.65417680357446 0 +301 -1.1000223670117 1.48226679975507 0 +302 -0.8058391952123173 1.321099363762176 0 +303 -1.200015925102734 1.654001339671436 0 +304 -1.300001661288311 1.4821862652492 0 +305 -1.49985875435587 1.827436697867106 0 +306 -1.603527600612199 1.670674566586002 0 +307 5.499999999998385 1.826794919242762 0 +308 5.603551141552696 1.670567603481963 0 +309 5.399999999998276 1.653589838485846 0 +310 5.509837670829932 1.481629091635636 0 +311 -1.400567323558871 1.656913091057229 0 +312 -1.515789822106204 1.484589124741159 0 +313 5.830914871827806 1.298110847621704 0 +314 5.665501031863917 1.201317489887759 0 +315 5.656804168447166 1.004395628983714 0 +316 5.480384757729246 1.099999999997537 0 +317 3.499999999990197 0.1732050807574295 0 +318 3.399999999989999 0.3446712196452 0 +319 3.595020984584157 0.3432457019849392 0 +320 3.19809859783542 0.3294471628575002 0 +321 3.502572509313167 0.5061206063945867 0 +322 3.699079419656267 0.1720147487760303 0 +323 3.807806239947121 0.3727714847803223 0 +324 3.099999999989636 0.1732050807579082 0 +325 4.287438443686776 0.1800726166462276 0 +326 4.094937372947683 0.173727112810578 0 +327 4.714702671454938 0.1804433491965457 0 +328 4.908969241554524 0.1751603403515806 0 +329 0.1000000000058252 1.826794919242632 0 +330 -0.09922598381920536 1.820540724211511 0 +331 -0.2064110937829766 1.649799302264625 0 +332 -0.4025907161289083 1.654144571069703 0 +333 -0.9000116530278772 1.826891618314632 0 +334 0.1732050807560454 0.7000000000004494 0 +335 -1.100016793623148 1.826977446800536 0 +336 0.3465342120659792 1.200214861858269 0 +337 5.022947341948033 0.4999999999863232 0 +338 0.8632776245558037 0.6996113213878521 0 +339 0.3471623787516567 1.399593784205329 0 +340 4.330125542649516 0.7000025569713202 0 +341 2.412933764763291 0.3683969956771279 0 +342 2.251666049837786 0.4999999999940998 0 +343 3.812162552631986 0.5931428276500698 0 +344 3.116942382359428 0.801990526300487 0 +345 3.114731248084986 0.613446379672641 0 +346 4.330127018920472 0.4999999999882476 0 +347 5.523404169412773 0.7876708026657232 0 +348 2.25166849748219 0.7000042394375872 0 +349 0.6923623598220716 0.7999352202301475 0 +350 1.03068584204855 0.7998131326163543 0 +351 0.8660427593763594 1.300095963086265 0 +352 2.420359064926746 1.583439179256986 0 +353 2.944486372868165 1.499999999992146 0 +354 1.909719070115393 1.493540023583126 0 +355 1.387347341110679 1.584783090323297 0 +356 -1.830617856523183 1.300113576604748 0 +357 1.732200382276675 1.200259070988611 0 +358 1.906054308780622 1.299017832448957 0 +359 1.038908192961318 1.399529644367572 0 +360 2.771281678332552 0.8000006689509324 0 +361 2.42487566404218 1.200007852151597 0 +362 2.424174706254478 1.396651692381812 0 +363 1.7069451926657 0.5788699886331586 0 +364 1.905255888324014 0.4999999999950787 0 +365 2.945114717616757 0.7037664395072412 0 +366 1.212464703248721 1.300065113367335 0 +367 3.807242115229021 0.7964493971652346 0 +368 5.022700195492663 0.7001832820186 0 +369 1.210268798745662 1.493246287708919 0 +370 2.944486429886274 1.300000098751386 0 +371 2.078527630653733 0.800115461224182 0 +372 5.188384454300706 1.405512631184641 0 +373 5.019607467069226 1.301768369238187 0 +374 5.369357503464008 1.299999999985368 0 +375 4.156973122880004 1.19991134550986 0 +376 4.849609282921735 0.8001895270761195 0 +377 3.464077761565158 1.2000413249214 0 +378 2.078494702955419 0.6000584287787054 0 +379 3.637268405683868 1.300066322977115 0 +380 3.637306695895707 1.499999999990189 0 +381 4.503240256885547 1.199841141812989 0 +382 2.07553849000577 0.4050813584976872 0 +383 1.906405646275527 0.3043419355987299 0 +384 2.085278805223078 0.194534435816369 0 +385 4.330088169994672 1.29969636466331 0 +386 3.464101615139101 1.599999999990666 0 +387 4.329338184039713 1.499292003499575 0 +388 4.32422419753062 1.699896762712545 0 +389 5.832210535285554 0.909051289042383 0 +390 3.700000000006382 1.826794919243588 0 +391 2.700000000009153 1.826794919243593 0 +392 3.900000000005778 1.826794919244212 0 +393 -1.30000000000208 1.173205080755686 0 +394 4.499379578224482 1.80704129798149 0 +395 0.1724607774438955 0.5012891711541198 0 +396 0.3455418076480208 0.4015040330127715 0 +397 0.3257580674725927 0.1842319739517271 0 +398 2.096542765454291 1.810997910731447 0 +399 1.708938321142409 1.822768918122662 0 +400 -0.3051964471274664 1.483275222142123 0 +401 3.28908283385365 0.515654550574122 0 +402 -0.3066178263334563 1.160554648292915 0 +403 5.811754837909366 1.683609464896652 0 +404 -1.684375567857481 1.1891186510725 0 +405 5.69019237886328 0.1675426480548286 0 +406 -1.811751474918005 1.68362474534067 0 +407 4.907208104747131 1.830352619829741 0 +408 1.028839573816394 0.6036687628588496 0 +409 -1.684116286406766 1.386340244868585 0 +410 5.677459108719409 1.386885790363416 0 +411 -0.9009731992036432 1.162448200375055 0 +412 -0.6028455809324694 1.321657870539575 0 +413 -1.392233613016618 1.32401370156327 0 +414 -1.199999999995738 1.326794919244458 0 +415 -0.7023325250021514 1.160985912800983 0 +416 -0.9999999999957386 1.326829224317669 0 +417 -0.4006742121601586 1.320796048600325 0 +418 2.999999999988044 0.3260175575665244 0 +419 2.712033120063154 0.1646097622315301 0 +420 2.492780539661815 0.1660605761080349 0 +421 1.287445451722304 0.1570626836232561 0 +422 1.762491580749971 0.3943382985636837 0 +423 1.217813847779629 1.677603286483487 0 +424 2.564509777687912 1.687018363635503 0 +425 3.995998173949687 1.67482481437328 0 +426 3.994294093967127 0.319970654152817 0 +427 5.009743497271311 0.3227692607221974 0 +428 0.1899704662793975 1.671196303076394 0 +429 5.114471064203338 1.834128597052813 0 +430 5.383412242583139 0.3298325302468211 0 +431 1.400777065824461 0.8262170480831602 0 +432 3.618231237977881 1.672946413277539 0 +433 1.207520506936367 0.7146307158203782 0 +434 1.500000000009709 1.860895300726176 0 +435 3.668497031892817 0.5010562147965769 0 +436 2.2340974415139 0.3304297222285669 0 +437 2.969687981070667 1.668666488179675 0 +438 2.301818574521045 1.858905485932389 0 +439 5.332682841707856 1.497778316214054 0 +440 4.300000000003051 1.865153937878145 0 +441 4.359026123276716 0.3332911243959581 0 +442 0.500000000006935 1.866025403784319 0 +443 0.5022663947261852 0.1335616222989274 0 +444 4.647666631734427 1.668829334957335 0 +445 4.647889906346783 0.3346492586031081 0 +446 3.298820439016035 1.864260717778096 0 +447 1.902453308386973 0.1376638065582401 0 +448 -0.1387866474061928 1.492255836565508 0 +449 1.920986368890174 1.672723820525664 0 +450 0.8999999999961846 0.1305125270622566 0 +451 0.8856248965360678 1.864636067908769 0 +452 5.538223072832562 1.298488650756656 0 +453 1.511131962955855 0.1661347467123421 0 +454 5.875475075227786 0.4999999999988243 0 +455 5.658678150342996 0.6919758183291449 0 +456 5.319411810450838 1.104217408063228 0 +457 1.518687188593869 0.5537771146992819 0 +458 0.1464101615141686 0.1464101615146301 0 +459 -1.853589838486143 1.146410161513634 0 +460 -1.853589838485723 1.853589838485872 0 +461 5.853589838485812 1.853589838485814 0 +462 -0.1464101615138812 1.146410161513663 0 +463 5.853589838486437 0.1464101615132628 0 +464 5.66437825603305 0.32010994793523 0 +465 3.110775565050634 0.4583452779703726 0 +466 5.67586250699043 0.845917921811135 0 +467 5.523063515482079 0.9418530826340643 0 +468 3.433183893948602 0.6535510639423447 0 +469 -1.547924112768604 1.313756895468756 0 +470 1.410352789291488 1.73849816706357 0 +471 5.689696971232865 1.541798135581961 0 +472 -1.689696971233646 1.541798135583045 0 +473 2.398979008517908 1.735642598100635 0 +$EndNodes +$Elements +944 +1 1 2 3 1 1 7 +2 1 2 3 1 7 8 +3 1 2 3 1 8 9 +4 1 2 3 1 9 10 +5 1 2 3 1 10 11 +6 1 2 3 1 11 12 +7 1 2 3 1 12 13 +8 1 2 3 1 13 14 +9 1 2 3 1 14 15 +10 1 2 3 1 15 2 +11 1 2 3 2 2 16 +12 1 2 3 2 16 17 +13 1 2 3 2 17 18 +14 1 2 3 2 18 19 +15 1 2 3 2 19 3 +16 1 2 3 3 3 20 +17 1 2 3 3 20 21 +18 1 2 3 3 21 22 +19 1 2 3 3 22 23 +20 1 2 3 3 23 24 +21 1 2 3 3 24 25 +22 1 2 3 3 25 26 +23 1 2 3 3 26 27 +24 1 2 3 3 27 28 +25 1 2 3 3 28 29 +26 1 2 3 3 29 30 +27 1 2 3 3 30 31 +28 1 2 3 3 31 32 +29 1 2 3 3 32 33 +30 1 2 3 3 33 34 +31 1 2 3 3 34 35 +32 1 2 3 3 35 36 +33 1 2 3 3 36 37 +34 1 2 3 3 37 38 +35 1 2 3 3 38 39 +36 1 2 3 3 39 40 +37 1 2 3 3 40 41 +38 1 2 3 3 41 42 +39 1 2 3 3 42 43 +40 1 2 3 3 43 44 +41 1 2 3 3 44 45 +42 1 2 3 3 45 46 +43 1 2 3 3 46 47 +44 1 2 3 3 47 48 +45 1 2 3 3 48 4 +46 1 2 2 4 4 49 +47 1 2 2 4 49 50 +48 1 2 2 4 50 51 +49 1 2 2 4 51 52 +50 1 2 2 4 52 53 +51 1 2 2 4 53 54 +52 1 2 2 4 54 55 +53 1 2 2 4 55 56 +54 1 2 2 4 56 57 +55 1 2 2 4 57 5 +56 1 2 3 5 5 58 +57 1 2 3 5 58 59 +58 1 2 3 5 59 60 +59 1 2 3 5 60 61 +60 1 2 3 5 61 62 +61 1 2 3 5 62 63 +62 1 2 3 5 63 64 +63 1 2 3 5 64 65 +64 1 2 3 5 65 66 +65 1 2 3 5 66 67 +66 1 2 3 5 67 68 +67 1 2 3 5 68 69 +68 1 2 3 5 69 70 +69 1 2 3 5 70 71 +70 1 2 3 5 71 72 +71 1 2 3 5 72 73 +72 1 2 3 5 73 74 +73 1 2 3 5 74 75 +74 1 2 3 5 75 76 +75 1 2 3 5 76 77 +76 1 2 3 5 77 78 +77 1 2 3 5 78 79 +78 1 2 3 5 79 80 +79 1 2 3 5 80 81 +80 1 2 3 5 81 82 +81 1 2 3 5 82 83 +82 1 2 3 5 83 84 +83 1 2 3 5 84 85 +84 1 2 3 5 85 86 +85 1 2 3 5 86 87 +86 1 2 3 5 87 88 +87 1 2 3 5 88 89 +88 1 2 3 5 89 90 +89 1 2 3 5 90 91 +90 1 2 3 5 91 92 +91 1 2 3 5 92 93 +92 1 2 3 5 93 94 +93 1 2 3 5 94 95 +94 1 2 3 5 95 96 +95 1 2 3 5 96 6 +96 1 2 1 6 6 97 +97 1 2 1 6 97 98 +98 1 2 1 6 98 99 +99 1 2 1 6 99 100 +100 1 2 1 6 100 1 +101 2 2 4 1 261 449 399 +102 2 2 4 1 81 291 110 +103 2 2 4 1 399 449 109 +104 2 2 4 1 327 328 251 +105 2 2 4 1 255 407 101 +106 2 2 4 1 253 326 325 +107 2 2 4 1 122 391 264 +108 2 2 4 1 270 388 102 +109 2 2 4 1 101 444 255 +110 2 2 4 1 251 445 327 +111 2 2 4 1 105 282 281 +112 2 2 4 1 287 453 106 +113 2 2 4 1 240 392 390 +114 2 2 4 1 265 419 104 +115 2 2 4 1 113 292 269 +116 2 2 4 1 329 330 237 +117 2 2 4 1 103 420 341 +118 2 2 4 1 104 418 265 +119 2 2 4 1 252 343 323 +120 2 2 4 1 325 441 253 +121 2 2 4 1 117 429 263 +122 2 2 4 1 128 293 276 +123 2 2 4 1 106 422 287 +124 2 2 4 1 388 440 102 +125 2 2 4 1 291 423 110 +126 2 2 4 1 268 386 107 +127 2 2 4 1 82 291 81 +128 2 2 4 1 143 433 431 +129 2 2 4 1 264 437 122 +130 2 2 4 1 281 450 105 +131 2 2 4 1 341 436 103 +132 2 2 4 1 263 309 117 +133 2 2 4 1 330 331 237 +134 2 2 4 1 26 453 421 +135 2 2 4 1 383 422 106 +136 2 2 4 1 21 397 20 +137 2 2 4 1 431 433 288 +138 2 2 4 1 286 453 287 +139 2 2 4 1 343 435 323 +140 2 2 4 1 103 384 30 +141 2 2 4 1 276 442 128 +142 2 2 4 1 390 432 240 +143 2 2 4 1 363 457 287 +144 2 2 4 1 287 422 363 +145 2 2 4 1 421 453 286 +146 2 2 4 1 263 429 257 +147 2 2 4 1 122 280 72 +148 2 2 4 1 105 283 282 +149 2 2 4 1 292 430 269 +150 2 2 4 1 30 384 29 +151 2 2 4 1 72 280 71 +152 2 2 4 1 237 428 329 +153 2 2 4 1 116 276 274 +154 2 2 4 1 277 279 123 +155 2 2 4 1 398 449 259 +156 2 2 4 1 269 427 113 +157 2 2 4 1 419 420 32 +158 2 2 4 1 407 429 62 +159 2 2 4 1 397 458 20 +160 2 2 4 1 341 420 267 +161 2 2 4 1 400 417 232 +162 2 2 4 1 240 425 392 +163 2 2 4 1 267 419 265 +164 2 2 4 1 123 281 277 +165 2 2 4 1 107 446 268 +166 2 2 4 1 249 394 388 +167 2 2 4 1 265 418 266 +168 2 2 4 1 274 289 116 +169 2 2 4 1 257 407 255 +170 2 2 4 1 234 406 306 +171 2 2 4 1 308 403 233 +172 2 2 4 1 106 447 383 +173 2 2 4 1 323 426 252 +174 2 2 4 1 261 399 262 +175 2 2 4 1 136 464 271 +176 2 2 4 1 260 398 259 +177 2 2 4 1 249 444 394 +178 2 2 4 1 112 445 273 +179 2 2 4 1 273 441 112 +180 2 2 4 1 391 424 264 +181 2 2 4 1 109 449 398 +182 2 2 4 1 404 409 356 +183 2 2 4 1 267 420 419 +184 2 2 4 1 257 429 407 +185 2 2 4 1 316 452 374 +186 2 2 4 1 406 472 306 +187 2 2 4 1 308 471 403 +188 2 2 4 1 277 281 278 +189 2 2 4 1 102 425 270 +190 2 2 4 1 253 426 326 +191 2 2 4 1 404 469 409 +192 2 2 4 1 232 417 402 +193 2 2 4 1 283 284 282 +194 2 2 4 1 328 427 251 +195 2 2 4 1 386 432 107 +196 2 2 4 1 128 428 293 +197 2 2 4 1 271 430 136 +198 2 2 4 1 281 282 278 +199 2 2 4 1 121 455 272 +200 2 2 4 1 272 454 121 +201 2 2 4 1 316 456 208 +202 2 2 4 1 374 456 316 +203 2 2 4 1 231 454 272 +204 2 2 4 1 126 331 330 +205 2 2 4 1 321 435 196 +206 2 2 4 1 103 436 384 +207 2 2 4 1 196 468 321 +208 2 2 4 1 126 332 331 +209 2 2 4 1 310 439 374 +210 2 2 4 1 394 444 101 +211 2 2 4 1 327 445 112 +212 2 2 4 1 12 415 411 +213 2 2 4 1 412 417 299 +214 2 2 4 1 283 285 284 +215 2 2 4 1 232 448 400 +216 2 2 4 1 116 442 276 +217 2 2 4 1 279 443 123 +218 2 2 4 1 122 437 280 +219 2 2 4 1 283 286 285 +220 2 2 4 1 285 288 284 +221 2 2 4 1 286 287 285 +222 2 2 4 1 295 332 294 +223 2 2 4 1 294 332 126 +224 2 2 4 1 299 332 295 +225 2 2 4 1 303 311 304 +226 2 2 4 1 133 311 303 +227 2 2 4 1 306 312 311 +228 2 2 4 1 111 323 322 +229 2 2 4 1 111 322 38 +230 2 2 4 1 311 312 304 +231 2 2 4 1 303 335 133 +232 2 2 4 1 35 324 34 +233 2 2 4 1 39 326 111 +234 2 2 4 1 34 324 104 +235 2 2 4 1 133 335 93 +236 2 2 4 1 136 292 46 +237 2 2 4 1 113 328 44 +238 2 2 4 1 38 322 37 +239 2 2 4 1 42 327 112 +240 2 2 4 1 236 339 293 +241 2 2 4 1 36 317 120 +242 2 2 4 1 59 307 233 +243 2 2 4 1 87 330 329 +244 2 2 4 1 89 294 126 +245 2 2 4 1 234 305 95 +246 2 2 4 1 333 335 300 +247 2 2 4 1 307 308 233 +248 2 2 4 1 317 318 120 +249 2 2 4 1 234 306 305 +250 2 2 4 1 221 339 236 +251 2 2 4 1 19 220 18 +252 2 2 4 1 86 329 128 +253 2 2 4 1 313 314 114 +254 2 2 4 1 300 335 303 +255 2 2 4 1 325 326 40 +256 2 2 4 1 45 292 113 +257 2 2 4 1 112 325 41 +258 2 2 4 1 43 328 327 +259 2 2 4 1 54 313 114 +260 2 2 4 1 127 333 296 +261 2 2 4 1 14 134 13 +262 2 2 4 1 37 322 317 +263 2 2 4 1 296 333 300 +264 2 2 4 1 36 120 35 +265 2 2 4 1 54 114 53 +266 2 2 4 1 40 326 39 +267 2 2 4 1 39 111 38 +268 2 2 4 1 42 112 41 +269 2 2 4 1 47 136 46 +270 2 2 4 1 37 317 36 +271 2 2 4 1 45 113 44 +272 2 2 4 1 43 327 42 +273 2 2 4 1 120 324 35 +274 2 2 4 1 305 311 133 +275 2 2 4 1 78 109 77 +276 2 2 4 1 126 330 88 +277 2 2 4 1 127 294 90 +278 2 2 4 1 94 305 133 +279 2 2 4 1 336 339 221 +280 2 2 4 1 135 313 55 +281 2 2 4 1 307 309 308 +282 2 2 4 1 61 117 60 +283 2 2 4 1 117 307 60 +284 2 2 4 1 309 310 308 +285 2 2 4 1 84 116 83 +286 2 2 4 1 86 128 85 +287 2 2 4 1 87 329 86 +288 2 2 4 1 89 126 88 +289 2 2 4 1 94 133 93 +290 2 2 4 1 91 127 90 +291 2 2 4 1 59 233 58 +292 2 2 4 1 81 110 80 +293 2 2 4 1 96 234 95 +294 2 2 4 1 56 135 55 +295 2 2 4 1 93 335 92 +296 2 2 4 1 127 295 294 +297 2 2 4 1 91 333 127 +298 2 2 4 1 117 309 307 +299 2 2 4 1 46 292 45 +300 2 2 4 1 41 325 40 +301 2 2 4 1 44 328 43 +302 2 2 4 1 92 335 333 +303 2 2 4 1 55 313 54 +304 2 2 4 1 95 305 94 +305 2 2 4 1 88 330 87 +306 2 2 4 1 16 115 2 +307 2 2 4 1 90 294 89 +308 2 2 4 1 60 307 59 +309 2 2 4 1 317 319 318 +310 2 2 4 1 318 320 120 +311 2 2 4 1 306 311 305 +312 2 2 4 1 53 389 52 +313 2 2 4 1 314 315 114 +314 2 2 4 1 92 333 91 +315 2 2 4 1 118 336 221 +316 2 2 4 1 199 336 131 +317 2 2 4 1 127 296 295 +318 2 2 4 1 199 339 336 +319 2 2 4 1 319 321 318 +320 2 2 4 1 317 322 319 +321 2 2 4 1 114 389 53 +322 2 2 4 1 320 324 120 +323 2 2 4 1 314 316 315 +324 2 2 4 1 131 336 119 +325 2 2 4 1 251 337 241 +326 2 2 4 1 253 346 224 +327 2 2 4 1 296 297 295 +328 2 2 4 1 192 368 217 +329 2 2 4 1 140 350 143 +330 2 2 4 1 129 334 244 +331 2 2 4 1 322 323 319 +332 2 2 4 1 338 350 139 +333 2 2 4 1 188 368 192 +334 2 2 4 1 139 349 338 +335 2 2 4 1 224 346 340 +336 2 2 4 1 345 365 266 +337 2 2 4 1 139 350 140 +338 2 2 4 1 293 339 229 +339 2 2 4 1 229 339 199 +340 2 2 4 1 344 365 345 +341 2 2 4 1 178 340 179 +342 2 2 4 1 223 367 343 +343 2 2 4 1 168 344 202 +344 2 2 4 1 343 367 196 +345 2 2 4 1 179 340 182 +346 2 2 4 1 224 340 178 +347 2 2 4 1 132 349 139 +348 2 2 4 1 166 370 163 +349 2 2 4 1 296 298 297 +350 2 2 4 1 344 345 202 +351 2 2 4 1 337 368 241 +352 2 2 4 1 297 299 295 +353 2 2 4 1 188 192 191 +354 2 2 4 1 217 368 337 +355 2 2 4 1 176 367 223 +356 2 2 4 1 226 360 195 +357 2 2 4 1 354 358 211 +358 2 2 4 1 165 344 168 +359 2 2 4 1 195 360 158 +360 2 2 4 1 130 349 132 +361 2 2 4 1 223 224 178 +362 2 2 4 1 200 381 184 +363 2 2 4 1 210 370 166 +364 2 2 4 1 185 376 188 +365 2 2 4 1 209 347 208 +366 2 2 4 1 216 341 267 +367 2 2 4 1 157 348 214 +368 2 2 4 1 271 272 258 +369 2 2 4 1 155 348 157 +370 2 2 4 1 242 358 354 +371 2 2 4 1 223 343 252 +372 2 2 4 1 254 353 210 +373 2 2 4 1 212 273 213 +374 2 2 4 1 211 358 153 +375 2 2 4 1 154 361 198 +376 2 2 4 1 298 302 297 +377 2 2 4 1 163 370 162 +378 2 2 4 1 296 300 298 +379 2 2 4 1 353 370 210 +380 2 2 4 1 168 202 170 +381 2 2 4 1 185 187 186 +382 2 2 4 1 247 248 200 +383 2 2 4 1 141 351 138 +384 2 2 4 1 153 358 151 +385 2 2 4 1 162 225 159 +386 2 2 4 1 160 164 163 +387 2 2 4 1 165 168 167 +388 2 2 4 1 170 196 171 +389 2 2 4 1 173 176 175 +390 2 2 4 1 175 375 204 +391 2 2 4 1 176 178 177 +392 2 2 4 1 179 182 181 +393 2 2 4 1 184 381 181 +394 2 2 4 1 185 188 187 +395 2 2 4 1 186 205 200 +396 2 2 4 1 193 373 187 +397 2 2 4 1 192 209 208 +398 2 2 4 1 201 377 243 +399 2 2 4 1 238 240 207 +400 2 2 4 1 230 254 210 +401 2 2 4 1 228 259 211 +402 2 2 4 1 213 251 241 +403 2 2 4 1 214 348 342 +404 2 2 4 1 218 369 366 +405 2 2 4 1 227 250 225 +406 2 2 4 1 228 362 352 +407 2 2 4 1 235 271 258 +408 2 2 4 1 338 349 246 +409 2 2 4 1 198 362 228 +410 2 2 4 1 259 354 211 +411 2 2 4 1 227 353 264 +412 2 2 4 1 138 351 137 +413 2 2 4 1 156 361 154 +414 2 2 4 1 173 367 176 +415 2 2 4 1 147 190 149 +416 2 2 4 1 155 371 348 +417 2 2 4 1 173 174 172 +418 2 2 4 1 189 376 183 +419 2 2 4 1 194 219 218 +420 2 2 4 1 265 266 226 +421 2 2 4 1 132 137 131 +422 2 2 4 1 158 160 159 +423 2 2 4 1 182 212 189 +424 2 2 4 1 256 257 255 +425 2 2 4 1 275 359 290 +426 2 2 4 1 246 349 203 +427 2 2 4 1 351 359 275 +428 2 2 4 1 137 215 199 +429 2 2 4 1 152 371 155 +430 2 2 4 1 182 340 212 +431 2 2 4 1 201 230 210 +432 2 2 4 1 150 151 148 +433 2 2 4 1 150 152 151 +434 2 2 4 1 214 341 216 +435 2 2 4 1 226 365 360 +436 2 2 4 1 366 369 359 +437 2 2 4 1 99 206 98 +438 2 2 4 1 154 198 153 +439 2 2 4 1 173 175 174 +440 2 2 4 1 198 228 211 +441 2 2 4 1 215 274 229 +442 2 2 4 1 216 267 265 +443 2 2 4 1 242 261 219 +444 2 2 4 1 227 264 250 +445 2 2 4 1 137 199 131 +446 2 2 4 1 132 139 138 +447 2 2 4 1 140 141 138 +448 2 2 4 1 140 143 142 +449 2 2 4 1 144 147 146 +450 2 2 4 1 194 218 145 +451 2 2 4 1 151 357 148 +452 2 2 4 1 149 161 150 +453 2 2 4 1 150 371 152 +454 2 2 4 1 162 227 225 +455 2 2 4 1 167 377 201 +456 2 2 4 1 182 183 181 +457 2 2 4 1 215 229 199 +458 2 2 4 1 228 352 260 +459 2 2 4 1 300 301 298 +460 2 2 4 1 303 304 301 +461 2 2 4 1 155 156 154 +462 2 2 4 1 225 361 159 +463 2 2 4 1 175 204 174 +464 2 2 4 1 245 279 277 +465 2 2 4 1 256 263 257 +466 2 2 4 1 141 359 351 +467 2 2 4 1 144 146 145 +468 2 2 4 1 147 148 146 +469 2 2 4 1 152 153 151 +470 2 2 4 1 160 162 159 +471 2 2 4 1 168 169 167 +472 2 2 4 1 192 217 209 +473 2 2 4 1 194 242 219 +474 2 2 4 1 245 246 203 +475 2 2 4 1 204 238 207 +476 2 2 4 1 258 347 209 +477 2 2 4 1 215 275 274 +478 2 2 4 1 223 252 224 +479 2 2 4 1 250 362 225 +480 2 2 4 1 277 278 246 +481 2 2 4 1 132 138 137 +482 2 2 4 1 144 145 142 +483 2 2 4 1 218 366 145 +484 2 2 4 1 171 172 169 +485 2 2 4 1 176 177 175 +486 2 2 4 1 188 191 187 +487 2 2 4 1 213 241 189 +488 2 2 4 1 204 239 238 +489 2 2 4 1 235 258 209 +490 2 2 4 1 300 303 301 +491 2 2 4 1 139 140 138 +492 2 2 4 1 148 357 146 +493 2 2 4 1 186 373 205 +494 2 2 4 1 216 226 195 +495 2 2 4 1 290 291 289 +496 2 2 4 1 100 356 99 +497 2 2 4 1 152 154 153 +498 2 2 4 1 157 158 156 +499 2 2 4 1 170 171 169 +500 2 2 4 1 361 362 198 +501 2 2 4 1 244 245 203 +502 2 2 4 1 215 351 275 +503 2 2 4 1 266 365 226 +504 2 2 4 1 274 276 229 +505 2 2 4 1 245 277 246 +506 2 2 4 1 278 338 246 +507 2 2 4 1 275 289 274 +508 2 2 4 1 137 351 215 +509 2 2 4 1 149 150 148 +510 2 2 4 1 157 214 195 +511 2 2 4 1 171 173 172 +512 2 2 4 1 239 270 238 +513 2 2 4 1 185 186 184 +514 2 2 4 1 236 237 222 +515 2 2 4 1 242 354 261 +516 2 2 4 1 275 290 289 +517 2 2 4 1 129 244 203 +518 2 2 4 1 203 349 130 +519 2 2 4 1 146 194 145 +520 2 2 4 1 160 163 162 +521 2 2 4 1 164 365 344 +522 2 2 4 1 214 342 341 +523 2 2 4 1 130 132 131 +524 2 2 4 1 99 356 206 +525 2 2 4 1 140 142 141 +526 2 2 4 1 147 149 148 +527 2 2 4 1 190 363 149 +528 2 2 4 1 152 155 154 +529 2 2 4 1 198 211 153 +530 2 2 4 1 155 157 156 +531 2 2 4 1 158 159 156 +532 2 2 4 1 159 361 156 +533 2 2 4 1 157 195 158 +534 2 2 4 1 363 364 161 +535 2 2 4 1 165 166 163 +536 2 2 4 1 164 344 165 +537 2 2 4 1 201 210 166 +538 2 2 4 1 168 170 169 +539 2 2 4 1 204 207 174 +540 2 2 4 1 176 223 178 +541 2 2 4 1 188 376 368 +542 2 2 4 1 194 357 242 +543 2 2 4 1 214 216 195 +544 2 2 4 1 204 375 239 +545 2 2 4 1 207 380 379 +546 2 2 4 1 216 265 226 +547 2 2 4 1 261 262 219 +548 2 2 4 1 252 253 224 +549 2 2 4 1 225 362 361 +550 2 2 4 1 227 370 353 +551 2 2 4 1 268 280 254 +552 2 2 4 1 383 384 382 +553 2 2 4 1 143 144 142 +554 2 2 4 1 165 167 166 +555 2 2 4 1 201 243 230 +556 2 2 4 1 276 293 229 +557 2 2 4 1 142 366 141 +558 2 2 4 1 149 363 161 +559 2 2 4 1 262 355 219 +560 2 2 4 1 179 180 177 +561 2 2 4 1 217 337 269 +562 2 2 4 1 355 369 218 +563 2 2 4 1 230 268 254 +564 2 2 4 1 178 179 177 +565 2 2 4 1 164 165 163 +566 2 2 4 1 192 208 191 +567 2 2 4 1 146 357 194 +568 2 2 4 1 179 181 180 +569 2 2 4 1 217 269 235 +570 2 2 4 1 352 362 250 +571 2 2 4 1 169 377 167 +572 2 2 4 1 167 201 166 +573 2 2 4 1 205 256 255 +574 2 2 4 1 141 366 359 +575 2 2 4 1 145 366 142 +576 2 2 4 1 183 184 181 +577 2 2 4 1 219 355 218 +578 2 2 4 1 371 378 348 +579 2 2 4 1 161 371 150 +580 2 2 4 1 187 373 186 +581 2 2 4 1 205 255 247 +582 2 2 4 1 359 369 290 +583 2 2 4 1 130 131 119 +584 2 2 4 1 129 203 130 +585 2 2 4 1 222 232 197 +586 2 2 4 1 228 260 259 +587 2 2 4 1 158 360 160 +588 2 2 4 1 161 378 371 +589 2 2 4 1 360 365 164 +590 2 2 4 1 182 189 183 +591 2 2 4 1 183 376 185 +592 2 2 4 1 186 200 184 +593 2 2 4 1 212 213 189 +594 2 2 4 1 205 247 200 +595 2 2 4 1 217 235 209 +596 2 2 4 1 221 236 222 +597 2 2 4 1 247 249 248 +598 2 2 4 1 151 358 357 +599 2 2 4 1 191 193 187 +600 2 2 4 1 379 380 243 +601 2 2 4 1 364 378 161 +602 2 2 4 1 162 370 227 +603 2 2 4 1 196 367 171 +604 2 2 4 1 372 373 193 +605 2 2 4 1 193 374 372 +606 2 2 4 1 357 358 242 +607 2 2 4 1 207 379 174 +608 2 2 4 1 183 185 184 +609 2 2 4 1 348 378 342 +610 2 2 4 1 177 375 175 +611 2 2 4 1 377 379 243 +612 2 2 4 1 160 360 164 +613 2 2 4 1 241 376 189 +614 2 2 4 1 240 380 207 +615 2 2 4 1 364 383 382 +616 2 2 4 1 172 377 169 +617 2 2 4 1 171 367 173 +618 2 2 4 1 248 381 200 +619 2 2 4 1 181 381 180 +620 2 2 4 1 174 379 172 +621 2 2 4 1 180 375 177 +622 2 2 4 1 205 373 256 +623 2 2 4 1 256 372 263 +624 2 2 4 1 256 373 372 +625 2 2 4 1 375 385 239 +626 2 2 4 1 378 382 342 +627 2 2 4 1 172 379 377 +628 2 2 4 1 340 346 212 +629 2 2 4 1 212 346 273 +630 2 2 4 1 364 382 378 +631 2 2 4 1 221 222 197 +632 2 2 4 1 129 130 119 +633 2 2 4 1 368 376 241 +634 2 2 4 1 239 387 270 +635 2 2 4 1 248 387 385 +636 2 2 4 1 248 385 381 +637 2 2 4 1 249 387 248 +638 2 2 4 1 385 387 239 +639 2 2 4 1 180 385 375 +640 2 2 4 1 230 386 268 +641 2 2 4 1 380 386 243 +642 2 2 4 1 243 386 230 +643 2 2 4 1 381 385 180 +644 2 2 4 1 115 129 119 +645 2 2 4 1 118 221 197 +646 2 2 4 1 119 336 118 +647 2 2 4 1 387 388 270 +648 2 2 4 1 249 388 387 +649 2 2 4 1 118 197 2 +650 2 2 4 1 115 119 118 +651 2 2 4 1 220 395 18 +652 2 2 4 1 52 389 121 +653 2 2 4 1 315 389 114 +654 2 2 4 1 69 390 68 +655 2 2 4 1 50 231 49 +656 2 2 4 1 70 107 69 +657 2 2 4 1 74 391 73 +658 2 2 4 1 73 122 72 +659 2 2 4 1 75 108 74 +660 2 2 4 1 52 121 51 +661 2 2 4 1 34 104 33 +662 2 2 4 1 68 392 67 +663 2 2 4 1 25 105 24 +664 2 2 4 1 28 106 27 +665 2 2 4 1 31 103 30 +666 2 2 4 1 23 123 22 +667 2 2 4 1 115 118 2 +668 2 2 4 1 220 396 395 +669 2 2 4 1 18 395 17 +670 2 2 4 1 115 334 129 +671 2 2 4 1 107 390 69 +672 2 2 4 1 73 391 122 +673 2 2 4 1 108 391 74 +674 2 2 4 1 16 334 115 +675 2 2 4 1 11 125 10 +676 2 2 4 1 10 393 9 +677 2 2 4 1 17 334 16 +678 2 2 4 1 9 124 8 +679 2 2 4 1 245 396 279 +680 2 2 4 1 67 102 66 +681 2 2 4 1 390 392 68 +682 2 2 4 1 334 395 244 +683 2 2 4 1 64 101 63 +684 2 2 4 1 65 394 64 +685 2 2 4 1 17 395 334 +686 2 2 4 1 67 392 102 +687 2 2 4 1 395 396 244 +688 2 2 4 1 125 393 10 +689 2 2 4 1 9 393 124 +690 2 2 4 1 244 396 245 +691 2 2 4 1 64 394 101 +692 2 2 4 1 396 397 279 +693 2 2 4 1 220 397 396 +694 2 2 4 1 77 398 76 +695 2 2 4 1 190 431 288 +696 2 2 4 1 405 464 136 +697 2 2 4 1 208 467 316 +698 2 2 4 1 109 398 77 +699 2 2 4 1 112 441 325 +700 2 2 4 1 355 470 423 +701 2 2 4 1 272 464 231 +702 2 2 4 1 79 399 78 +703 2 2 4 1 345 401 202 +704 2 2 4 1 284 408 282 +705 2 2 4 1 282 408 278 +706 2 2 4 1 321 401 318 +707 2 2 4 1 345 465 401 +708 2 2 4 1 15 402 14 +709 2 2 4 1 48 405 47 +710 2 2 4 1 57 403 56 +711 2 2 4 1 98 406 97 +712 2 2 4 1 8 404 7 +713 2 2 4 1 288 457 190 +714 2 2 4 1 394 440 388 +715 2 2 4 1 78 399 109 +716 2 2 4 1 411 415 302 +717 2 2 4 1 63 407 62 +718 2 2 4 1 134 417 412 +719 2 2 4 1 314 452 316 +720 2 2 4 1 190 457 363 +721 2 2 4 1 332 400 331 +722 2 2 4 1 299 400 332 +723 2 2 4 1 423 470 110 +724 2 2 4 1 338 408 350 +725 2 2 4 1 313 410 314 +726 2 2 4 1 280 446 71 +727 2 2 4 1 12 411 11 +728 2 2 4 1 297 412 299 +729 2 2 4 1 134 415 13 +730 2 2 4 1 298 416 302 +731 2 2 4 1 320 418 324 +732 2 2 4 1 33 419 32 +733 2 2 4 1 32 420 31 +734 2 2 4 1 79 434 399 +735 2 2 4 1 374 452 310 +736 2 2 4 1 123 450 281 +737 2 2 4 1 363 422 364 +738 2 2 4 1 290 423 291 +739 2 2 4 1 355 423 369 +740 2 2 4 1 318 401 320 +741 2 2 4 1 238 425 240 +742 2 2 4 1 384 447 29 +743 2 2 4 1 401 465 320 +744 2 2 4 1 14 402 134 +745 2 2 4 1 47 405 136 +746 2 2 4 1 56 403 135 +747 2 2 4 1 206 406 98 +748 2 2 4 1 124 404 8 +749 2 2 4 1 337 427 269 +750 2 2 4 1 105 421 283 +751 2 2 4 1 236 428 237 +752 2 2 4 1 62 429 61 +753 2 2 4 1 269 430 235 +754 2 2 4 1 147 431 190 +755 2 2 4 1 143 431 144 +756 2 2 4 1 350 433 143 +757 2 2 4 1 240 432 380 +758 2 2 4 1 398 438 76 +759 2 2 4 1 101 407 63 +760 2 2 4 1 399 434 262 +761 2 2 4 1 110 434 80 +762 2 2 4 1 111 426 323 +763 2 2 4 1 323 435 319 +764 2 2 4 1 342 436 341 +765 2 2 4 1 384 436 382 +766 2 2 4 1 353 437 264 +767 2 2 4 1 280 437 254 +768 2 2 4 1 75 438 108 +769 2 2 4 1 260 438 398 +770 2 2 4 1 424 473 352 +771 2 2 4 1 102 440 66 +772 2 2 4 1 65 440 394 +773 2 2 4 1 253 441 346 +774 2 2 4 1 128 442 85 +775 2 2 4 1 84 442 116 +776 2 2 4 1 123 443 22 +777 2 2 4 1 255 444 247 +778 2 2 4 1 213 445 251 +779 2 2 4 1 70 446 107 +780 2 2 4 1 28 447 106 +781 2 2 4 1 374 439 372 +782 2 2 4 1 354 449 261 +783 2 2 4 1 105 450 24 +784 2 2 4 1 23 450 123 +785 2 2 4 1 196 435 343 +786 2 2 4 1 278 408 338 +787 2 2 4 1 356 409 206 +788 2 2 4 1 135 410 313 +789 2 2 4 1 26 421 25 +790 2 2 4 1 27 453 26 +791 2 2 4 1 272 455 258 +792 2 2 4 1 193 456 374 +793 2 2 4 1 208 456 191 +794 2 2 4 1 50 454 231 +795 2 2 4 1 121 454 51 +796 2 2 4 1 3 458 19 +797 2 2 4 1 20 458 3 +798 2 2 4 1 7 459 1 +799 2 2 4 1 6 460 96 +800 2 2 4 1 97 460 6 +801 2 2 4 1 58 461 5 +802 2 2 4 1 5 461 57 +803 2 2 4 1 2 462 15 +804 2 2 4 1 197 462 2 +805 2 2 4 1 1 459 100 +806 2 2 4 1 49 463 4 +807 2 2 4 1 4 463 48 +808 2 2 4 1 220 458 397 +809 2 2 4 1 107 432 390 +810 2 2 4 1 250 424 352 +811 2 2 4 1 389 466 121 +812 2 2 4 1 170 468 196 +813 2 2 4 1 285 457 288 +814 2 2 4 1 289 451 116 +815 2 2 4 1 108 473 424 +816 2 2 4 1 302 416 411 +817 2 2 4 1 412 415 134 +818 2 2 4 1 11 411 125 +819 2 2 4 1 302 412 297 +820 2 2 4 1 25 421 105 +821 2 2 4 1 264 424 250 +822 2 2 4 1 136 430 292 +823 2 2 4 1 125 416 414 +824 2 2 4 1 393 414 413 +825 2 2 4 1 125 414 393 +826 2 2 4 1 393 413 124 +827 2 2 4 1 312 413 304 +828 2 2 4 1 304 414 301 +829 2 2 4 1 413 414 304 +830 2 2 4 1 414 416 301 +831 2 2 4 1 13 415 12 +832 2 2 4 1 301 416 298 +833 2 2 4 1 324 418 104 +834 2 2 4 1 329 428 128 +835 2 2 4 1 222 448 232 +836 2 2 4 1 104 419 33 +837 2 2 4 1 268 446 280 +838 2 2 4 1 113 427 328 +839 2 2 4 1 31 420 103 +840 2 2 4 1 283 421 286 +841 2 2 4 1 392 425 102 +842 2 2 4 1 326 426 111 +843 2 2 4 1 299 417 400 +844 2 2 4 1 263 439 309 +845 2 2 4 1 383 447 384 +846 2 2 4 1 364 422 383 +847 2 2 4 1 369 423 290 +848 2 2 4 1 408 433 350 +849 2 2 4 1 108 424 391 +850 2 2 4 1 116 451 83 +851 2 2 4 1 270 425 238 +852 2 2 4 1 252 426 253 +853 2 2 4 1 402 417 134 +854 2 2 4 1 321 468 401 +855 2 2 4 1 251 427 337 +856 2 2 4 1 331 448 237 +857 2 2 4 1 306 472 312 +858 2 2 4 1 310 471 308 +859 2 2 4 1 293 428 236 +860 2 2 4 1 124 469 404 +861 2 2 4 1 403 471 135 +862 2 2 4 1 206 472 406 +863 2 2 4 1 61 429 117 +864 2 2 4 1 235 430 271 +865 2 2 4 1 21 443 397 +866 2 2 4 1 144 431 147 +867 2 2 4 1 82 451 291 +868 2 2 4 1 380 432 386 +869 2 2 4 1 411 416 125 +870 2 2 4 1 302 415 412 +871 2 2 4 1 288 433 284 +872 2 2 4 1 80 434 79 +873 2 2 4 1 237 448 222 +874 2 2 4 1 319 435 321 +875 2 2 4 1 83 451 82 +876 2 2 4 1 382 436 342 +877 2 2 4 1 271 464 272 +878 2 2 4 1 121 466 455 +879 2 2 4 1 401 468 202 +880 2 2 4 1 254 437 353 +881 2 2 4 1 76 438 75 +882 2 2 4 1 372 439 263 +883 2 2 4 1 284 433 408 +884 2 2 4 1 309 439 310 +885 2 2 4 1 66 440 65 +886 2 2 4 1 346 441 273 +887 2 2 4 1 85 442 84 +888 2 2 4 1 22 443 21 +889 2 2 4 1 397 443 279 +890 2 2 4 1 312 472 409 +891 2 2 4 1 410 471 310 +892 2 2 4 1 247 444 249 +893 2 2 4 1 273 445 213 +894 2 2 4 1 71 446 70 +895 2 2 4 1 29 447 28 +896 2 2 4 1 259 449 354 +897 2 2 4 1 24 450 23 +898 2 2 4 1 410 452 314 +899 2 2 4 1 291 451 289 +900 2 2 4 1 106 453 27 +901 2 2 4 1 400 448 331 +902 2 2 4 1 48 463 405 +903 2 2 4 1 15 462 402 +904 2 2 4 1 57 461 403 +905 2 2 4 1 406 460 97 +906 2 2 4 1 404 459 7 +907 2 2 4 1 258 455 347 +908 2 2 4 1 51 454 50 +909 2 2 4 1 191 456 193 +910 2 2 4 1 287 457 285 +911 2 2 4 1 438 473 108 +912 2 2 4 1 233 461 58 +913 2 2 4 1 96 460 234 +914 2 2 4 1 19 458 220 +915 2 2 4 1 232 462 197 +916 2 2 4 1 100 459 356 +917 2 2 4 1 231 463 49 +918 2 2 4 1 110 470 434 +919 2 2 4 1 320 465 418 +920 2 2 4 1 310 452 410 +921 2 2 4 1 316 467 315 +922 2 2 4 1 266 465 345 +923 2 2 4 1 347 467 208 +924 2 2 4 1 315 466 389 +925 2 2 4 1 202 468 170 +926 2 2 4 1 262 470 355 +927 2 2 4 1 352 473 260 +928 2 2 4 1 402 462 232 +929 2 2 4 1 405 463 231 +930 2 2 4 1 403 461 233 +931 2 2 4 1 234 460 406 +932 2 2 4 1 356 459 404 +933 2 2 4 1 466 467 347 +934 2 2 4 1 315 467 466 +935 2 2 4 1 409 469 312 +936 2 2 4 1 312 469 413 +937 2 2 4 1 231 464 405 +938 2 2 4 1 418 465 266 +939 2 2 4 1 409 472 206 +940 2 2 4 1 135 471 410 +941 2 2 4 1 413 469 124 +942 2 2 4 1 260 473 438 +943 2 2 4 1 434 470 262 +944 2 2 4 1 455 466 347 +$EndElements diff --git a/examples/backward_facing_step/bc_dir/bc_config b/examples/backward_facing_step/bc_dir/bc_config new file mode 100644 index 0000000..d0aa223 --- /dev/null +++ b/examples/backward_facing_step/bc_dir/bc_config @@ -0,0 +1,12 @@ +[DIRICHLET] +u = inlet -> [1.0, 0.0] + wall -> [0.0, 0.0] + +[STRESS] +u = outlet -> [0.0, 0.0] + +[NEUMANN] +k = outlet -> 0.0 + wall -> 0.0 +epsilon = outlet -> 0.0 + wall -> 0.0 diff --git a/examples/backward_facing_step/config b/examples/backward_facing_step/config new file mode 100644 index 0000000..496efed --- /dev/null +++ b/examples/backward_facing_step/config @@ -0,0 +1,46 @@ +[MESH] +filename = backward_facing_step.msh +curved_elements = False + +[FINITE ELEMENT SPACE] +elements = u -> HDiv + p -> L2 + k -> L2 + epsilon -> L2 +interpolant_order = 3 + +[DG] +DG = True +interior_penalty_coefficient = 10.0 + +[SOLVER] +linear_solver = direct +preconditioner = default +linearization_method = Oseen +nonlinear_solver = NoMixing +nonlinear_tolerance = relative -> 1e-4 + absolute -> 1e-4 +nonlinear_max_iterations = 200 +relaxation_factors = 0.5, 0.5, 0.3, 0.3 + + +[TRANSIENT] +transient = False + + +[VISUALIZATION] +save_to_file = True +save_type = .vtu +subdivision = 3 + + +[OTHER] +model = KEpsilonINS +run_dir = . +num_threads = 6 + +; Optional k-epsilon switches +wall_function = True +wall_boundary = wall +production_limiter = True +auto_turbulence_inlet = inlet diff --git a/examples/backward_facing_step/config_IC b/examples/backward_facing_step/config_IC new file mode 100644 index 0000000..abc13b9 --- /dev/null +++ b/examples/backward_facing_step/config_IC @@ -0,0 +1,26 @@ +[MESH] +filename = backward_facing_step.msh +curved_elements = False + +[FINITE ELEMENT SPACE] +elements = u -> HDiv + p -> L2 +interpolant_order = 3 + +[DG] +DG = True +interior_penalty_coefficient = 10.0 + +[SOLVER] +linear_solver = default +preconditioner = default + +[VISUALIZATION] +save_to_file = True +save_type = .sol +split_components = True + +[OTHER] +num_threads = 1 +model = Stokes +run_dir = . diff --git a/examples/backward_facing_step/ic_dir/ic_config b/examples/backward_facing_step/ic_dir/ic_config new file mode 100644 index 0000000..c738cbe --- /dev/null +++ b/examples/backward_facing_step/ic_dir/ic_config @@ -0,0 +1,6 @@ +[STOKES] +all = all -> None + +[KEpsilonINS] +u = all -> output/components_sol/u.sol +p = all -> output/components_sol/p.sol diff --git a/examples/backward_facing_step/model_dir/model_config b/examples/backward_facing_step/model_dir/model_config new file mode 100644 index 0000000..19a243c --- /dev/null +++ b/examples/backward_facing_step/model_dir/model_config @@ -0,0 +1,28 @@ +[PARAMETERS] +kinematic_viscosity = all -> 0.00002 + +; Standard high-Re k-epsilon closure constants. These are the built-in defaults, +; listed here to show what can be modified; omitting any of them changes nothing. +c_mu = all -> 0.09 +c_1 = all -> 1.44 +c_2 = all -> 1.92 +sigma_k = all -> 1.0 +sigma_epsilon = all -> 1.3 +; kappa and e_log are only read when wall_function = True. +kappa = all -> 0.4187 +e_log = all -> 9.793 + +; Numerical safeguards +k_floor = all -> 1e-8 +epsilon_floor = all -> 1e-6 +max_viscosity_ratio = all -> 200000 +production_limit_coefficient = all -> 10 +max_epsilon_k_ratio = all -> 10 + +turbulence_hydraulic_diameter = all -> 2.0 +turbulence_length_scale_ratio = all -> 0.07 + +[FUNCTIONS] +source = u -> [0.0, 0.0] + k -> 0.0 + epsilon -> 0.0 diff --git a/examples/tutorial_6/config b/examples/tutorial_6/config index 804a253..71f336f 100644 --- a/examples/tutorial_6/config +++ b/examples/tutorial_6/config @@ -23,15 +23,16 @@ nonlinear_max_iterations = 3 [TRANSIENT] transient = True scheme = implicit euler -time_range = 0.0, 1.0 +time_range = 0.0, 3.0 dt = 1e-2 [VISUALIZATION] save_to_file = True save_type = .vtu -save_frequency = 0.1, time +save_frequency = 0.01, time [OTHER] model = INS run_dir = . num_threads = 6 +resume_from_previous = False diff --git a/opencmp/helpers/dg.py b/opencmp/helpers/dg.py index edd0f5c..3ad41c6 100644 --- a/opencmp/helpers/dg.py +++ b/opencmp/helpers/dg.py @@ -81,3 +81,45 @@ def grad_avg(q: CoefficientFunction) -> CoefficientFunction: return 0.5 * (Grad(q) + Grad(q.Other())) else: return 0.5 * (Grad(q) + Grad(q).Other()) + + +def weighted_grad_avg(q: CoefficientFunction, c: CoefficientFunction) -> CoefficientFunction: + """ + Returns the average of the gradient of a field weighted by a (possibly discontinuous) coefficient. + + Args: + q: The field. + c: The coefficient weighting the gradient on each side of the facet. + + Returns: + The average of c * Grad(q) at every facet of the mesh. + """ + + # Grad must be called differently if q is a trial or testfunction instead of a coefficientfunction/gridfunction. + if isinstance(q, ProxyFunction): + return 0.5 * (c * Grad(q) + c.Other() * Grad(q.Other())) + else: + return 0.5 * (c * Grad(q) + c.Other() * Grad(q).Other()) + + +def weighted_div_avg(q: CoefficientFunction, c: CoefficientFunction) -> CoefficientFunction: + """ + Returns the average of the divergence of a field weighted by a (possibly discontinuous) coefficient. + + Args: + q: The field. + c: The coefficient weighting the divergence on each side of the facet. + + Returns: + The average of c * div(q) at every facet of the mesh. + """ + + # Grad must be called differently if q is a trial or testfunction instead of a coefficientfunction/gridfunction. + if isinstance(q, ProxyFunction): + div_q = sum(Grad(q)[i, i] for i in range(q.dim)) + div_q_other = sum(Grad(q.Other())[i, i] for i in range(q.dim)) + else: + div_q = sum(Grad(q)[i, i] for i in range(q.dim)) + div_q_other = div_q.Other() + + return 0.5 * (div_q * c + div_q_other * c.Other()) diff --git a/opencmp/helpers/limiter.py b/opencmp/helpers/limiter.py new file mode 100644 index 0000000..312f745 --- /dev/null +++ b/opencmp/helpers/limiter.py @@ -0,0 +1,298 @@ +######################################################################################################################## +# Copyright 2021 the authors (see AUTHORS file for full list). # +# # +# This file is part of OpenCMP. # +# # +# OpenCMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public # +# License as published by the Free Software Foundation, either version 2.1 of the License, or (at your option) any # +# later version. # +# # +# OpenCMP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied # +# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # +# details. # +# # +# You should have received a copy of the GNU Lesser General Public License along with OpenCMP. If not, see # +# . # +######################################################################################################################## + +""" +Bound-preserving scaling limiters for L2 DG scalar GridFunctions. Provides: + + - p1_vertex_bound : vertex-based scaling limiter (P1 Dunbar basis) + - bezier_bound : Bernstein/Bezier maximum-principle-preserving limiter -- + GUARANTEES the per-element polynomial stays in bounds + everywhere (any order). +""" + +import ngsolve as ngs +import numpy as np + + +class Limiter: + def __init__(self, mesh): + self.mesh = mesh + self._bezier_cache = {} # (id(fes), order) -> built BezierBoundLimiter + + # ── Bezier bound limiter ───────────────────────────────────────────────── + # Thin accessor over BezierBoundLimiter (defined below). The built limiter + # precomputes its change-of-basis matrix, so it is cached and reused across + # calls -- keep the Limiter instance alive to avoid rebuilding. + + def _bezier_limiter(self, fes, order): + key = (id(fes), order) + lim = self._bezier_cache.get(key) + if lim is None: + lim = BezierBoundLimiter(self.mesh, fes, order) + self._bezier_cache[key] = lim + return lim + + def bezier_bound(self, gfu, fes, order, bounds=(0.0, 1.0)): + '''Bound-preserving scaling limiter via the Bernstein/Bezier convex-hull + property. Unlike node-sampling limiters, this GUARANTEES the per-element + polynomial stays in `bounds` everywhere (no between-node leakage). + Returns number of modified elements.''' + return self._bezier_limiter(fes, order).apply(gfu, bounds) + + def ref_element_vertices_val(self, gfu: ngs.GridFunction, vertices: np.ndarray, + element_index: int, element_type: str) -> np.ndarray: + ''' + Evaluates the value of gfu at vertices of a specified element (by gfu_index). + The coefficients from gfu and Dunbar basis functions (L2) space are used to do the calculation. + ''' + if element_type == "TRIG": + [a, b, c] = gfu.vec[element_index: element_index + 3] # coefficients + val = (a - b - c) + (3 * b + c) * vertices[:, 0] + (2 * c) * vertices[:, 1] + + if element_type == "TET": + [a, b, c, d] = gfu.vec[element_index: element_index + 4] # coefficients + val = ((a - b - 2 * c - 4 * d) + (4 * b + 2 * c + 4 * d) * vertices[:, 0] + (6 * c + 4 * d) * + vertices[:, 1] + 8 * d * vertices[:, 2]) + return val + + def vertices_gfu_val(self, gfu: ngs.GridFunction, element_type: str = "TRIG") -> np.ndarray: + """ + Evaluates gfu at the vertices of every mesh cell. + """ + dof_per_element = int(len(gfu.vec) / self.mesh.ne) # dof per element + gfu_vertices_val = np.zeros((self.mesh.ne, dof_per_element)) + + if element_type == "TRIG": + vertices = np.array([(0, 0), (1, 0), (0, 1)]) + elif element_type == "TET": + vertices = np.array([(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)]) + else: + raise ValueError("Bound limiter is only implemented for element_type TRIG and TET.") + + for i in range(self.mesh.ne): # iterate over every mesh element + gfu_index = i * dof_per_element # index pointer for coefficients corresponding to the element + gfu_vertices_val[i, :] = self.ref_element_vertices_val(gfu, vertices, gfu_index, element_type) + + return gfu_vertices_val + + def p1_vertex_bound(self, gfu: ngs.GridFunction, bounds: tuple) -> None: + ''' + Vertex-based scaling limiter for P1 fields: scales the grid function within each + element if the max/min value at the mesh cell violates the upper/lower bound. + (P1 Dunbar reconstruction -- prefer bezier_bound for order >= 2.) + ''' + (r1, r2) = bounds # lower and upper bounds + if self.mesh.dim == 2: + element_type = "TRIG" + elif self.mesh.dim == 3: + element_type = "TET" + + number_of_elements = self.mesh.ne # number of mesh elements + dof_per_element = int(len(gfu.vec) / number_of_elements) # dofs per element + quad_val_gfu = self.vertices_gfu_val(gfu, element_type) # gf at the vertices + quad_min_val, quad_max_val = quad_val_gfu.min(axis=1), quad_val_gfu.max(axis=1) + theta = np.ones(number_of_elements, dtype=float) # scaling coefficients + + for i in range(number_of_elements): + nn = dof_per_element * i # index of the cell averaged value + if gfu.vec[nn] < r1: + gfu.vec[nn] = r1 + elif gfu.vec[nn] > r2: + gfu.vec[nn] = r2 + if (quad_min_val[i] < r1): + theta[i] = (gfu.vec[nn] - r1) / (gfu.vec[nn] - quad_min_val[i]) + if (quad_max_val[i] > r2): + theta2 = (gfu.vec[nn] - r2) / (gfu.vec[nn] - quad_max_val[i]) + theta[i] = min(theta[i], theta2) + if theta[i] < 1: + for k in range(1, dof_per_element): + gfu.vec[nn + k] = theta[i] * gfu.vec[nn + k] + + +# ── Reference-element evaluation helpers ───────────────────────────────────── + + +def _ndof_el(order: int, dim: int) -> int: + """DOFs per element on the reference simplex.""" + if dim == 2: + return (order + 1) * (order + 2) // 2 + return (order + 1) * (order + 2) * (order + 3) // 6 + + +def _lagrange_nodes(p: int, dim: int) -> np.ndarray: + """Uniform Lagrange nodes on the reference simplex for degree p.""" + if p == 0: + return np.ones((1, dim)) / (dim + 1) # centroid only + if dim == 2: + nodes = [(i/p, j/p) + for i in range(p + 1) + for j in range(p + 1 - i)] + else: + nodes = [(i/p, j/p, k/p) + for i in range(p + 1) + for j in range(p + 1 - i) + for k in range(p + 1 - i - j)] + return np.array(nodes, dtype=float) + + +def _build_eval_matrix_ngs(fes: ngs.FESpace, order: int, dim: int, + nodes: np.ndarray = None) -> np.ndarray: + """Build evaluation matrix via NGSolve's FiniteElement.CalcShape.""" + if nodes is None: + nodes = _lagrange_nodes(order, dim) + else: + nodes = np.asarray(nodes, dtype=float) + ndof_el = _ndof_el(order, dim) + M = np.zeros((len(nodes), ndof_el)) + fe = fes.GetFE(ngs.ElementId(ngs.VOL, 0)) + for i, node in enumerate(nodes): + # CalcShape(x, y, z) evaluates the REFERENCE shape functions directly — + # geometry-independent, so this is correct for every element regardless + # of mesh anisotropy/curvature. (Raw-coord overload; pad 2D with z=0.) + coords = (*node, 0.0) if dim == 2 else tuple(node) + M[i, :] = np.asarray(fe.CalcShape(*coords)) + return M + + +def _build_eval_matrix_gf(mesh: ngs.Mesh, fes: ngs.FESpace, order: int, dim: int, + nodes: np.ndarray = None) -> np.ndarray: + """Fallback: build evaluation matrix by probing individual basis vectors. + + Uses a centroid-ward nudge (EPS=1e-10) to keep points strictly inside + element 0, avoiding ambiguous lookups at shared vertices/faces. + """ + ndof_el = _ndof_el(order, dim) + if nodes is None: + nodes = _lagrange_nodes(order, dim) + else: + nodes = np.asarray(nodes, dtype=float) + M = np.zeros((len(nodes), ndof_el)) + + el0 = list(mesh.Elements(ngs.VOL))[0] + verts = [np.array(list(mesh[v].point)[:dim]) for v in el0.vertices] + dofs = fes.GetDofNrs(ngs.ElementId(ngs.VOL, 0)) + gf = ngs.GridFunction(fes) + gf_np = gf.vec.FV().NumPy() + ctr = np.ones(dim) / (dim + 1) # reference centroid + EPS = 1e-10 + + for j in range(ndof_el): + gf_np[:] = 0.0 + gf_np[int(dofs[j])] = 1.0 + for i, node in enumerate(nodes): + node_n = node + EPS * (ctr - node) + bary = np.concatenate(([1.0 - node_n.sum()], node_n)) + phys = sum(b * v for b, v in zip(bary, verts)) + M[i, j] = float(gf(mesh(*phys))) + return M + + +def _build_eval_matrix(mesh, fes, order, dim, nodes=None): + """Try CalcShape first, fall back to GridFunction probing.""" + try: + return _build_eval_matrix_ngs(fes, order, dim, nodes) + except Exception: + return _build_eval_matrix_gf(mesh, fes, order, dim, nodes) + + +# ── Bound-preserving Bezier/Bernstein limiter ───────────────────────────────── + +def _multiindices(p: int, n: int): + """All length-n nonneg integer tuples summing to p.""" + if n == 1: + yield (p,) + return + for i in range(p + 1): + for rest in _multiindices(p - i, n - 1): + yield (i,) + rest + + +def _bernstein_matrix(nodes: np.ndarray, p: int, dim: int) -> np.ndarray: + """Degree-p Bernstein basis values at `nodes` on the reference simplex. + B_alpha(lambda) = (p!/prod alpha_i!) * prod lambda_i^alpha_i.""" + from math import factorial + idx = list(_multiindices(p, dim + 1)) + M = np.zeros((len(nodes), len(idx))) + pf = factorial(p) + for k, node in enumerate(nodes): + lam = np.concatenate(([1.0 - node.sum()], node)) # barycentric coords + for j, a in enumerate(idx): + term = pf + for d in range(dim + 1): + term *= lam[d] ** a[d] / factorial(a[d]) + M[k, j] = term + return M + + +class BezierBoundLimiter: + """Maximum-principle-preserving scaling limiter for NGSolve DG on simplices, + any order. Scales each element's polynomial about its cell mean by a single + theta computed from the Bernstein/Bezier ordinates. Because the polynomial + lies within the convex hull of its ordinates, bounding the ordinates bounds + the polynomial EVERYWHERE — not just at sample nodes (the failure mode of + node-sampling limiters on high-order, high-curvature near-wall cells).""" + + def __init__(self, mesh: ngs.Mesh, fes: ngs.FESpace, order: int = 1): + self.mesh = mesh + self.order = order + self.dim = mesh.dim + self.ndof_el = _ndof_el(order, self.dim) + # Change-of-basis: L2 element dofs -> Bezier ordinates. Built once on the + # reference element. b = T @ c, where v = Lmat@c = Bmat@b at unisolvent + # degree-p Lagrange nodes, so T = Bmat^{-1} @ Lmat. + ref_nodes = _lagrange_nodes(order, self.dim) # ndof_el points + Lmat = _build_eval_matrix(mesh, fes, order, self.dim, nodes=ref_nodes) + Bmat = _bernstein_matrix(ref_nodes, order, self.dim) + self._to_bezier = np.linalg.solve(Bmat, Lmat) + self.dof_starts = np.array([ + fes.GetDofNrs(ngs.ElementId(ngs.VOL, i))[0] + for i in range(mesh.ne) + ], dtype=np.intp) + + def apply(self, gfu: ngs.GridFunction, bounds: tuple = (0.0, 1.0)) -> int: + r1, r2 = bounds + nd = self.ndof_el + vec = gfu.vec.FV().NumPy() + n_lim = 0 + + for base in self.dof_starts: + c = vec[base : base + nd] + ubar = c[0] # cell mean (NGSolve L2: phi0 == 1) + + # bound the average first + if ubar < r1 or ubar > r2: + ubar = float(np.clip(ubar, r1, r2)) + vec[base] = ubar + vec[base + 1 : base + nd] = 0.0 + n_lim += 1 + continue + + b = self._to_bezier @ c # Bezier ordinates: poly in [b.min, b.max] + bmin, bmax = float(b.min()), float(b.max()) + + theta = 1.0 + if bmax > r2 + 1e-14: + theta = min(theta, (r2 - ubar) / (bmax - ubar)) + if bmin < r1 - 1e-14: + theta = min(theta, (r1 - ubar) / (bmin - ubar)) + theta = max(0.0, min(1.0, theta)) + + if theta < 1.0 - 1e-14: + vec[base + 1 : base + nd] *= theta + n_lim += 1 + + return n_lim diff --git a/opencmp/helpers/wall_func.py b/opencmp/helpers/wall_func.py new file mode 100644 index 0000000..afe2416 --- /dev/null +++ b/opencmp/helpers/wall_func.py @@ -0,0 +1,393 @@ +######################################################################################################################## +# Copyright 2021 the authors (see AUTHORS file for full list). # +# # +# This file is part of OpenCMP. # +# # +# OpenCMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public # +# License as published by the Free Software Foundation, either version 2.1 of the License, or (at your option) any # +# later version. # +# # +# OpenCMP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied # +# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # +# details. # +# # +# You should have received a copy of the GNU Lesser General Public License along with OpenCMP. If not, see # +# . # +######################################################################################################################## + +"""Geometry-independent wall functions for the k-epsilon turbulence model.""" + +import logging +import numpy as np +import ngsolve as ngs + + +class KEpsilonWallFunction: + """Wall-layer eddy viscosity and dissipation for high-Re k-epsilon. + + The wall layer is found from mesh topology: cells owning a facet on + ``wall_boundary``. ``eval_nu_t`` puts those cells on the wall law and every + other cell on bulk + ``C_mu k^2/epsilon``. + + Friction velocity is either ``C_mu**0.25 * sqrt(k)`` or the square root of + the resolved tangential wall traction per unit density. Both are stored as + P0 data on wall-facet owners. ``update`` refreshes them once per Picard + iteration. + + Wall distance comes from a regularized-Eikonal solve: the continuous field + for coefficients integrated across a cell (``y_plus_field``, + ``eval_nu_wall``), and its cell average for per-cell quantities + (``y_plus_cell``, ``epsilon_wall_cell``). Turbulence production is not + handled here -- ``KEpsilonINS`` owns it. + """ + + #: y+ below which the viscous sublayer is assumed (no wall eddy viscosity). + YPLUS_VISCOUS = 11.25 + #: Recommended range for the equilibrium log-law wall treatment. + YPLUS_RECOMMENDED_MIN = 30.0 + YPLUS_RECOMMENDED_MAX = 300.0 + + def __init__(self, mesh: ngs.comp.Mesh, nu: float, C_mu: float, kappa: float, + E_log: float = 9.8, wall_boundary: str = "wall", + dist_order: int = 2, dist_relax: float = 0.1, + u_tau_method: int = 0) -> None: + self.mesh = mesh + self.nu = nu # laminar kinematic viscosity + self.C_mu = C_mu + self.kappa = kappa + self.E_log = E_log # log-law roughness constant E + self.wall_boundary = wall_boundary + if u_tau_method not in (0, 1): + raise ValueError('u_tau_method must be 0 (k-based) or 1 (velocity-based).') + self.u_tau_method = u_tau_method + self._warned_yplus_low = False + self._warned_yplus_high = False + self.h = ngs.specialcf.mesh_size + + # Piecewise-constant space shared by the masks, the cell distance and u_tau. + self._fes0 = ngs.L2(mesh, order=0) + + self._mark_wall_cells() + + self._dist_gf = self._compute_distance_field(dist_order, dist_relax) + self._dist_cell = ngs.GridFunction(self._fes0) + self._dist_cell.Set(self._dist_gf) + # y+ and epsilon_wall_cell both divide by this; guard against a Newton + # undershoot on a degenerate cell producing a negative distance. + distance = self._dist_cell.vec.FV().NumPy() + distance[:] = np.maximum(distance, 1e-12) + + # Zero until the first update(); eval_nu_t is then simply bulk everywhere. + self.u_tau_cell = ngs.GridFunction(self._fes0) + self.wall_nu_t_cell = ngs.GridFunction(self._fes0) + self.wall_shear_cell = ngs.GridFunction(self._fes0) + self._y_plus_cell_gf = ngs.GridFunction(self._fes0) + + # ------------------------------------------------------------------ + # Construction helpers + # ------------------------------------------------------------------ + + def _mark_wall_cells(self) -> None: + """Find the wall layer from mesh topology. + + Sets ``_wall_measure``, ``_marked`` / ``wall_facet_cell_mask`` and + ``mask`` to the cells that own a physical wall facet. + """ + # Assumes a simplicial mesh (one L2(0) DOF per cell). Quads/hexes parse + # but the resulting layer is untested, so refuse rather than go silently wrong. + element_types = {element.type for element in self.mesh.Elements(ngs.VOL)} + if not element_types <= {ngs.ET.TRIG, ngs.ET.TET}: + raise NotImplementedError( + 'KEpsilonWallFunction supports simplicial meshes only (TRIG in 2D, ' + f'TET in 3D); this mesh contains {sorted(t.name for t in element_types)}.') + + # An L2 basis function has no boundary trace, so plain ds() assembles to + # zero; skeleton=True evaluates it on the facet instead. + lf = ngs.LinearForm(self._fes0) + lf += self._fes0.TestFunction() * ngs.ds( + definedon=self.mesh.Boundaries(self.wall_boundary), skeleton=True) + lf.Assemble() + self._wall_measure = np.array(lf.vec).ravel() + + # Wall-facet owners: the only cells where u_tau can be evaluated directly, + # and where the algebraic epsilon condition is anchored. + self._marked = self._wall_measure > 0.0 + if not self._marked.any(): + raise ValueError( + f"KEpsilonWallFunction: boundary marker '{self.wall_boundary}' has " + f"no boundary elements. Available boundaries: " + f"{self.mesh.GetBoundaries()}.") + self.wall_facet_cell_mask = ngs.GridFunction(self._fes0) + self.wall_facet_cell_mask.vec.FV().NumPy()[:] = self._marked.astype(float) + + # A wall function is a first-cell treatment. Cells that do not own a + # physical wall facet use the bulk closure, even when facet-adjacent to + # a wall owner. + self.mask = ngs.GridFunction(self._fes0) + self.mask.vec.FV().NumPy()[:] = self._marked.astype(float) + self._wall_layer_sources = {} + self._element_dofs = { + element.nr: tuple(element.dofs) + for element in self._fes0.Elements(ngs.VOL) + } + self._wall_facets = self._find_physical_wall_facets() + + def _find_physical_wall_facets(self): + """Map each wall facet to its owner, direction and geometric distance. + + At a corner, the two facets remain separate. Their friction velocities + are combined only after projecting velocity with each facet's own normal. + """ + owner_by_vertices = {} + volume_elements = list(self.mesh.Elements(ngs.VOL)) + element_by_number = {element.nr: element for element in volume_elements} + for element in volume_elements: + for facet_id in element.facets: + numbers = tuple(sorted( + vertex.nr for vertex in self.mesh[facet_id].vertices)) + owner_by_vertices.setdefault(numbers, []).append(element.nr) + + result = [] + for facet in self.mesh.Elements(ngs.BND): + if facet.mat != self.wall_boundary: + continue + numbers = tuple(sorted(vertex.nr for vertex in facet.vertices)) + owners = owner_by_vertices.get(numbers, ()) + if len(owners) != 1: + raise ValueError( + f'Wall facet {numbers} has {len(owners)} volume owners; ' + 'expected exactly one.') + points = [np.asarray(self.mesh.vertices[number].point[:self.mesh.dim], + dtype=float) for number in numbers] + owner = element_by_number[owners[0]] + centroid = np.mean([ + np.asarray(self.mesh.vertices[vertex.nr].point[:self.mesh.dim], + dtype=float) + for vertex in owner.vertices + ], axis=0) + if self.mesh.dim == 2: + direction = points[1] - points[0] + direction /= np.linalg.norm(direction) # unit tangent + normal = np.asarray((-direction[1], direction[0])) + distance = abs(float(np.dot(centroid - points[0], normal))) + else: + direction = np.cross(points[1] - points[0], points[2] - points[0]) + direction /= np.linalg.norm(direction) # unit normal + distance = abs(float(np.dot(centroid - points[0], direction))) + result.append((owners[0], direction, max(distance, 1e-12))) + return result + + def _compute_distance_field(self, order: int, relax: float) -> ngs.GridFunction: + """Distance to ``wall_boundary`` from a regularized Eikonal solve.""" + eps = relax * self.h + fes = ngs.H1(self.mesh, order=order, dirichlet=self.wall_boundary) + u, v = fes.TnT() + y = ngs.GridFunction(fes) + + a = ngs.BilinearForm(fes) + a += ngs.grad(u) * ngs.grad(v) * ngs.dx + f = ngs.LinearForm(fes) + f += 1.0 * v * ngs.dx + a.Assemble() + f.Assemble() + y.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec + + gu = ngs.grad(u) + residual = ngs.BilinearForm(fes) + residual += ( + ngs.sqrt(gu * gu + 1e-12) * v - v + + eps * gu * ngs.grad(v) + ) * ngs.dx + ngs.solvers.Newton(residual, y, printing=False) + return y + + # ------------------------------------------------------------------ + # Per-iteration update + # ------------------------------------------------------------------ + + def _resolved_wall_shear(self, U) -> np.ndarray: + """Facet-average molecular tangential traction for each wall cell. + + The tangential projection removes pressure and the isotropic part of + the deviatoric stress. Taking the magnitude before facet integration + prevents cancellation at corners while retaining a local value instead + of imposing a streamwise/global average. + """ + n = ngs.specialcf.normal(self.mesh.dim) + grad_u = ngs.grad(U) + traction = self.nu * (grad_u + grad_u.trans) * n + tangential = traction - (traction * n) * n + + test = self._fes0.TestFunction() + form = ngs.LinearForm(self._fes0) + form += test * ngs.Norm(tangential) * ngs.ds( + definedon=self.mesh.Boundaries(self.wall_boundary), skeleton=True) + form.Assemble() + + integrated = np.asarray(form.vec).ravel() + shear = np.zeros(self._fes0.ndof) + shear[self._marked] = (integrated[self._marked] + / self._wall_measure[self._marked]) + shear = np.maximum(shear, 0.0) + scale = max(1.0, float(np.max(shear[self._marked]))) + shear[shear < 1e-14 * scale] = 0.0 + return shear + + def _wall_viscosity_from_yplus(self, y_plus: np.ndarray) -> np.ndarray: + """Log-law wall viscosity from resolved ``u_tau`` and cell distance.""" + wall_nu_t = np.zeros_like(y_plus) + active = self.mask.vec.FV().NumPy() > 0.5 + log_cells = active & (y_plus > self.YPLUS_VISCOUS) + if np.any(log_cells): + u_plus = np.log(self.E_log * y_plus[log_cells]) / self.kappa + wall_nu_t[log_cells] = self.nu * ( + y_plus[log_cells] / u_plus - 1.0) + return np.maximum(wall_nu_t, 0.0) + + def update(self, K, U=None) -> None: + """Refresh wall quantities using the configured friction-velocity method. + + Call once per Picard iteration, BEFORE ``eval_nu_t``. + """ + u_tau = np.zeros(self._fes0.ndof) + wall_nu_t = np.zeros(self._fes0.ndof) + wall_shear = np.zeros(self._fes0.ndof) + y_plus_cell = np.zeros(self._fes0.ndof) + if self.u_tau_method == 0: + projected = ngs.GridFunction(self._fes0) + projected.Set(K) + k_values = np.maximum(projected.vec.FV().NumPy(), 0.0) + u_tau[self._marked] = self.C_mu ** 0.25 * np.sqrt( + k_values[self._marked]) + y_plus_cell = (self._dist_cell.vec.FV().NumPy() + * u_tau / self.nu) + else: + if U is None: + raise ValueError('Velocity-based u_tau requires the velocity iterate.') + wall_shear = self._resolved_wall_shear(U) + u_tau[self._marked] = np.sqrt(wall_shear[self._marked]) + y_plus_cell = (self._dist_cell.vec.FV().NumPy() + * u_tau / self.nu) + wall_nu_t = self._wall_viscosity_from_yplus(y_plus_cell) + + wall_nu_t = self._wall_viscosity_from_yplus(y_plus_cell) + + self.u_tau_cell.vec.FV().NumPy()[:] = u_tau + self.wall_nu_t_cell.vec.FV().NumPy()[:] = wall_nu_t + self.wall_shear_cell.vec.FV().NumPy()[:] = wall_shear + self._y_plus_cell_gf.vec.FV().NumPy()[:] = y_plus_cell + self._warn_if_yplus_outside_recommended_range(y_plus_cell) + + def _warn_if_yplus_outside_recommended_range(self, y_plus: np.ndarray) -> None: + """Warn once for each side of the recommended wall-function band.""" + wall_values = y_plus[self._marked] + total = wall_values.size + observed_min = float(np.min(wall_values)) + observed_max = float(np.max(wall_values)) + + low_count = int(np.count_nonzero( + wall_values < self.YPLUS_RECOMMENDED_MIN)) + if low_count and not self._warned_yplus_low: + logging.warning( + 'Wall-function validity: %d/%d wall cells (%.1f%%) have y+ < %.0f; ' + 'the equilibrium log-law treatment is recommended for %.0f <= y+ <= %.0f. ' + 'Observed wall-cell range: %.6g <= y+ <= %.6g.', + low_count, total, 100.0 * low_count / total, + self.YPLUS_RECOMMENDED_MIN, self.YPLUS_RECOMMENDED_MIN, + self.YPLUS_RECOMMENDED_MAX, observed_min, observed_max) + self._warned_yplus_low = True + + high_count = int(np.count_nonzero( + wall_values > self.YPLUS_RECOMMENDED_MAX)) + if high_count and not self._warned_yplus_high: + logging.warning( + 'Wall-function validity: %d/%d wall cells (%.1f%%) have y+ > %.0f; ' + 'the equilibrium log-law treatment is recommended for %.0f <= y+ <= %.0f. ' + 'Observed wall-cell range: %.6g <= y+ <= %.6g.', + high_count, total, 100.0 * high_count / total, + self.YPLUS_RECOMMENDED_MAX, self.YPLUS_RECOMMENDED_MIN, + self.YPLUS_RECOMMENDED_MAX, observed_min, observed_max) + self._warned_yplus_high = True + + # ------------------------------------------------------------------ + # Diagnostics + # ------------------------------------------------------------------ + + def near_wall_mask(self) -> ngs.GridFunction: + """1 only on cells that own a physical wall facet.""" + return self.mask + + def wall_facet_mask(self) -> ngs.GridFunction: + """1 only on cells owning a physical wall facet.""" + return self.wall_facet_cell_mask + + def wall_distance_cell(self) -> ngs.GridFunction: + return self._dist_cell + + def wall_distance_field(self) -> ngs.GridFunction: + """Continuous Eikonal wall-distance field.""" + return self._dist_gf + + def y_plus_cell(self, K=None) -> ngs.CoefficientFunction: + """Cellwise y+, from the cell-averaged wall distance.""" + if K is not None and self.u_tau_method == 0: + self.update(K) + return self._y_plus_cell_gf + + def y_plus_field(self, K=None) -> ngs.CoefficientFunction: + """Pointwise y+, from the continuous Eikonal distance. + + Use this wherever y+ feeds a coefficient the weak form integrates over a + volume, so the result varies across the wall cell instead of being one + number per cell. + """ + if K is not None and self.u_tau_method == 0: + self.update(K) + return self._dist_gf * self.u_tau_cell / self.nu + + def epsilon_wall_cell(self, K) -> ngs.CoefficientFunction: + """Wall dissipation selected from one cell-averaged y+ value.""" + k_nonnegative = ngs.IfPos(K, K, 0.0) + log_layer = (self.C_mu ** 0.75 * k_nonnegative ** 1.5 + / (self.kappa * self._dist_cell)) + viscous_sublayer = (2.0 * self.nu * k_nonnegative + / self._dist_cell ** 2) + return ngs.IfPos(self.y_plus_cell(K) - self.YPLUS_VISCOUS, + log_layer, viscous_sublayer) + + # ------------------------------------------------------------------ + # Eddy viscosity + # ------------------------------------------------------------------ + + def _wall_law(self, yplus, K, epsilon) -> ngs.CoefficientFunction: + """Wall-only eddy viscosity as a function of the y+ handed in:: + + nu_t = 0 y+ < 11.25 (sublayer) + nu_t = nu*(y+ / ((1/kappa)*ln(E*y+)) - 1) y+ >= 11.25 (log law) + + Wall-owner cells never fall back to the bulk k-epsilon viscosity. + """ + # ln(E*y+) vanishes at y+ = 1/E; pinning y+ to YPLUS_VISCOUS below the + # sublayer threshold avoids that root (the clamp below then zeroes it). + yplus_log = ngs.IfPos(yplus - self.YPLUS_VISCOUS, yplus, self.YPLUS_VISCOUS) + u_plus = ngs.log(self.E_log * yplus_log) / self.kappa + log_law = self.nu * (yplus_log / u_plus - 1.0) + log_law = ngs.IfPos(log_law, log_law, 0.0) + + return ngs.IfPos(yplus - self.YPLUS_VISCOUS, log_law, 0.0) + + def eval_nu_wall(self, K, epsilon) -> ngs.CoefficientFunction: + """One wall eddy viscosity branch per wall cell, selected by P0 y+.""" + return self._wall_law(self.y_plus_cell(K), K, epsilon) + + def eval_nu_t(self, K, E) -> ngs.CoefficientFunction: + """Wall law inside the mask, bulk ``C_mu k^2/epsilon`` outside it. + + Requires :meth:`update` to have run this iteration. + """ + # Cell-based, not pointwise: the Eikonal distance is zero ON the wall, so + # a pointwise wall law collapses there and jumps to bulk one cell in. + nu_t_bulk = self.C_mu * K ** 2 / E + nu_t = self.mask * self.eval_nu_wall(K, E) + (1 - self.mask) * nu_t_bulk + return ngs.IfPos(nu_t, nu_t, 0) # turbulent viscosity cannot be negative diff --git a/opencmp/models/__init__.py b/opencmp/models/__init__.py index 63b8790..f3999b2 100644 --- a/opencmp/models/__init__.py +++ b/opencmp/models/__init__.py @@ -26,6 +26,8 @@ from .stokes import Stokes from .stokes_dim import StokesDIM from .multi_component_ins import MultiComponentINS +from .k_epsilon import KEpsilonINS + models_dict = {"INS": INS, "INS-DIM": INSDIM, @@ -33,7 +35,8 @@ "Poisson-DIM": PoissonDIM, "Stokes": Stokes, "Stokes-DIM": StokesDIM, - "MultiComponentINS": MultiComponentINS} + "MultiComponentINS": MultiComponentINS, + "KEpsilonINS": KEpsilonINS} # Helper functions from .misc import get_model_class diff --git a/opencmp/models/ins.py b/opencmp/models/ins.py index 602b4d5..0c01469 100644 --- a/opencmp/models/ins.py +++ b/opencmp/models/ins.py @@ -25,7 +25,7 @@ Preconditioner, div, dx from ..helpers.ngsolve_ import get_special_functions -from ..helpers.dg import avg, jump, grad_avg +from ..helpers.dg import avg, jump, weighted_grad_avg from . import Model from ..helpers.error import norm, mean @@ -92,6 +92,14 @@ def _set_model_parameters(self) -> None: # Remove the source term for the conservation of momentum if it's not being solved. self.f.pop('u') + def _get_effective_viscosity(self, time_step: int): + """ + Return the viscosity used by the momentum weak form. + + Molecular viscosity here; turbulence models add their eddy viscosity. + """ + return self.kv[time_step] + def _construct_fes(self) -> FESpace: return FESpace(self._construct_fes_helper(), dgjumps=self.DG) @@ -179,6 +187,7 @@ def construct_bilinear_time_ODE(self, U: Union[List[ProxyFunction], List[GridFun dt: Parameter = Parameter(1.0), time_step: int = 0) -> List[BilinearForm]: w = self._get_wind(U, time_step) + kv = self._get_effective_viscosity(time_step) # Define the special DG functions n, _, alpha, I_mat = get_special_functions(self.mesh, self.nu) @@ -188,7 +197,7 @@ def construct_bilinear_time_ODE(self, U: Union[List[ProxyFunction], List[GridFun v = V[self.model_components['u']] # Domain integrals. Newtonian Stress - a = dt * (self.kv[time_step] * InnerProduct(Grad(u), Grad(v))) * dx + a = dt * (kv * InnerProduct(Grad(u), Grad(v))) * dx if self.linearize == 'Oseen': # Linearized convection term. @@ -198,9 +207,9 @@ def construct_bilinear_time_ODE(self, U: Union[List[ProxyFunction], List[GridFun # Penalty for dirichlet BCs if self.dirichlet_names.get('u', None) is not None: a += dt * ( - self.kv[time_step] * alpha * u * v # 1/2 of penalty term for u=g on 𝚪_D from ∇u^ - - self.kv[time_step] * InnerProduct(Grad(u), OuterProduct(v, n)) # ∇u^ = ∇u - - self.kv[time_step] * InnerProduct(Grad(v), OuterProduct(u, n)) # 1/2 of penalty for u=g on 𝚪_D + kv * alpha * u * v # 1/2 of penalty term for u=g on 𝚪_D from ∇u^ + - kv * InnerProduct(Grad(u), OuterProduct(v, n)) # ∇u^ = ∇u + - kv * InnerProduct(Grad(v), OuterProduct(u, n)) # 1/2 of penalty for u=g on 𝚪_D ) * self._ds(self.dirichlet_names['u']) if self.linearize == 'Oseen': @@ -228,6 +237,7 @@ def construct_bilinear_time_coefficient(self, U: List[ProxyFunction], V: List[Pr time_step: int) -> List[BilinearForm]: w = self._get_wind(U, time_step) + kv = self._get_effective_viscosity(time_step) # Define the special DG functions. n, _, alpha, I_mat = get_special_functions(self.mesh, self.nu) @@ -246,10 +256,12 @@ def construct_bilinear_time_coefficient(self, U: List[ProxyFunction], V: List[Pr if self.DG: avg_u = avg(u) jump_u = jump(u) - avg_grad_u = grad_avg(u) + facet_kv = ngsolve.CoefficientFunction(kv) + avg_kv = avg(facet_kv) + avg_grad_u = weighted_grad_avg(u, facet_kv) jump_v = jump(v) - avg_grad_v = grad_avg(v) + avg_grad_v = weighted_grad_avg(v, facet_kv) # Penalty for discontinuities # TODO: Why are these in here? @@ -258,9 +270,9 @@ def construct_bilinear_time_coefficient(self, U: List[ProxyFunction], V: List[Pr # E.g. p is an unknown, grad(p) would be fixed in terms of a previous value of p and could NOT be # solved for. a += dt * ( - self.kv[time_step] * alpha * InnerProduct(jump_u, jump_v) # Penalty term for u+=u- on 𝚪_I from ∇u^ - - self.kv[time_step] * InnerProduct(avg_grad_u, OuterProduct(jump_v, n)) # Stress - - self.kv[time_step] * InnerProduct(avg_grad_v, OuterProduct(jump_u, n)) # U + avg_kv * alpha * InnerProduct(jump_u, jump_v) # Penalty term for u+=u- on 𝚪_I from ∇u^ + - InnerProduct(avg_grad_u, OuterProduct(jump_v, n)) # Stress + - InnerProduct(avg_grad_v, OuterProduct(jump_u, n)) # U ) * dx(skeleton=True) if self.linearize == 'Oseen': @@ -273,6 +285,7 @@ def construct_linear(self, V: List[ProxyFunction], gfu_0: Optional[List[GridFunc dt: Parameter, time_step: int) -> List[LinearForm]: w = self._get_wind(gfu_0, time_step) + kv = self._get_effective_viscosity(time_step) # Define the special DG functions. n, h, alpha, I_mat = get_special_functions(self.mesh, self.nu) @@ -288,8 +301,8 @@ def construct_linear(self, V: List[ProxyFunction], gfu_0: Optional[List[GridFunc for marker in self.BC.get('dirichlet', {}).get('u', {}): g = self.BC['dirichlet']['u'][marker][time_step] L += dt * ( - self.kv[time_step] * alpha * g * v # 1/2 of penalty for u=g from ∇u^ on 𝚪_D - - self.kv[time_step] * InnerProduct(Grad(v), OuterProduct(g, n)) # 1/2 of penalty for u=g + kv * alpha * g * v # 1/2 of penalty for u=g from ∇u^ on 𝚪_D + - kv * InnerProduct(Grad(v), OuterProduct(g, n)) # 1/2 of penalty for u=g ) * self._ds(marker) if self.linearize == 'Oseen': diff --git a/opencmp/models/k_epsilon.py b/opencmp/models/k_epsilon.py new file mode 100644 index 0000000..3ebdcf5 --- /dev/null +++ b/opencmp/models/k_epsilon.py @@ -0,0 +1,666 @@ +######################################################################################################################## +# Copyright 2021 the authors (see AUTHORS file for full list). # +# # +# This file is part of OpenCMP. # +# # +# OpenCMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public # +# License as published by the Free Software Foundation, either version 2.1 of the License, or (at your option) any # +# later version. # +# # +# OpenCMP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied # +# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # +# details. # +# # +# You should have received a copy of the GNU Lesser General Public License along with OpenCMP. If not, see # +# . # +######################################################################################################################## + +"""Single-phase incompressible Navier--Stokes with a k-epsilon closure.""" + +import logging +from typing import Dict, List, Optional, Union + +import numpy as np +import ngsolve as ngs +from ngsolve import (BilinearForm, FESpace, GridFunction, LinearForm, + Parameter, Preconditioner) +from ngsolve.comp import ProxyFunction + +from .ins import INS +from ..helpers.dg import avg, jump, weighted_grad_avg +from ..helpers.limiter import Limiter +from ..helpers.math import Max +from ..helpers.ngsolve_ import get_special_functions +from ..helpers.wall_func import KEpsilonWallFunction + + +class KEpsilonINS(INS): + """RANS INS model closed with the standard high-Reynolds-number k-epsilon equations. + + Adds transport equations for turbulent kinetic energy and dissipation to INS. + """ + + DEFAULT_PARAMETERS = { + 'c_mu': 0.09, + 'c_1': 1.44, + 'c_2': 1.92, + 'sigma_k': 1.0, + 'sigma_epsilon': 1.3, + 'kappa': 0.4187, + 'e_log': 9.793, + 'k_floor': 1e-10, + 'epsilon_floor': 1e-10, + 'max_viscosity_ratio': 1e5, + 'production_limit_coefficient': 10.0, + 'max_epsilon_k_ratio': 10.0, + 'realizability_coefficient': 1.0, + 'wall_u_tau_method': 0.0, + 'turbulence_hydraulic_diameter': 0.0, + 'turbulence_length_scale_ratio': 0.07, + } + + #: Optional ``[OTHER]`` switches + DEFAULT_OPTIONS = { + 'wall_function': True, + 'wall_boundary': 'wall', + 'production_limiter': True, + 'realizability_limiter': True, + # Boundary marker whose prescribed velocity generates default inlet and + # initial k-epsilon data. Empty means that all values remain user supplied. + 'auto_turbulence_inlet': '', + } + + def _parameter(self, parameters: Dict, name: str) -> float: + """Config value if present, else the standard constant from DEFAULT_PARAMETERS. + + These are all constants, so the per-time-level list the config parser + returns is collapsed to its first entry. + """ + if name not in parameters: + return self.DEFAULT_PARAMETERS[name] + return parameters[name]['all'][0] + + def _optional_config(self, key: str): + """One optional [OTHER] switch, from DEFAULT_OPTIONS if the config omits it.""" + default = self.DEFAULT_OPTIONS[key] + try: + value = self.config.get_item(['OTHER', key], type(default), quiet=True) + except Exception: + return default + return default if value is None else value + + def _pre_init(self) -> None: + for key in self.DEFAULT_OPTIONS: + setattr(self, key, self._optional_config(key)) + + def _define_model_components(self) -> Dict[str, Optional[int]]: + return {'u': 0, 'p': 1, 'k': 2, 'epsilon': 3} + + def _define_model_local_error_components(self) -> Dict[str, bool]: + return {'u': True, 'p': False, 'k': True, 'epsilon': True} + + def _define_time_derivative_components(self) -> List[Dict[str, bool]]: + return [{'u': True, 'p': False, 'k': True, 'epsilon': True}] + + def _define_bc_types(self) -> List[str]: + return super()._define_bc_types() + ['neumann'] + + def _construct_fes(self) -> FESpace: + spaces = self._construct_fes_helper() + scalar_order = max(self.interp_ord - 1, 0) + + for component in ('k', 'epsilon'): + element = self.element[component] + kwargs = { + 'mesh': self.mesh, + 'order': scalar_order, + 'dgjumps': self.DG, + } + if element != 'L2': + kwargs['dirichlet'] = self.dirichlet_names.get(component, '') + spaces.append(getattr(ngs, element)(**kwargs)) + + return FESpace(spaces, dgjumps=self.DG) + + + def _bound_turbulence(self, gfu: GridFunction) -> None: + """Floor and slope-limit the stored k and epsilon. + + Not optional: an unfloored k or epsilon zeroes the turbulent viscosity and + with it the production term, an absorbing state with no way back above + zero. Bezier bounds the DG polynomial everywhere (not just at sample + nodes), so a positive cell mean can't hide a negative value inside the cell. + """ + if self._limiter is None: + return + comp = self.model_components + for component, floor in (('k', self.k_floor), + ('epsilon', self.epsilon_floor)): + field = gfu.components[comp[component]] + # Pass the stable FES, not the per-call component, so the limiter's + # cache keys match across iterations. + space = self.fes.components[comp[component]] + self._limiter.bezier_bound(field, space, space.globalorder, + (floor, 1e20)) + + def _set_model_parameters(self) -> None: + super()._set_model_parameters() + parameters = self.model_functions.model_parameters_dict + self.C_mu = self._parameter(parameters, 'c_mu') + self.C_1 = self._parameter(parameters, 'c_1') + self.C_2 = self._parameter(parameters, 'c_2') + self.sigma_k = self._parameter(parameters, 'sigma_k') + self.sigma_epsilon = self._parameter(parameters, 'sigma_epsilon') + self.kappa = self._parameter(parameters, 'kappa') + self.E_log = self._parameter(parameters, 'e_log') + + # Floors and the viscosity ratio cap are numerical safeguards, not + # replacements for physically meaningful ICs and BCs. + self.k_floor = self._parameter(parameters, 'k_floor') + self.epsilon_floor = self._parameter(parameters, 'epsilon_floor') + self.max_viscosity_ratio = self._parameter(parameters, 'max_viscosity_ratio') + self.production_limit_coefficient = self._parameter( + parameters, 'production_limit_coefficient') + self.max_epsilon_k_ratio = self._parameter(parameters, 'max_epsilon_k_ratio') + self.realizability_coefficient = self._parameter( + parameters, 'realizability_coefficient') + self.wall_u_tau_method = int(self._parameter( + parameters, 'wall_u_tau_method')) + self.turbulence_hydraulic_diameter = self._parameter( + parameters, 'turbulence_hydraulic_diameter') + self.turbulence_length_scale_ratio = self._parameter( + parameters, 'turbulence_length_scale_ratio') + + if self.production_limit_coefficient <= 0.0: + raise ValueError('production_limit_coefficient must be positive.') + if self.max_epsilon_k_ratio <= 0.0: + raise ValueError('max_epsilon_k_ratio must be positive.') + if self.wall_u_tau_method not in (0, 1): + raise ValueError( + 'wall_u_tau_method must be 0 (k-based) or 1 (velocity-based).') + if self.auto_turbulence_inlet: + if self.turbulence_hydraulic_diameter <= 0.0: + raise ValueError( + 'turbulence_hydraulic_diameter must be positive when ' + 'auto_turbulence_inlet is enabled.') + if self.turbulence_length_scale_ratio <= 0.0: + raise ValueError( + 'turbulence_length_scale_ratio must be positive when ' + 'auto_turbulence_inlet is enabled.') + + def _post_init(self) -> None: + super()._post_init() + if self.linearize != 'Oseen': + raise NotImplementedError('KEpsilonINS currently supports Oseen linearization only.') + + if self.auto_turbulence_inlet: + self._apply_auto_turbulence_defaults() + + self.UIter = ngs.GridFunction(self.fes) + self.UIter.vec.data = self.IC.vec + + try: + relaxation = self.config.get_list(['SOLVER', 'relaxation_factors'], float) + except Exception: + relaxation = [] + self.relaxation_factors = ( + relaxation if len(relaxation) == len(self.model_components) + else [1.0] * len(self.model_components) + ) + # Bounding needs a discontinuous L2 space; with CG/H1 turbulence spaces + # the coefficient floors in _regularized_turbulence are the only safeguard. + self._bounded = (self.DG and self.element['k'] == 'L2' + and self.element['epsilon'] == 'L2') + self._limiter = Limiter(self.mesh) if self._bounded else None + # The IC is user-supplied and may sit below the floors (e.g. zero fields). + self._bound_turbulence(self.UIter) + # nu_t is built straight from the DG solution, no H1 recovery. These are + # live views into UIter, so they track it with no update step. + self._k_cell = self.UIter.components[self.model_components['k']] + self._epsilon_cell = self.UIter.components[self.model_components['epsilon']] + self._wallf = None + if self.wall_function: + if not self.DG or self.element['epsilon'] != 'L2': + raise NotImplementedError( + 'The wall-law epsilon Dirichlet condition requires a ' + 'discontinuous L2 epsilon space.') + self._wallf = KEpsilonWallFunction( + self.mesh, + nu=self.kv[0], + C_mu=self.C_mu, + kappa=self.kappa, + E_log=self.E_log, + wall_boundary=self.wall_boundary, + u_tau_method=self.wall_u_tau_method, + ) + self._update_wall_function() + + # Build and compile nu_t once per time level so every form reuses the + # same optimized evaluation tree. + self._turbulent_viscosity = [] + for time_step in range(len(self.t_param)): + nu_t = self._build_turbulent_viscosity(time_step) + self._turbulent_viscosity.append(nu_t.Compile()) + + self._normal, _, self._penalty, _ = get_special_functions(self.mesh, self.nu) + + def _apply_auto_turbulence_defaults(self) -> None: + """Fill missing inlet and initial k-epsilon data from inlet velocity. + + Explicit user values always take precedence. This first implementation + supports one steady inlet marker. + """ + marker = self.auto_turbulence_inlet + velocity_markers = self.BC.get('dirichlet', {}).get('u', {}) + if marker not in velocity_markers: + raise ValueError( + f"auto_turbulence_inlet '{marker}' has no prescribed Dirichlet " + 'velocity boundary condition.') + + # TODO: Support multiple inlet markers and form the uniform IC from + # flow-rate-weighted inlet turbulence quantities. + # TODO: Re-evaluate automatic inlet values for time-dependent velocity + # data; the present implementation intentionally uses only t=0. + inlet_velocity = velocity_markers[marker][0] + boundary = self.mesh.Boundaries(marker) + area = float(ngs.Integrate(1.0, self.mesh, definedon=boundary)) + if area <= 0.0: + raise ValueError( + f"auto_turbulence_inlet '{marker}' has zero boundary measure.") + + normal = ngs.specialcf.normal(self.mesh.dim) + outward_flux = float(ngs.Integrate( + inlet_velocity * normal, self.mesh, definedon=boundary)) + flux_tolerance = 1e-14 * max(1.0, area) + if outward_flux >= -flux_tolerance: + raise ValueError( + f"auto_turbulence_inlet '{marker}' must have nonzero inward net " + f'flux; computed outward flux is {outward_flux:.6g}.') + + bulk_velocity = -outward_flux / area + reynolds = (bulk_velocity * self.turbulence_hydraulic_diameter + / self.kv[0]) + intensity = 0.16 * reynolds ** (-1.0 / 8.0) + k_value = 1.5 * (bulk_velocity * intensity) ** 2 + length_scale = (self.turbulence_length_scale_ratio + * self.turbulence_hydraulic_diameter) + epsilon_value = (self.C_mu ** 0.75 * k_value ** 1.5 + / length_scale) + + dirichlet = self.BC.setdefault('dirichlet', {}) + for component, value in (('k', k_value), ('epsilon', epsilon_value)): + component_boundaries = dirichlet.setdefault(component, {}) + if marker not in component_boundaries: + component_boundaries[marker] = [value] * len(self.t_param) + existing = self.dirichlet_names.get(component, '') + names = [name for name in existing.split('|') if name] + if marker not in names: + names.append(marker) + self.dirichlet_names[component] = '|'.join(names) + + ic_dict = self.ic_functions.ic_dict.get(self.name(), {}) + explicit_all = 'all' in ic_dict + if not explicit_all and 'k' not in ic_dict: + self.IC.components[self.model_components_ic['k']].Set(k_value) + if not explicit_all and 'epsilon' not in ic_dict: + self.IC.components[self.model_components_ic['epsilon']].Set( + epsilon_value) + + def _regularized_turbulence(self, time_step: int): + """Floor-bounded k and epsilon, for use in ratios like k**2/epsilon. + + Guards against division by zero and negative coefficients during + nonlinear iterations. Does not modify the stored ``UIter`` fields. + """ + comp = self.model_components + k_safe = Max(self._k_cell, ngs.CoefficientFunction(self.k_floor)) + epsilon_safe = Max( + self._epsilon_cell, ngs.CoefficientFunction(self.epsilon_floor)) + # epsilon >= C_mu k**2 / (max_viscosity_ratio * nu), tied to local k so + # C_mu k**2/epsilon stays under the ratio by construction. A scalar + # limiter bound can't express this, so it's applied here instead. + epsilon_safe = Max( + epsilon_safe, + self.C_mu * k_safe ** 2 + / (self.max_viscosity_ratio * self.kv[time_step])) + return k_safe, epsilon_safe + + def _update_wall_function(self) -> None: + """Refresh wall data from the current lagged turbulence and velocity.""" + comp = self.model_components + self._wallf.update( + self.UIter.components[comp['k']], + self.UIter.components[comp['u']]) + + def _build_turbulent_viscosity(self, time_step: int): + k, epsilon = self._regularized_turbulence(time_step) + comp = self.model_components + k_raw = self.UIter.components[comp['k']] + epsilon_raw = self.UIter.components[comp['epsilon']] + if self._bounded: + # An epsilon at its floor is a clamped undershoot, not a physical + # state; trusting k**2/epsilon there would inflate nu_t to the cap. + # Fade it out smoothly instead (see _epsilon_trust). + bulk_valid = self._epsilon_trust() + else: + bulk_valid = ngs.IfPos( + k_raw, ngs.IfPos(epsilon_raw, 1.0, 0.0), 0.0) + + if self._wallf is not None: + nu_t = self._wallf.eval_nu_t(k, epsilon) + wall_mask = self._wallf.near_wall_mask() + # Wall law stays active throughout the wall layer; in the bulk, + # bulk_valid switches turbulence off instead of trusting the floor. + nu_t = wall_mask * nu_t + (1.0 - wall_mask) * bulk_valid * nu_t + else: + nu_t = bulk_valid * self.C_mu * k ** 2 / epsilon + + if self.realizability_limiter: + nu_t = self._realizable_nu_t(nu_t, k) + + cap = self.max_viscosity_ratio * self.kv[time_step] + return ngs.IfPos(nu_t, ngs.IfPos(cap - nu_t, nu_t, cap), 0.0) + + def _strain_magnitude(self): + """S = sqrt(2 S_ij S_ij) from the lagged velocity.""" + velocity = self.UIter.components[self.model_components['u']] + strain = 0.5 * (ngs.grad(velocity) + ngs.grad(velocity).trans) + return ngs.sqrt(2.0 * ngs.InnerProduct(strain, strain) + 1e-30) + + def _realizable_nu_t(self, nu_t, k): + """Durbin (1996) realizability bound on the turbulent time scale. + + T = min(k/epsilon, a / (sqrt(6) C_mu S)) with nu_t = C_mu k T, i.e. + nu_t <= a k / (sqrt(6) S). Follows from positivity of the normal Reynolds + stresses, so it caps nu_t by the local strain instead of letting k**2 + over a collapsing epsilon run to the viscosity ratio. + """ + bound = (self.realizability_coefficient * k + / (np.sqrt(6.0) * self._strain_magnitude())) + return ngs.IfPos(bound - nu_t, nu_t, bound) + + def _get_turbulent_viscosity(self, time_step: int): + return self._turbulent_viscosity[time_step] + + def _get_effective_viscosity(self, time_step: int): + return self.kv[time_step] + self._get_turbulent_viscosity(time_step) + + def _epsilon_trust(self): + """Smooth 0..1 weight that vanishes as epsilon nears its floor (a clamped + undershoot there, not a physical state). Healthy cells get weight ~1. + ponytail: 100 is a trust margin, not physics. + """ + epsilon_raw = self._epsilon_cell + return epsilon_raw / (epsilon_raw + 100.0 * self.epsilon_floor) + + def _epsilon_k_ratio(self, k, epsilon): + """Smooth, upper-bounded epsilon/k reaction coefficient. + + Raw epsilon/k diverges as k -> 0; this form tends to max_epsilon_k_ratio + instead. The trust weight also removes the artificial O(1) sink that + forms where both k and epsilon sit at their floors. + """ + ratio = epsilon / (k + epsilon / self.max_epsilon_k_ratio) + if self._bounded: + ratio = ratio * self._epsilon_trust() + return ratio + + def _limit_production(self, production, epsilon): + """Apply the configured bulk production cap uniformly in every cell.""" + limit = self.production_limit_coefficient * epsilon + return ngs.IfPos(limit - production, production, limit) + + def _production(self, time_step: int): + velocity = self.UIter.components[self.model_components['u']] + nu_t = self._get_turbulent_viscosity(time_step) + strain = 0.5 * (ngs.grad(velocity) + ngs.grad(velocity).trans) + production = 2.0 * nu_t * ngs.InnerProduct(strain, strain) + production = ngs.IfPos(production, production, 0.0) + + if self.production_limiter: + _, epsilon = self._regularized_turbulence(time_step) + production = self._limit_production(production, epsilon) + + return production + + def _neumann_markers(self, component: str) -> str: + markers = self.BC.get('neumann', {}).get(component, {}).keys() + if component == 'epsilon' and self._wallf is not None: + markers = (marker for marker in markers + if marker != self.wall_boundary) + return '|'.join(markers) + + def _epsilon_dirichlet_markers(self) -> str: + """Configured epsilon boundaries plus the active wall-function wall.""" + markers = list(self.BC.get('dirichlet', {}).get('epsilon', {}).keys()) + if self._wallf is not None and self.wall_boundary not in markers: + markers.append(self.wall_boundary) + return '|'.join(markers) + + def _add_scalar_transport(self, form, scalar, test, wind, diffusivity, + dirichlet_markers: str, neumann_markers: str, + dt: Parameter): + n = self._normal + form += -dt * scalar * (wind * ngs.grad(test)) * ngs.dx + form += dt * diffusivity * ngs.grad(scalar) * ngs.grad(test) * ngs.dx + + if self.DG: + wind_n = wind * n + flux = avg(scalar) * wind_n + 0.5 * ngs.Norm(wind_n) * jump(scalar) + form += dt * jump(test) * flux * ngs.dx(skeleton=True) + facet_diffusivity = ngs.CoefficientFunction(diffusivity) + avg_diffusivity = avg(facet_diffusivity) + avg_diffusive_grad_test = weighted_grad_avg( + test, facet_diffusivity) + avg_diffusive_grad_scalar = weighted_grad_avg( + scalar, facet_diffusivity) + form += -dt * (n * avg_diffusive_grad_test) * jump( + scalar) * ngs.dx(skeleton=True) + form += dt * ( + avg_diffusivity * self._penalty * jump(scalar) + - avg_diffusive_grad_scalar * n + ) * jump(test) * ngs.dx(skeleton=True) + + if dirichlet_markers: + form += dt * test * ( + 0.5 * scalar * wind_n + 0.5 * scalar * ngs.Norm(wind_n) + ) * self._ds(dirichlet_markers) + form += -dt * diffusivity * scalar * ( + ngs.grad(test) * n) * self._ds(dirichlet_markers) + form += dt * diffusivity * ( + self._penalty * scalar - ngs.grad(scalar) * n + ) * test * self._ds(dirichlet_markers) + + if neumann_markers: + form += dt * test * scalar * Max( + wind_n, ngs.CoefficientFunction(0.0) + ) * self._ds(neumann_markers) + + return form + + def construct_bilinear_time_ODE( + self, U: Union[List[ProxyFunction], List[GridFunction]], + V: List[ProxyFunction], dt: Parameter = Parameter(1.0), + time_step: int = 0) -> List[BilinearForm]: + forms = super().construct_bilinear_time_ODE(U, V, dt, time_step) + comp = self.model_components + wind = self._get_wind(U, time_step) + nu_t = self._get_turbulent_viscosity(time_step) + k = U[comp['k']] + epsilon = U[comp['epsilon']] + zeta = V[comp['k']] + psi = V[comp['epsilon']] + + d_k = self.kv[time_step] + nu_t / self.sigma_k + d_epsilon = self.kv[time_step] + nu_t / self.sigma_epsilon + k_previous, epsilon_previous = self._regularized_turbulence(time_step) + form = forms[0] + form = self._add_scalar_transport( + form, k, zeta, wind, d_k, + self.dirichlet_names.get('k', ''), self._neumann_markers('k'), dt) + form = self._add_scalar_transport( + form, epsilon, psi, wind, d_epsilon, + self._epsilon_dirichlet_markers(), + self._neumann_markers('epsilon'), dt) + # Sink terms go on the implicit side with lagged coefficients only + # (standard segregated k-epsilon linearization); on the RHS they make + # the Picard iteration unstable when k or epsilon changes rapidly. + epsilon_k_ratio = self._epsilon_k_ratio(k_previous, epsilon_previous) + form += dt * epsilon_k_ratio * k * zeta * ngs.dx + form += dt * self.C_2 * ( + epsilon_k_ratio) * epsilon * psi * ngs.dx + forms[0] = form + return forms + + def construct_linear(self, V: List[ProxyFunction], + gfu_0: Optional[List[GridFunction]], dt: Parameter, + time_step: int) -> List[LinearForm]: + forms = super().construct_linear(V, gfu_0, dt, time_step) + comp = self.model_components + zeta = V[comp['k']] + psi = V[comp['epsilon']] + wind = self._get_wind(gfu_0, time_step) + k, epsilon = self._regularized_turbulence(time_step) + nu_t = self._get_turbulent_viscosity(time_step) + production = self._production(time_step) + form = forms[0] + + source_k = self.f.get('k', [0.0] * len(self.t_param))[time_step] + source_epsilon = self.f.get( + 'epsilon', [0.0] * len(self.t_param))[time_step] + form += dt * ( + production + source_k + ) * zeta * ngs.dx + form += dt * ( + self.C_1 * self._epsilon_k_ratio(k, epsilon) * production + + source_epsilon + ) * psi * ngs.dx + + if self.DG: + n = self._normal + d_k = self.kv[time_step] + nu_t / self.sigma_k + d_epsilon = ( + self.kv[time_step] + nu_t / self.sigma_epsilon) + for component, test, diffusivity in ( + ('k', zeta, d_k), ('epsilon', psi, d_epsilon)): + for marker, values in self.BC.get( + 'dirichlet', {}).get(component, {}).items(): + value = values[time_step] + wind_n = wind * n + form += -dt * test * ( + 0.5 * value * wind_n + - 0.5 * value * ngs.Norm(wind_n) + ) * self._ds(marker) + form += dt * diffusivity * self._penalty * value * test * self._ds(marker) + form += -dt * diffusivity * value * ( + ngs.grad(test) * n) * self._ds(marker) + + if self._wallf is not None: + # Impose the wall-law epsilon through the same weak DG/Nitsche + # terms as configured Dirichlet data, using Picard-lagged k. + value = self._wallf.epsilon_wall_cell(k) + wind_n = wind * n + form += -dt * psi * ( + 0.5 * value * wind_n + - 0.5 * value * ngs.Norm(wind_n) + ) * self._ds(self.wall_boundary) + form += (dt * d_epsilon * self._penalty * value * psi + * self._ds(self.wall_boundary)) + form += -dt * d_epsilon * value * ( + ngs.grad(psi) * n) * self._ds(self.wall_boundary) + + for component, test in (('k', zeta), ('epsilon', psi)): + for marker, values in self.BC.get( + 'neumann', {}).get(component, {}).items(): + if (component == 'epsilon' and self._wallf is not None + and marker == self.wall_boundary): + continue + form += -dt * test * values[time_step] * self._ds(marker) + + forms[0] = form + return forms + + def solve_single_step(self, a_lst: List[BilinearForm], + L_lst: List[LinearForm], + precond_lst: List[Preconditioner], + gfu: GridFunction, time_step: int = 0) -> None: + comp = self.model_components + if (gfu.components[comp['k']].vec.Norm() == 0.0 + or gfu.components[comp['epsilon']].vec.Norm() == 0.0): + gfu.vec.data = self.IC.vec + + previous = ngs.GridFunction(self.fes) + + for iteration in range(self.nonlinear_max_iters): + previous.vec.data = gfu.vec + self.UIter.vec.data = gfu.vec + self.W[0].vec.data = gfu.components[comp['u']].vec + + if self._wallf is not None: + self._update_wall_function() + + self.apply_dirichlet_bcs_to(gfu, time_step) + a_lst[0].Assemble() + L_lst[0].Assemble() + if precond_lst[0] is not None: + precond_lst[0].Update() + self.linear_solve(a_lst[0], L_lst[0], precond_lst[0], gfu) + + for index, factor in enumerate(self.relaxation_factors): + if factor < 1.0: + gfu.components[index].vec.data = ( + factor * gfu.components[index].vec + + (1.0 - factor) * previous.components[index].vec) + + self._bound_turbulence(gfu) + + difference = gfu.vec.CreateVector() + difference.data = gfu.vec - previous.vec + tolerance = ( + self.abs_nonlinear_tolerance + + self.rel_nonlinear_tolerance * gfu.vec.Norm()) + if difference.Norm() < tolerance: + logging.info( + 'KEpsilonINS converged in %d nonlinear iteration(s).', + iteration + 1) + break + else: + logging.warning( + 'KEpsilonINS did not converge within %d nonlinear iterations.', + self.nonlinear_max_iters) + + self.UIter.vec.data = gfu.vec + self.W[0].vec.data = gfu.components[comp['u']].vec + + def linearized_solve(self, a_assembled: BilinearForm, L_assembled: LinearForm, + precond: Preconditioner, gfu: GridFunction): + """Perform one stationary Picard update.""" + super().linearized_solve(a_assembled, L_assembled, precond, gfu) + + for index, factor in enumerate(self.relaxation_factors): + if factor < 1.0: + gfu.components[index].vec.data = ( + factor * gfu.components[index].vec + + (1.0 - factor) * self.UIter.components[index].vec) + + self._bound_turbulence(gfu) + + comp = self.model_components + error_squared = 0.0 + norm_squared = 0.0 + for component in ('u', 'k', 'epsilon'): + index = comp[component] + difference = gfu.components[index].vec.CreateVector() + difference.data = ( + gfu.components[index].vec - self.UIter.components[index].vec) + error_squared += difference.Norm() ** 2 + norm_squared += gfu.components[index].vec.Norm() ** 2 + + return error_squared ** 0.5, norm_squared ** 0.5 + + def update_linearization(self, gfu: GridFunction) -> None: + self._bound_turbulence(gfu) + super().update_linearization(gfu) + self.UIter.vec.data = gfu.vec + if self._wallf is not None: + self._update_wall_function() diff --git a/opencmp/models/stokes.py b/opencmp/models/stokes.py index d54ce91..8761be0 100644 --- a/opencmp/models/stokes.py +++ b/opencmp/models/stokes.py @@ -31,7 +31,9 @@ class Stokes(INS): """ def _define_bc_types(self) -> List[str]: - return ['dirichlet', 'stress', 'pinned'] + # Accept NEUMANN sections in shared multiphysics BC files. Entries for + # variables not present in Stokes (for example k and epsilon) are ignored. + return ['dirichlet', 'stress', 'pinned', 'neumann'] def _post_init(self) -> None: # TODO: see if this override is still needed when transient tests are added diff --git a/opencmp/solvers/nonlinear_mixing.py b/opencmp/solvers/nonlinear_mixing.py index d37d61f..66f634b 100644 --- a/opencmp/solvers/nonlinear_mixing.py +++ b/opencmp/solvers/nonlinear_mixing.py @@ -6,9 +6,11 @@ _solve() is responsible for everything else (BC application, re-assembly, convergence checks, gfu/x_curr updates). -All three schemes follow the same first-step convention as the original code: on iteration 2 -(the first real mixing step) every scheme performs linear mixing and initialises alpha from -the dominant-eigenvalue estimate. Scheme-specific logic takes over from iteration 3 onward. +LinearMixing, DiagBroyden, and Anderson follow the same first-step convention as the +original code: on iteration 2 (the first real mixing step) each initialises alpha from the +dominant-eigenvalue estimate. Scheme-specific logic takes over from iteration 3 onward. +NoMixing has no such state; it passes every iteration's update through unchanged, for models +(e.g. KEpsilonINS) that already apply their own per-component relaxation. """ import numpy as np @@ -124,10 +126,21 @@ def step(self, f_vec: BaseVector, x_prev_vec: BaseVector, num_iterations: int) - return dx_np +class NoMixing: + """dx = f (accept the model's own update, no extra damping).""" + + def __init__(self, **_) -> None: + pass + + def step(self, f_vec: BaseVector, x_prev_vec: BaseVector, num_iterations: int) -> np.ndarray: + return f_vec.FV().NumPy().copy() + + _SCHEMES = { 'LinearMixing': LinearMixing, 'DiagBroyden': DiagBroyden, 'Anderson': Anderson, + 'NoMixing': NoMixing, 'default': Anderson, } diff --git a/pytests/full_system/k_epsilon/bc_dir/bc_config b/pytests/full_system/k_epsilon/bc_dir/bc_config new file mode 100644 index 0000000..1f27a9c --- /dev/null +++ b/pytests/full_system/k_epsilon/bc_dir/bc_config @@ -0,0 +1,13 @@ +[DIRICHLET] +u = top -> [0.0, 0.0] + bottom -> [0.0, 0.0] + right -> [0.0, 0.0] + left -> [0.0, 0.0] +k = top -> 0.001 + bottom -> 0.001 + right -> 0.001 + left -> 0.001 +epsilon = top -> 0.001 + bottom -> 0.001 + right -> 0.001 + left -> 0.001 diff --git a/pytests/full_system/k_epsilon/config b/pytests/full_system/k_epsilon/config new file mode 100644 index 0000000..40fef94 --- /dev/null +++ b/pytests/full_system/k_epsilon/config @@ -0,0 +1,47 @@ +[MESH] +filename = pytests/mesh_files/coarse_large_square_4bcs.vol + +[DIM] +diffuse_interface_method = False + +[FINITE ELEMENT SPACE] +elements = u -> VectorH1 + p -> H1 + k -> H1 + epsilon -> H1 +interpolant_order = 2 + +[DG] +DG = False +interior_penalty_coefficient = 10.0 + +[SOLVER] +linear_solver = default +preconditioner = default +linearization_method = Oseen +nonlinear_solver = default +nonlinear_tolerance = relative -> 1e-6 + absolute -> 1e-10 +nonlinear_max_iterations = 5 +relaxation_factors = 1.0, 1.0, 1.0, 1.0 + +[TRANSIENT] +transient = True +scheme = implicit euler +time_range = 0.0, 0.001 +dt = 0.001 +dt_range = 1e-12, 0.001 +dt_tolerance = relative -> 1e-6 + absolute -> 1e-6 + +[VISUALIZATION] +save_to_file = False + +[ERROR ANALYSIS] +check_error = False + +[OTHER] +num_threads = 1 +model = KEpsilonINS +run_dir = pytests/full_system/k_epsilon +wall_function = False diff --git a/pytests/full_system/k_epsilon/ic_dir/ic_config b/pytests/full_system/k_epsilon/ic_dir/ic_config new file mode 100644 index 0000000..0e45816 --- /dev/null +++ b/pytests/full_system/k_epsilon/ic_dir/ic_config @@ -0,0 +1,5 @@ +[KEpsilonINS] +u = all -> [0.0, 0.0] +p = all -> 0.0 +k = all -> 0.001 +epsilon = all -> 0.001 diff --git a/pytests/full_system/k_epsilon/model_dir/model_config b/pytests/full_system/k_epsilon/model_dir/model_config new file mode 100644 index 0000000..da06bae --- /dev/null +++ b/pytests/full_system/k_epsilon/model_dir/model_config @@ -0,0 +1,19 @@ +[PARAMETERS] +kinematic_viscosity = all -> 1.0 +c_mu = all -> 0.09 +c_1 = all -> 1.44 +c_2 = all -> 1.92 +sigma_k = all -> 1.0 +sigma_epsilon = all -> 1.3 +k_floor = all -> 1e-10 +epsilon_floor = all -> 1e-10 +max_viscosity_ratio = all -> 1e5 +production_limit_coefficient = all -> 10.0 +kappa = all -> 0.41 + +[FUNCTIONS] +; A body force, so the no-slip cavity has a non-trivial flow and the wall +; function sees real shear. The smoke tests do not check against a reference. +source = u -> [1.0, 0.0] + k -> 0.0 + epsilon -> 0.0 diff --git a/pytests/full_system/k_epsilon/ref_sol_dir/ref_sol_config b/pytests/full_system/k_epsilon/ref_sol_dir/ref_sol_config new file mode 100644 index 0000000..b0ece3e --- /dev/null +++ b/pytests/full_system/k_epsilon/ref_sol_dir/ref_sol_config @@ -0,0 +1,3 @@ +[REFERENCE SOLUTIONS] + +[METRICS] diff --git a/pytests/full_system/k_epsilon/test_k_epsilon_full_system.py b/pytests/full_system/k_epsilon/test_k_epsilon_full_system.py new file mode 100644 index 0000000..10dec72 --- /dev/null +++ b/pytests/full_system/k_epsilon/test_k_epsilon_full_system.py @@ -0,0 +1,35 @@ +"""Full-system smoke coverage for the single-phase k-epsilon model.""" + +from opencmp.config_functions import ConfigParser +from opencmp.helpers.testing import run_example + + +def test_transient_cg_smoke() -> None: + config = ConfigParser('pytests/full_system/k_epsilon/config') + run_example(config) + + +def test_transient_dg_smoke() -> None: + config = ConfigParser('pytests/full_system/k_epsilon/config') + config['DG']['DG'] = 'True' + config['FINITE ELEMENT SPACE']['elements'] = ( + 'u -> HDiv\n' + 'p -> L2\n' + 'k -> L2\n' + 'epsilon -> L2' + ) + run_example(config) + + +def test_wall_function_smoke() -> None: + config = ConfigParser('pytests/full_system/k_epsilon/config') + config['DG']['DG'] = 'True' + config['FINITE ELEMENT SPACE']['elements'] = ( + 'u -> HDiv\n' + 'p -> L2\n' + 'k -> L2\n' + 'epsilon -> L2' + ) + config['OTHER']['wall_function'] = 'True' + config['OTHER']['wall_boundary'] = 'bottom' + run_example(config) diff --git a/pytests/helpers/test_k_epsilon_wall_func.py b/pytests/helpers/test_k_epsilon_wall_func.py new file mode 100644 index 0000000..1739402 --- /dev/null +++ b/pytests/helpers/test_k_epsilon_wall_func.py @@ -0,0 +1,508 @@ +"""Tests for the geometry-independent k-epsilon wall function. + +The point of ``KEpsilonWallFunction`` is that it needs no cylinder radius and no +hand-labelled near-wall/core materials, so most of these tests check that the +near-wall layer, the wall distance and the wall law come out right on meshes +that carry none of that information. +""" + +import numpy as np +import ngsolve as ngs +import pytest +from netgen.csg import unit_cube +from ngsolve.meshes import MakeStructured2DMesh, MakeStructured3DMesh + +from opencmp.helpers.wall_func import KEpsilonWallFunction + +SQUARE = 'pytests/mesh_files/unit_square_coarse.vol' +CHANNEL = 'pytests/mesh_files/channel_3bcs.vol' + +NU = 1.5 +E_LOG = 9.8 + + +def build(mesh, wall='bottom', **kwargs): + return KEpsilonWallFunction(mesh, nu=NU, C_mu=0.09, kappa=0.41, E_log=E_LOG, + wall_boundary=wall, **kwargs) + + +def velocity_field(mesh, value): + field = ngs.GridFunction(ngs.VectorH1(mesh, order=1)) + field.Set(ngs.CoefficientFunction(value)) + return field + + +def marked_element_numbers(wf): + """Element numbers of the marked cells, via the space's element->DOF map. + + Deliberately does not assume DOF number == element number. + """ + values = wf.mask.vec.FV().NumPy() + return {el.nr for el in wf._fes0.Elements(ngs.VOL) + if any(values[dof] > 0.5 for dof in el.dofs)} + + +def elements_owning_a_facet_on(mesh, wall): + """Ground truth: an element owns a wall facet when at least ``dim`` of its + vertices lie on that boundary (``dim - 1`` vertices is a vertex/edge touch).""" + wall_vertices = set() + for sel in mesh.Elements(ngs.BND): + if sel.mat == wall: + wall_vertices.update(v.nr for v in sel.vertices) + return {el.nr for el in mesh.Elements(ngs.VOL) + if sum(1 for v in el.vertices if v.nr in wall_vertices) >= mesh.dim} + + +def add_one_face_connected_layer(mesh, element_numbers): + """Ground truth for one topological dilation through volume-cell facets.""" + elements = list(mesh.Elements(ngs.VOL)) + facet_to_elements = {} + element_by_number = {element.nr: element for element in elements} + for element in elements: + for facet in element.facets: + facet_to_elements.setdefault(facet.nr, set()).add(element.nr) + + expanded = set(element_numbers) + for element_number in element_numbers: + for facet in element_by_number[element_number].facets: + expanded.update(facet_to_elements[facet.nr]) + return expanded + + +# ---------------------------------------------------------------------- +# Topology: the near-wall mask +# ---------------------------------------------------------------------- + +@pytest.mark.parametrize('meshfile, wall', [(SQUARE, 'bottom'), (CHANNEL, 'wall')]) +def test_mask_contains_only_wall_facet_owners(meshfile, wall): + mesh = ngs.Mesh(meshfile) + wf = build(mesh, wall) + wall_cells = elements_owning_a_facet_on(mesh, wall) + assert marked_element_numbers(wf) == wall_cells + + +def test_near_wall_mask_matches_wall_facet_owners(): + mesh = ngs.Mesh(CHANNEL) + wf = build(mesh, 'wall') + values = wf.near_wall_mask().vec.FV().NumPy() + actual = {element.nr for element in wf._fes0.Elements(ngs.VOL) + if any(values[dof] > 0.5 for dof in element.dofs)} + wall_cells = elements_owning_a_facet_on(mesh, 'wall') + assert actual == wall_cells + + +def test_wall_facet_mask_retains_the_true_boundary_owners(): + mesh = ngs.Mesh(CHANNEL) + wf = build(mesh, 'wall') + values = wf.wall_facet_mask().vec.FV().NumPy() + actual = {element.nr for element in wf._fes0.Elements(ngs.VOL) + if any(values[dof] > 0.5 for dof in element.dofs)} + assert actual == elements_owning_a_facet_on(mesh, 'wall') + + +def test_mask_is_not_empty_and_is_a_strict_subset(): + mesh = ngs.Mesh(CHANNEL) + wf = build(mesh, 'wall') + marked = marked_element_numbers(wf) + assert 0 < len(marked) < mesh.ne + + +def test_wall_measure_sums_to_the_exact_boundary_measure(): + """Every wall facet is attributed to exactly one cell -- no double counting, + none dropped.""" + mesh = ngs.Mesh(CHANNEL) + wf = build(mesh, 'wall') + exact = ngs.Integrate(ngs.CoefficientFunction(1.0), mesh, + definedon=mesh.Boundaries('wall')) + assert wf._wall_measure.sum() == pytest.approx(exact, rel=1e-12) + + +def test_cell_touching_the_wall_only_at_a_vertex_is_not_marked(): + """A vertex-only touch is not a wall-function cell.""" + mesh = ngs.Mesh(SQUARE) + wf = build(mesh, 'bottom') + marked = marked_element_numbers(wf) + + bottom_vertices = set() + for sel in mesh.Elements(ngs.BND): + if sel.mat == 'bottom': + bottom_vertices.update(v.nr for v in sel.vertices) + + vertex_only = {el.nr for el in mesh.Elements(ngs.VOL) + if sum(1 for v in el.vertices if v.nr in bottom_vertices) == 1} + assert vertex_only, 'mesh exercises no vertex-only touch; test is vacuous' + assert not vertex_only & marked + + +def test_missing_wall_marker_raises_a_clear_error(): + mesh = ngs.Mesh(SQUARE) + with pytest.raises(ValueError, match='no boundary elements'): + build(mesh, 'not_a_boundary') + + +# ---------------------------------------------------------------------- +# Wall distance +# ---------------------------------------------------------------------- + +def test_wall_distance_is_zero_on_the_wall_and_grows_inward(): + mesh = ngs.Mesh(SQUARE) + wf = build(mesh, 'bottom') + dist = wf.wall_distance_field() + + assert abs(dist(mesh(0.5, 0.0))) < 1e-8 + samples = [dist(mesh(0.5, y)) for y in (0.1, 0.3, 0.6, 0.9)] + assert all(b > a for a, b in zip(samples, samples[1:])) + assert samples[0] == pytest.approx(0.1, abs=0.05) + + +def _epsilon_wall_values(wf, k): + epsilon = ngs.GridFunction(wf._fes0) + epsilon.Set(wf.epsilon_wall_cell(ngs.CoefficientFunction(k))) + return epsilon.vec.FV().NumPy() + + +def _k_for_yplus(wf, yplus): + """k giving the requested y+ in the *first* cell, from y+ = y*Cmu^0.25*sqrt(k)/nu.""" + y = float(wf.wall_distance_cell().vec.FV().NumPy().min()) + return (yplus * NU / (0.09 ** 0.25 * y)) ** 2 + + +def test_epsilon_wall_uses_high_re_equilibrium_relation_in_the_log_layer(): + mesh = ngs.Mesh(SQUARE) + wf = build(mesh, 'bottom') + k = _k_for_yplus(wf, 5.0 * wf.YPLUS_VISCOUS) + distance = wf.wall_distance_cell().vec.FV().NumPy() + expected = 0.09 ** 0.75 * k ** 1.5 / (0.41 * distance) + + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + assert _epsilon_wall_values(wf, k)[owners] == pytest.approx( + expected[owners], rel=1e-9) + + +def test_epsilon_wall_switches_to_viscous_dissipation_below_yplus_lam(): + """Below y+_lam use the viscous-sublayer dissipation relation.""" + mesh = ngs.Mesh(SQUARE) + wf = build(mesh, 'bottom') + k = _k_for_yplus(wf, 0.1 * wf.YPLUS_VISCOUS) + distance = wf.wall_distance_cell().vec.FV().NumPy() + expected = 2.0 * NU * k / distance ** 2 + + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + assert _epsilon_wall_values(wf, k)[owners] == pytest.approx( + expected[owners], rel=1e-9) + + +def wall_law_nu_t(yplus): + """Log-law wall viscosity: nu*(y+/u+ - 1) with u+ = ln(E*y+)/kappa.""" + return NU * (yplus / (np.log(E_LOG * yplus) / 0.41) - 1.0) + + +def test_marked_cells_are_the_closest_cells_to_the_wall(): + mesh = ngs.Mesh(CHANNEL) + wf = build(mesh, 'wall') + dist = wf.wall_distance_cell().vec.FV().NumPy() + marked = wf.mask.vec.FV().NumPy() > 0.5 + assert dist[marked].max() < dist[~marked].max() + + +# ---------------------------------------------------------------------- +# Viscosity selection +# ---------------------------------------------------------------------- + +def force_yplus(wf, target): + """Return cellwise k values that give target y+ on every marked cell.""" + dist = wf.wall_distance_cell().vec.FV().NumPy() + marked = wf.mask.vec.FV().NumPy() > 0.5 + k = ngs.GridFunction(wf._fes0) + values = k.vec.FV().NumPy() + values[:] = 1.0 + values[marked] = (target * NU / (wf.C_mu ** 0.25 * dist[marked])) ** 2 + return marked, k + + +def sample_nu_t(wf, mesh, k, eps): + """nu_t as a cellwise array, via an L2(0) projection.""" + K = ngs.CoefficientFunction(k) + E = ngs.CoefficientFunction(eps) + out = ngs.GridFunction(wf._fes0) + out.Set(wf.eval_nu_t(K, E)) + return out.vec.FV().NumPy() + + +def test_compiled_viscosity_matches_symbolic_viscosity(): + mesh = ngs.Mesh(CHANNEL) + wf = build(mesh, 'wall') + _, k = force_yplus(wf, 50.0) + epsilon = ngs.CoefficientFunction(1.0) + symbolic = wf.eval_nu_t(k, epsilon) + compiled = symbolic.Compile() + symbolic_values = ngs.GridFunction(wf._fes0) + compiled_values = ngs.GridFunction(wf._fes0) + symbolic_values.Set(symbolic) + compiled_values.Set(compiled) + + assert compiled_values.vec.FV().NumPy() == pytest.approx( + symbolic_values.vec.FV().NumPy(), rel=1e-12, abs=1e-12) + + +@pytest.fixture +def channel_wf(): + mesh = ngs.Mesh(CHANNEL) + return mesh, build(mesh, 'wall') + + +def test_low_yplus_selects_zero_wall_viscosity(channel_wf): + mesh, wf = channel_wf + marked, k = force_yplus(wf, 5.0) # below 11.25 + nu_t = sample_nu_t(wf, mesh, k=k, eps=1.0) + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + assert nu_t[owners] == pytest.approx(0.0, abs=1e-12) + + +def test_log_range_yplus_selects_the_log_law(channel_wf): + mesh, wf = channel_wf + yplus = 50.0 + marked, k = force_yplus(wf, yplus) + nu_t = sample_nu_t(wf, mesh, k=k, eps=1.0) + expected = wall_law_nu_t(yplus) + assert expected > 0 + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + assert nu_t[owners] == pytest.approx(expected, rel=1e-6) + + +def test_wall_viscosity_is_constant_within_each_wall_cell(channel_wf): + """One P0 y+ selects one wall-law branch and viscosity per wall cell.""" + mesh, wf = channel_wf + _, k = force_yplus(wf, 60.0) + nu_t = wf.eval_nu_wall(ngs.CoefficientFunction(k), ngs.CoefficientFunction(1.0)) + + cellwise_mean = ngs.GridFunction(wf._fes0) + cellwise_mean.Set(nu_t) + spread = ngs.Integrate((nu_t - cellwise_mean) ** 2, mesh) + magnitude = ngs.Integrate(cellwise_mean ** 2, mesh) + assert np.sqrt(spread / magnitude) < 1e-12 + + +def test_wall_viscosity_is_continuous_at_the_sublayer_threshold(channel_wf): + """The log law is <= 0 at YPLUS_VISCOUS and clamped to zero, so the sublayer + branch joins it continuously.""" + mesh, wf = channel_wf + threshold = wf.YPLUS_VISCOUS + samples = [] + for yplus in (threshold - 1e-5, threshold + 1e-5): + marked, k = force_yplus(wf, yplus) + values = sample_nu_t(wf, mesh, k=k, eps=wf.epsilon_wall_cell(k)) + samples.append(values[marked]) + + scale = np.maximum(np.maximum(np.abs(samples[0]), np.abs(samples[1])), 1.0) + assert np.max(np.abs(samples[1] - samples[0]) / scale) < 1e-4 + + +def test_yplus_above_200_continues_to_use_the_log_law(channel_wf): + mesh, wf = channel_wf + marked, k = force_yplus(wf, 5.0 * wf.YPLUS_RECOMMENDED_MAX) + eps = 2.5 + nu_t = sample_nu_t(wf, mesh, k=k, eps=eps) + + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + assert nu_t[owners] == pytest.approx( + wall_law_nu_t(5.0 * wf.YPLUS_RECOMMENDED_MAX), rel=1e-6) + + +def test_warns_once_when_wall_yplus_is_above_recommended_band(channel_wf, caplog): + _, wf = channel_wf + _, k = force_yplus(wf, 400.0) + + wf.update(k) + wf.update(k) + + messages = [record.message for record in caplog.records + if 'Wall-function validity' in record.message] + assert len(messages) == 1 + assert 'y+ > 300' in messages[0] + assert '100.0%' in messages[0] + + +def test_wall_law_holds_across_the_whole_log_layer(channel_wf): + mesh, wf = channel_wf + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + for yplus in (20.0, 60.0, 150.0, 199.0): + _, k = force_yplus(wf, yplus) + nu_t = sample_nu_t(wf, mesh, k=k, eps=2.5) + assert nu_t[owners] == pytest.approx(wall_law_nu_t(yplus), rel=1e-6), yplus + + +def test_unmarked_cells_use_bulk_even_when_their_yplus_is_low(channel_wf): + """Every non-wall-owner cell uses bulk k-epsilon regardless of y+.""" + mesh, wf = channel_wf + marked, k = force_yplus(wf, 5.0) # wall term would be 0 + k.vec.FV().NumPy()[~marked] = 1e-12 # also low y+ outside mask + eps = 2.5 + nu_t = sample_nu_t(wf, mesh, k=k, eps=eps) + + kvals = k.vec.FV().NumPy() + bulk = 0.09 * kvals ** 2 / eps + assert nu_t[~marked] == pytest.approx(bulk[~marked], rel=1e-6) + owners = wf.wall_facet_mask().vec.FV().NumPy() > 0.5 + assert nu_t[owners] == pytest.approx(0.0, abs=1e-12) + + yplus = ngs.GridFunction(wf._fes0) + yplus.Set(wf.y_plus_cell(k)) + assert yplus.vec.FV().NumPy()[~marked].min() < 11.25, \ + 'no low-y+ unmarked cell; test is vacuous' + + +def test_bulk_viscosity_matches_the_plain_k_epsilon_formula(channel_wf): + mesh, wf = channel_wf + k, eps = 0.8, 1.3 + nu_t = sample_nu_t(wf, mesh, k=k, eps=eps) + marked = wf.mask.vec.FV().NumPy() > 0.5 + assert nu_t[~marked] == pytest.approx(0.09 * k ** 2 / eps, rel=1e-6) + + +def test_wall_owner_neighbours_are_not_in_the_wall_function_mask(channel_wf): + mesh, wf = channel_wf + owners = elements_owning_a_facet_on(mesh, 'wall') + expanded = add_one_face_connected_layer(mesh, owners) + neighbours = expanded - owners + assert neighbours + + k = ngs.GridFunction(wf._fes0) + k.vec.FV().NumPy()[:] = 1.0 + wf.update(k) + mask = wf.near_wall_mask().vec.FV().NumPy() + u_tau = wf.u_tau_cell.vec.FV().NumPy() + for element in wf._fes0.Elements(ngs.VOL): + if element.nr in neighbours: + assert mask[list(element.dofs)] == pytest.approx(0.0, abs=1e-14) + assert u_tau[list(element.dofs)] == pytest.approx(0.0, abs=1e-14) + + +def test_velocity_method_uses_molecular_tangential_not_normal_stress(): + mesh = ngs.Mesh(SQUARE) + wf = build(mesh, 'bottom', u_tau_method=1) + + wf.update(ngs.CoefficientFunction(1.0), + velocity_field(mesh, (0.0, ngs.y))) + assert wf.u_tau_cell.vec.FV().NumPy()[wf._marked] == pytest.approx( + 0.0, abs=1e-14) + + wf.update(ngs.CoefficientFunction(1.0), + velocity_field(mesh, (ngs.y, 0.0))) + owners = wf._marked + assert wf.wall_shear_cell.vec.FV().NumPy()[owners] == pytest.approx( + NU, rel=1e-12) + assert wf.u_tau_cell.vec.FV().NumPy()[owners] == pytest.approx( + np.sqrt(NU), rel=1e-12) + + +def test_velocity_method_is_invariant_when_wall_and_velocity_are_rotated(): + mesh = ngs.Mesh(SQUARE) + horizontal = build(mesh, 'bottom', u_tau_method=1) + vertical = build(mesh, 'left', u_tau_method=1) + k = ngs.CoefficientFunction(1.0) + + horizontal.update(k, velocity_field(mesh, (ngs.y, 0.0))) + vertical.update(k, velocity_field(mesh, (0.0, ngs.x))) + + assert horizontal.u_tau_cell.vec.FV().NumPy()[horizontal._marked] == \ + pytest.approx(vertical.u_tau_cell.vec.FV().NumPy()[vertical._marked], + rel=1e-12) + + +def test_velocity_method_requires_a_velocity_iterate(): + mesh = ngs.Mesh(SQUARE) + wf = build(mesh, 'bottom', u_tau_method=1) + with pytest.raises(ValueError, match='requires the velocity'): + wf.update(ngs.CoefficientFunction(1.0)) + + +def test_negative_viscosity_is_clamped_to_zero(channel_wf): + mesh, wf = channel_wf + nu_t = sample_nu_t(wf, mesh, k=1.0, eps=-1.0) # negative bulk + assert nu_t.min() >= 0.0 + + +# ---------------------------------------------------------------------- +# Geometry independence +# ---------------------------------------------------------------------- + +def test_runs_on_a_mesh_with_no_named_regions_or_cylinder_radius(): + """channel_3bcs has a single 'default' material and no near-wall/core + labelling; the legacy class could not have been built on it.""" + mesh = ngs.Mesh(CHANNEL) + assert set(mesh.GetMaterials()) == {'default'} + wf = build(mesh, 'wall') + wf.update(ngs.CoefficientFunction(1.0)) + assert wf.u_tau_cell.vec.FV().NumPy().max() > 0 + + +# ---------------------------------------------------------------------- +# Element type: simplices only +# ---------------------------------------------------------------------- + +@pytest.fixture +def cube_wf(): + """Unit cube of tets, wall on the z = 0 face.""" + mesh = ngs.Mesh(unit_cube.GenerateMesh(maxh=0.25)) + return mesh, build(mesh, 'bottom') + + +def test_tet_mesh_attributes_every_wall_face_to_exactly_one_cell(cube_wf): + mesh, wf = cube_wf + exact = ngs.Integrate(ngs.CoefficientFunction(1.0), mesh, + definedon=mesh.Boundaries('bottom')) + assert wf._wall_measure.sum() == pytest.approx(exact, rel=1e-12) + + +def test_tet_mask_contains_only_wall_face_owners(cube_wf): + mesh, wf = cube_wf + wall_cells = elements_owning_a_facet_on(mesh, 'bottom') + assert wall_cells, 'no cell owns a wall face; test is vacuous' + assert marked_element_numbers(wf) == wall_cells + + +def test_tet_marked_cells_are_the_closest_cells_to_the_wall(cube_wf): + _, wf = cube_wf + dist = wf.wall_distance_cell().vec.FV().NumPy() + marked = wf.mask.vec.FV().NumPy() > 0.5 + assert dist[marked].max() < dist[~marked].max() + + +def test_tet_wall_distance_is_zero_on_the_wall_and_grows_inward(cube_wf): + mesh, wf = cube_wf + dist = wf.wall_distance_field() + assert abs(dist(mesh(0.5, 0.5, 0.0))) < 1e-8 + samples = [dist(mesh(0.5, 0.5, z)) for z in (0.1, 0.2, 0.3)] + assert all(b > a for a, b in zip(samples, samples[1:])) + assert samples[0] == pytest.approx(0.1, abs=0.05) + + +def test_tet_friction_velocity_is_the_equilibrium_value_on_owners(cube_wf): + _, wf = cube_wf + k = 0.5 + wf.update(ngs.CoefficientFunction(k)) + u_tau = wf.u_tau_cell.vec.FV().NumPy() + assert u_tau[wf._marked] == pytest.approx(0.09 ** 0.25 * np.sqrt(k), rel=1e-9) + assert u_tau[wf.mask.vec.FV().NumPy() < 0.5].max() == 0.0 + + +def test_tet_cells_outside_the_mask_use_bulk_viscosity(cube_wf): + mesh, wf = cube_wf + k, eps = 0.5, 1.0 + wf.update(ngs.CoefficientFunction(k)) + nu_t = ngs.GridFunction(wf._fes0) + nu_t.Set(wf.eval_nu_t(ngs.CoefficientFunction(k), ngs.CoefficientFunction(eps))) + outside = wf.mask.vec.FV().NumPy() < 0.5 + assert nu_t.vec.FV().NumPy()[outside] == pytest.approx(0.09 * k ** 2 / eps, rel=1e-9) + + +@pytest.mark.parametrize('mesh_factory, kind', [ + (lambda: MakeStructured2DMesh(quads=True, nx=4, ny=4), 'QUAD'), + (lambda: MakeStructured3DMesh(hexes=True, nx=3, ny=3, nz=3), 'HEX'), +]) +def test_non_simplicial_meshes_are_refused(mesh_factory, kind): + """Wall functions are not implemented for quads or hexes.""" + mesh = mesh_factory() + with pytest.raises(NotImplementedError, match=kind): + build(mesh, 'bottom') diff --git a/pytests/models/test_k_epsilon.py b/pytests/models/test_k_epsilon.py new file mode 100644 index 0000000..98bff6b --- /dev/null +++ b/pytests/models/test_k_epsilon.py @@ -0,0 +1,355 @@ +"""Structural regression tests for the k-epsilon INS model.""" + +import inspect +import numpy as np +import re +from types import SimpleNamespace + +import ngsolve as ngs +import pytest +from netgen.geom2d import unit_square + +from opencmp.helpers.limiter import Limiter +from opencmp.helpers.wall_func import KEpsilonWallFunction +from opencmp.models import KEpsilonINS, models_dict +from opencmp.models.ins import INS + + +def _auto_inlet_model(explicit_ic=(), explicit_bc=()) -> KEpsilonINS: + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + model = object.__new__(KEpsilonINS) + model.mesh = mesh + model.auto_turbulence_inlet = 'left' + model.turbulence_hydraulic_diameter = 2.0 + model.turbulence_length_scale_ratio = 0.07 + model.kv = [0.1] + model.C_mu = 0.09 + model.t_param = [ngs.Parameter(0.0)] + model.model_components_ic = {'u': 0, 'p': 1, 'k': 2, 'epsilon': 3} + model.BC = {'dirichlet': { + 'u': {'left': [ngs.CoefficientFunction((1.0, 0.0))]}, + }} + for component in explicit_bc: + model.BC['dirichlet'][component] = {'left': [9.0]} + model.dirichlet_names = {'u': 'left'} + model.ic_functions = SimpleNamespace(ic_dict={ + model.name(): {component: {'all': [9.0]} for component in explicit_ic} + }) + spaces = [ngs.L2(mesh, order=0) for _ in range(4)] + model.IC = ngs.GridFunction(ngs.FESpace(spaces)) + for component in explicit_ic: + model.IC.components[model.model_components_ic[component]].Set(9.0) + return model + + +def _turbulence_stub(k_value: float, epsilon_value: float, + ratio: float = 2000.0) -> KEpsilonINS: + """Bare model carrying only what _regularized_turbulence reads.""" + model = object.__new__(KEpsilonINS) + model.model_components = {'u': 0, 'p': 1, 'k': 2, 'epsilon': 3} + model._bounded = False + model.k_floor = 1e-8 + model.epsilon_floor = 1e-4 + model.C_mu = 0.09 + model.kv = [2e-5] + model.max_viscosity_ratio = ratio + model.UIter = SimpleNamespace(components=[ + None, None, + ngs.CoefficientFunction(k_value), + ngs.CoefficientFunction(epsilon_value)]) + model._k_cell = ngs.CoefficientFunction(k_value) + model._epsilon_cell = ngs.CoefficientFunction(epsilon_value) + # The stub carries no velocity, so the strain-based realizability bound + # cannot be evaluated; test_realizability_* builds its own model for that. + model.realizability_limiter = False + return model + + +def test_k_epsilon_ins_is_registered_as_an_ins_model() -> None: + assert models_dict['KEpsilonINS'] is KEpsilonINS + assert issubclass(KEpsilonINS, INS) + + +def test_k_epsilon_ins_component_contract() -> None: + model = object.__new__(KEpsilonINS) + + assert model._define_model_components() == { + 'u': 0, + 'p': 1, + 'k': 2, + 'epsilon': 3, + } + assert model._define_time_derivative_components() == [{ + 'u': True, + 'p': False, + 'k': True, + 'epsilon': True, + }] + + +def test_k_epsilon_supports_scalar_neumann_boundaries() -> None: + model = object.__new__(KEpsilonINS) + + assert 'neumann' in model._define_bc_types() + + +def test_k_epsilon_reuses_cached_turbulent_viscosity() -> None: + model = object.__new__(KEpsilonINS) + cached_viscosity = object() + model._turbulent_viscosity = [cached_viscosity] + + assert model._get_turbulent_viscosity(0) is cached_viscosity + + +def test_epsilon_bound_caps_the_eddy_viscosity_ratio() -> None: + """The epsilon bound keeps nu_t/nu below the configured ratio.""" + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + k, epsilon = _turbulence_stub(0.0149, 1e-4)._regularized_turbulence(0) + ratio = ngs.Integrate(0.09 * k ** 2 / epsilon, mesh) / 2e-5 + assert ratio == pytest.approx(2000.0, rel=1e-9) + + k, epsilon = _turbulence_stub(0.0149, 3e-3)._regularized_turbulence(0) + assert ngs.Integrate(epsilon, mesh) == pytest.approx(3e-3, rel=1e-9) + + +def test_epsilon_k_ratio_is_smoothly_bounded() -> None: + model = object.__new__(KEpsilonINS) + model.max_epsilon_k_ratio = 10.0 + model._bounded = False + + regular = model._epsilon_k_ratio(1.0, 0.1) + floor_state = model._epsilon_k_ratio(1e-8, 1e-4) + + assert regular == pytest.approx(0.1 / 1.01) + assert floor_state < 10.0 + assert floor_state == pytest.approx(1e-4 / 1.001e-5) + + +def test_ins_viscosity_hook_preserves_laminar_behavior() -> None: + model = object.__new__(INS) + model.kv = [1.25, 2.5] + + assert model._get_effective_viscosity(0) == 1.25 + assert model._get_effective_viscosity(1) == 2.5 + + +def test_turbulence_constants_fall_back_to_standard_values() -> None: + """Omitting a constant from [PARAMETERS] must use the prescribed value.""" + model = object.__new__(KEpsilonINS) + + for name, default in KEpsilonINS.DEFAULT_PARAMETERS.items(): + assert model._parameter({}, name) == default + assert model._parameter({'c_mu': {'all': [0.1, 0.1]}}, 'c_mu') == 0.1 + + +def test_every_constant_the_model_reads_has_a_default() -> None: + """_set_model_parameters must not KeyError on a config that omits constants.""" + source = inspect.getsource(KEpsilonINS._set_model_parameters) + requested = set(re.findall(r"_parameter\(\s*parameters,\s*'([a-z_]+)'", source)) + + assert requested, 'parameter reads no longer match the expected pattern' + assert requested <= set(KEpsilonINS.DEFAULT_PARAMETERS) + + +def test_k_based_wall_friction_velocity_remains_the_default() -> None: + assert KEpsilonINS.DEFAULT_PARAMETERS['wall_u_tau_method'] == 0.0 + + +def test_auto_turbulence_defaults_use_inlet_bulk_velocity() -> None: + model = _auto_inlet_model() + model._apply_auto_turbulence_defaults() + + reynolds = 1.0 * 2.0 / 0.1 + intensity = 0.16 * reynolds ** (-1.0 / 8.0) + expected_k = 1.5 * intensity ** 2 + expected_epsilon = (0.09 ** 0.75 * expected_k ** 1.5 + / (0.07 * 2.0)) + + assert model.BC['dirichlet']['k']['left'] == pytest.approx([expected_k]) + assert model.BC['dirichlet']['epsilon']['left'] == pytest.approx( + [expected_epsilon]) + assert model.dirichlet_names['k'] == 'left' + assert model.dirichlet_names['epsilon'] == 'left' + assert np.asarray(model.IC.components[2].vec).mean() == pytest.approx( + expected_k) + assert np.asarray(model.IC.components[3].vec).mean() == pytest.approx( + expected_epsilon) + + +def test_auto_turbulence_defaults_preserve_explicit_values() -> None: + model = _auto_inlet_model( + explicit_ic=('k', 'epsilon'), explicit_bc=('k', 'epsilon')) + model._apply_auto_turbulence_defaults() + + assert model.BC['dirichlet']['k']['left'] == [9.0] + assert model.BC['dirichlet']['epsilon']['left'] == [9.0] + assert np.asarray(model.IC.components[2].vec).mean() == pytest.approx(9.0) + assert np.asarray(model.IC.components[3].vec).mean() == pytest.approx(9.0) + + +def test_auto_turbulence_defaults_reject_non_inward_flux() -> None: + model = _auto_inlet_model() + model.BC['dirichlet']['u']['left'] = [ + ngs.CoefficientFunction((-1.0, 0.0))] + + with pytest.raises(ValueError, match='nonzero inward net flux'): + model._apply_auto_turbulence_defaults() + + +def test_k_epsilon_wall_function_is_enabled_by_default() -> None: + class EmptyConfig: + def get_item(self, *_args, **_kwargs): + raise KeyError + + model = object.__new__(KEpsilonINS) + model.config = EmptyConfig() + model._pre_init() + + assert model.wall_function is True + assert model.wall_boundary == 'wall' + + +def test_pre_init_sets_every_documented_option() -> None: + """_pre_init drives itself from DEFAULT_OPTIONS, so the dict is the contract.""" + class EmptyConfig: + def get_item(self, *_args, **_kwargs): + raise KeyError + + model = object.__new__(KEpsilonINS) + model.config = EmptyConfig() + model._pre_init() + + for key, default in KEpsilonINS.DEFAULT_OPTIONS.items(): + assert getattr(model, key) == default + + +def _limiter_stub(bounded: bool) -> KEpsilonINS: + """Bare model carrying only what _regularized_turbulence reads.""" + model = _turbulence_stub(0.5, 1.0) + model._bounded = bounded + return model + + +def test_recovered_coefficient_floors_apply_when_solution_is_bounded() -> None: + """Recovered closure fields retain a final coefficient-level safeguard.""" + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + model = _limiter_stub(True) + model._k_cell = ngs.CoefficientFunction(-1.0) + + k, _ = model._regularized_turbulence(0) + assert ngs.Integrate(k, mesh) == pytest.approx(model.k_floor, rel=1e-9) + + +def test_coefficient_floors_still_apply_on_unbounded_spaces() -> None: + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + model = _limiter_stub(False) + model._k_cell = ngs.CoefficientFunction(-1.0) + + k, _ = model._regularized_turbulence(0) + assert ngs.Integrate(k, mesh) == pytest.approx(model.k_floor, rel=1e-9) + + +def test_bound_epsilon_applies_in_both_modes() -> None: + """boundEpsilon depends on the local k, so no limiter can impose it; it must + survive whether or not the slope limiter is on.""" + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + for bounded in (True, False): + model = _limiter_stub(bounded) + _, epsilon = model._regularized_turbulence(0) + ratio = ngs.Integrate(0.09 * 0.5 ** 2 / epsilon, mesh) / 2e-5 + assert ratio <= 2000.0 + 1e-6, bounded + + +def test_bezier_bound_holds_k_above_its_floor_inside_the_element() -> None: + """What the dropped coefficient guards relied on: the limited polynomial is above + the floor everywhere in the element, not just at the DOFs.""" + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + fes = ngs.L2(mesh, order=1) + field = ngs.GridFunction(fes) + field.Set(ngs.x - 0.5) # negative over half the domain + floor = 1e-8 + + Limiter(mesh).bezier_bound(field, fes, fes.globalorder, (floor, 1e20)) + + # Sample the limited polynomial densely rather than trusting the DOFs. + below = ngs.Integrate(ngs.IfPos(floor - field, 1.0, 0.0), mesh) + assert below == pytest.approx(0.0, abs=1e-12) + + +def test_floored_epsilon_does_not_inflate_bulk_viscosity_to_the_cap() -> None: + """A bounded epsilon sitting at its floor is a clamped undershoot, not physics. + + Trusting C_mu*k**2/epsilon there saturates nu_t at the viscosity cap and paints + cap-to-zero jumps against neighbouring cells (the backward-facing-step scar). + The bulk viscosity must fade out instead. + """ + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + model = _turbulence_stub(2.6e-3, 1e-8) # moderate k, epsilon at its floor + model._bounded = True + model._wallf = None + model.epsilon_floor = 1e-8 + + nu_t = ngs.Integrate(model._build_turbulent_viscosity(0), mesh) + cap = model.max_viscosity_ratio * model.kv[0] + assert nu_t < 0.05 * cap + + # A healthy epsilon must be essentially untouched by the trust factor. + model = _turbulence_stub(2.6e-3, 3e-4) + model._bounded = True + model._wallf = None + model.epsilon_floor = 1e-8 + nu_t = ngs.Integrate(model._build_turbulent_viscosity(0), mesh) + expected = 0.09 * 2.6e-3 ** 2 / 3e-4 + assert nu_t == pytest.approx(expected, rel=0.01) + + +def _realizability_model(k_value: float, epsilon_value: float, + mesh: ngs.Mesh, a: float = 1.0) -> KEpsilonINS: + """Stub with a real velocity field so the strain-rate bound can be evaluated.""" + model = _turbulence_stub(k_value, epsilon_value) + model.realizability_limiter = True + model.realizability_coefficient = a + model._wallf = None + model._bounded = False + fes = ngs.VectorH1(mesh, order=1) + u = ngs.GridFunction(fes) + u.Set(ngs.CoefficientFunction((ngs.y, 0.0))) # grad u = [[0,1],[0,0]] -> S = 1 + model.UIter = SimpleNamespace(components=[u, None, + ngs.CoefficientFunction(k_value), + ngs.CoefficientFunction(epsilon_value)]) + return model + + +def test_realizability_bound_caps_nu_t_by_the_strain_rate() -> None: + """Durbin: nu_t <= a*k/(sqrt(6)*S). With u=(y,0) the strain magnitude S is 1, + so a collapsing epsilon must not drive nu_t past a*k/sqrt(6).""" + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + k_value = 1e-2 + model = _realizability_model(k_value, 1e-12, mesh) # epsilon -> 0 + area = ngs.Integrate(ngs.CoefficientFunction(1.0), mesh) + nu_t = ngs.Integrate(model._build_turbulent_viscosity(0), mesh) / area + + expected = k_value / np.sqrt(6.0) + assert nu_t == pytest.approx(expected, rel=1e-6) + + # Without the bound this same state saturates the viscosity-ratio cap, + # which is the 0-to-cap jump that wrecks the high-order solve. + unlimited = _realizability_model(k_value, 1e-12, mesh) + unlimited.realizability_limiter = False + cap = unlimited.max_viscosity_ratio * unlimited.kv[0] + nu_t_unlimited = ngs.Integrate( + unlimited._build_turbulent_viscosity(0), mesh) / area + + assert nu_t_unlimited == pytest.approx(cap, rel=1e-6) + assert nu_t < nu_t_unlimited + + +def test_realizability_bound_is_inactive_when_bulk_nu_t_is_small() -> None: + """The bound is a max-of-two min; a healthy k-epsilon state must pass through.""" + mesh = ngs.Mesh(unit_square.GenerateMesh(maxh=0.5)) + k_value, epsilon_value = 1e-3, 1.0 # C_mu k^2/eps = 9e-11, tiny + model = _realizability_model(k_value, epsilon_value, mesh) + area = ngs.Integrate(ngs.CoefficientFunction(1.0), mesh) + nu_t = ngs.Integrate(model._build_turbulent_viscosity(0), mesh) / area + + assert nu_t == pytest.approx(0.09 * k_value ** 2 / epsilon_value, rel=1e-6)