diff --git a/docs/_static/tfm_anim.mp4 b/docs/_static/tfm_anim.mp4 new file mode 100644 index 0000000..ac8c718 Binary files /dev/null and b/docs/_static/tfm_anim.mp4 differ diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index 32200d0..d59d573 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -24,3 +24,4 @@ Before attempting the tutorials, make sure to review :ref:`getting_started` and, tutorial_9 tutorial_10 tutorial_11 + tutorial_12 diff --git a/docs/tutorials/tutorial_12.rst b/docs/tutorials/tutorial_12.rst new file mode 100644 index 0000000..53dd485 --- /dev/null +++ b/docs/tutorials/tutorial_12.rst @@ -0,0 +1,234 @@ +.. Contains the twelfth tutorial. +.. _tutorial_12: + +Tutorial 12 - Simulating a Two-Fluid-Model Bubble Column (2D) +================================================================ + +The files for this tutorial can be found in ``examples/TFM``. + +This tutorial demonstrates OpenCMP's transient ``TwoFluidModel`` using a +two-dimensional air--water bubble column. It is a reduced-dimensional version +of the three-dimensional configuration studied by Fazeli, Rhebergen, and +Abukhdeir [1]_. The equations, material properties, interphase-momentum +closures, and injection concept follow that work. + +Governing Equations +------------------- + +The continuous liquid phase and dispersed gas phase are treated as +interpenetrating continua. Their volume fractions satisfy +:math:`\alpha_c+\alpha_d=1`. Writing :math:`\alpha_d=1-\alpha_c`, the +constant-density phase mass balances can be expressed as + +.. math:: + \partial_t\alpha_c + +\nabla\cdot\bigl((\alpha_c-1)\boldsymbol{u}_d\bigr) + &= -S_d, \\ + \nabla\cdot\bigl(\alpha_c\boldsymbol{u}_c + +(1-\alpha_c)\boldsymbol{u}_d\bigr) + &= S_d. + +Here :math:`S_d` is nonzero only in the gas-injection region. The interphase +momentum exchange is equal and opposite between the phases. This example +includes drag, lift, virtual-mass, and laminar-dispersion contributions: + +.. math:: + \boldsymbol{\mathcal{M}}'_d + =\boldsymbol{F}_D+\boldsymbol{F}_L+\boldsymbol{F}_{VM} + +\boldsymbol{F}_{LD}, \qquad + \boldsymbol{\mathcal{M}}'_c=-\boldsymbol{\mathcal{M}}'_d. + +The effect of the laminar-dispersion force is investigated in [1]_. + +OpenCMP provides two canonical forms through ``canonical_form``. They use the +same mass balances, pressure, gravity, and interphase momentum closures, but +differ in their treatment of molecular momentum fluxes. The B-TFM formulation +is due to Brennen [2]_; its momentum equations in primitive form are + +.. math:: + \partial_t\boldsymbol{u}_c + +\boldsymbol{u}_c\cdot\nabla\boldsymbol{u}_c + &= -\frac{\nabla p}{\rho_c} + +\frac{1}{\alpha_c\rho_c}\nabla\cdot\boldsymbol{\tau}^{\mathrm{eff}}_c + +\boldsymbol{g} + -\frac{\boldsymbol{\mathcal{M}}'_d}{\alpha_c\rho_c}, \\ + \partial_t\boldsymbol{u}_d + +\boldsymbol{u}_d\cdot\nabla\boldsymbol{u}_d + &= -\frac{\nabla p}{\rho_d} + +\boldsymbol{g} + +\frac{\boldsymbol{\mathcal{M}}'_d}{\alpha_d\rho_d}. + +The classical C-TFM formulation described by Ishii and Hibiki [3]_ instead +retains phase-weighted molecular stress terms in both phase momentum equations: + +.. math:: + \partial_t\boldsymbol{u}_c + +\boldsymbol{u}_c\cdot\nabla\boldsymbol{u}_c + &= -\frac{\nabla p}{\rho_c} + +\frac{1}{\alpha_c\rho_c} + \nabla\cdot(\alpha_c\boldsymbol{\tau}^{\mathrm{eff}}_c) + +\boldsymbol{g} + -\frac{\boldsymbol{\mathcal{M}}'_d}{\alpha_c\rho_c}, \\ + \partial_t\boldsymbol{u}_d + +\boldsymbol{u}_d\cdot\nabla\boldsymbol{u}_d + &= -\frac{\nabla p}{\rho_d} + +\frac{1}{\alpha_d\rho_d} + \nabla\cdot(\alpha_d\boldsymbol{\tau}_d) + +\boldsymbol{g} + +\frac{\boldsymbol{\mathcal{M}}'_d}{\alpha_d\rho_d}. + +The principal distinction is therefore the dispersed-phase molecular stress +:math:`\nabla\cdot(\alpha_d\boldsymbol{\tau}_d)`, which is present in C-TFM +and absent in B-TFM. The continuous-phase stress is also averaged differently: +C-TFM uses :math:`\nabla\cdot(\alpha_c\boldsymbol{\tau}^{\mathrm{eff}}_c)`, +whereas B-TFM uses :math:`\nabla\cdot\boldsymbol{\tau}^{\mathrm{eff}}_c`. +For this tutorial, ``canonical_form = C-TFM`` selects the classical form. + +Geometry and Injection Configuration +------------------------------------ + +The Gmsh geometry in ``2D_sample.msh`` is a :math:`0.20\,\mathrm{m}` wide by +:math:`0.45\,\mathrm{m}` high column. A :math:`0.04\,\mathrm{m}` by +:math:`0.02\,\mathrm{m}` rectangular injection region is centred immediately +above the bottom boundary. The physical surface named ``injection`` activates +the dispersed-phase mass and momentum sources; the remainder is named +``surface``. Boundary curves are marked ``wall``, ``bottom``, and ``outlet``. + +The model configuration supplies the source parameters:: + + [INJECTION] + region = injection + mass_flow_rate = 0.3 + velocity = 0.2 + +Thus gas is introduced only in ``injection`` with an upward injection velocity +of :math:`0.2\,\mathrm{m/s}`. + +The Main Configuration File +--------------------------- + +The velocity spaces are H(div)-conforming, while pressure and continuous-phase +volume fraction use discontinuous L2 spaces:: + + [FINITE ELEMENT SPACE] + elements = u_c -> HDiv + u_d -> HDiv + p -> L2 + alpha_c -> L2 + interpolant_order = 3 + + [DG] + DG = True + interior_penalty_coefficient = 10.0 + +The nonlinear problem is advanced with implicit Euler and Picard iteration:: + + [SOLVER] + linear_solver = direct + linearization_method = Picard + nonlinear_tolerance = relative -> 1e-5 + absolute -> 1e-5 + nonlinear_max_iterations = 10 + + [TRANSIENT] + transient = True + scheme = implicit euler + time_range = 0.0, 3.0 + dt = 2e-3 + +``slope_limiter`` bounds the transported volume fraction, while +``diffusion_switch`` enables the configured artificial diffusion. The TFM +section selects the canonical form and interphase closures:: + + [OTHER] + model = TwoFluidModel + slope_limiter = True + diffusion_switch = True + + [TFM] + canonical_form = C-TFM + IME = drag -> Tomiyama + lift -> Tomiyama + virtual_mass -> ConstantCoefficient + laminar_dispersion -> ConstantCoefficient + lift_wall_deactivation = True + lift_wall_boundaries = wall|bottom + +The Model Configuration File +---------------------------- + +The physical properties represent :math:`3\,\mathrm{mm}` air bubbles in water:: + + [PARAMETERS] + rho_c = all -> 998.2 + rho_d = all -> 1.204 + nu_c = all -> 1.0038e-6 + nu_d = all -> 1.5158e-5 + sigma_c = all -> 0.072 + dp = all -> 3e-3 + c_vm = all -> 0.5 + cdis = all -> 4.545 + d_artificial = all -> 1e-5 + +Boundary and Initial Conditions +------------------------------- + +The liquid satisfies no slip on the side and bottom walls. The dispersed +phase has a slip condition there, constraining its normal velocity while +leaving tangential motion free. Both velocities have zero stress at the open +top, and ``ZERO_BACKFLOW`` supplies the outflow treatment for :math:`\alpha_c`:: + + [DIRICHLET] + u_c = wall -> [0.0, 0.0] + bottom -> [0.0, 0.0] + + [SLIP] + u_d = wall|bottom + + [ZERO_STRESS] + u_c = outlet + u_d = outlet + + [ZERO_BACKFLOW] + alpha_c = outlet|wall|bottom + +The column initially contains quiescent liquid and no gas:: + + [TwoFluidModel] + u_c = all -> [0.0, 0.0] + u_d = all -> [0.0, 0.0] + alpha_c = all -> 1.0 + +Running and Visualizing the Simulation +-------------------------------------- + +From ``examples/TFM``, run: + +.. code-block:: console + + python3 -m opencmp config + +The supplied configuration currently has ``resume_from_previous = True`` and +``restart_from = LATEST``. Keep these settings to continue an existing run; +disable ``resume_from_previous`` when starting from the initial conditions in a +clean output directory. Open the generated output in ParaView to inspect +:math:`\alpha_c`, :math:`\boldsymbol{u}_c`, and :math:`\boldsymbol{u}_d`. + +The following animation shows these three solution fields for the example: + +.. raw:: html + + + +.. [1] A. Fazeli, S. Rhebergen, and N. M. Abukhdeir, "Laminar dispersion + force effects on two-fluid modelling and simulation of bubble column + hydrodynamics," *International Journal of Multiphase Flow*, 105590, 2026. + +.. [2] C. E. Brennen, *Fundamentals of Multiphase Flow*, Cambridge University + Press, 2005. + +.. [3] M. Ishii and T. Hibiki, *Thermo-Fluid Dynamics of Two-Phase Flow*, + Springer Science & Business Media, 2010. doi:10.1007/978-1-4419-7985-8. diff --git a/examples/TFM/bc_2d.msh b/examples/TFM/bc_2d.msh new file mode 100644 index 0000000..0f076f6 --- /dev/null +++ b/examples/TFM/bc_2d.msh @@ -0,0 +1,862 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +6 +1 14 "wall" +1 15 "outlet" +1 16 "bottom" +1 17 "empty" +2 12 "injection" +2 13 "surface" +$EndPhysicalNames +$Nodes +269 +1 0 0 0 +2 0.2 0 0 +3 0.2 0.45 0 +4 0 0.45 0 +5 0.08 0.01 0 +6 0.12 0.01 0 +7 0.12 0.03 0 +8 0.08 0.03 0 +9 0.02857142857142862 0 0 +10 0.05714285714285725 0 0 +11 0.0857142857142858 0 0 +12 0.1142857142857144 0 0 +13 0.1428571428571429 0 0 +14 0.1714285714285715 0 0 +15 0.06000000000000001 0.007500000000000001 0 +16 0.04000000000000004 0.005000000000000005 0 +17 0.02000000000000005 0.002500000000000006 0 +18 0.09333333333333332 0.01 0 +19 0.1066666666666667 0.01 0 +20 0.14 0.007500000000000001 0 +21 0.16 0.005000000000000005 0 +22 0.18 0.002500000000000006 0 +23 0.12 0.02333333333333334 0 +24 0.12 0.01666666666666667 0 +25 0.09333333333333332 0.03 0 +26 0.1066666666666667 0.03 0 +27 0.08 0.02333333333333334 0 +28 0.08 0.01666666666666667 0 +29 0 0.02368421052631575 0 +30 0 0.04736842105263149 0 +31 0 0.07105263157894723 0 +32 0 0.09473684210526293 0 +33 0 0.1184210526315787 0 +34 0 0.1421052631578944 0 +35 0 0.1657894736842101 0 +36 0 0.1894736842105258 0 +37 0 0.2131578947368415 0 +38 0 0.2368421052631572 0 +39 0 0.260526315789473 0 +40 0 0.2842105263157889 0 +41 0 0.3078947368421047 0 +42 0 0.3315789473684206 0 +43 0 0.3552631578947364 0 +44 0 0.3789473684210524 0 +45 0 0.4026315789473683 0 +46 0 0.4263157894736841 0 +47 0.1714285714285714 0.45 0 +48 0.1428571428571428 0.45 0 +49 0.1142857142857142 0.45 0 +50 0.08571428571428563 0.45 0 +51 0.05714285714285708 0.45 0 +52 0.02857142857142855 0.45 0 +53 0.2 0.4263157894736843 0 +54 0.2 0.4026315789473685 0 +55 0.2 0.3789473684210528 0 +56 0.2 0.3552631578947371 0 +57 0.2 0.3315789473684213 0 +58 0.2 0.3078947368421057 0 +59 0.2 0.2842105263157899 0 +60 0.2 0.2605263157894742 0 +61 0.2 0.2368421052631585 0 +62 0.2 0.2131578947368428 0 +63 0.2 0.189473684210527 0 +64 0.2 0.1657894736842112 0 +65 0.2 0.1421052631578953 0 +66 0.2 0.1184210526315794 0 +67 0.2 0.09473684210526356 0 +68 0.2 0.07105263157894764 0 +69 0.2 0.04736842105263173 0 +70 0.2 0.02368421052631592 0 +71 0.1957894736842105 0.4278947368421052 0 +72 0.191578947368421 0.4057894736842105 0 +73 0.1873684210526316 0.3836842105263157 0 +74 0.1831578947368421 0.3615789473684209 0 +75 0.1789473684210526 0.3394736842105261 0 +76 0.1747368421052631 0.3173684210526314 0 +77 0.1705263157894736 0.2952631578947366 0 +78 0.1663157894736842 0.2731578947368418 0 +79 0.1621052631578947 0.251052631578947 0 +80 0.1578947368421052 0.2289473684210522 0 +81 0.1536842105263157 0.2068421052631575 0 +82 0.1494736842105263 0.1847368421052628 0 +83 0.1452631578947368 0.1626315789473681 0 +84 0.1410526315789473 0.1405263157894734 0 +85 0.1368421052631579 0.1184210526315787 0 +86 0.1326315789473684 0.09631578947368402 0 +87 0.1284210526315789 0.07421052631578934 0 +88 0.1242105263157895 0.05210526315789465 0 +89 0.08577350269189626 0.02 0 +90 0.1142264973081037 0.02 0 +91 0.09999999999999999 0.02 0 +92 0.1588171431305086 0.3329009797964643 0 +93 0.1454190081291423 0.2676389470200093 0 +94 0.02056538775687013 0.3204121794150639 0 +95 0.01993223615173632 0.2247069203979927 0 +96 0.02021717712721847 0.08298598788501069 0 +97 0.1262321237084461 0.1569491534372543 0 +98 0.08716935821711924 0.03863475769956584 0 +99 0.07714214004921592 0.03891303427093282 0 +100 0.08247158180091055 0.04786991268613369 0 +101 0.0994348944073236 0.4205684129101692 0 +102 0.02406520564652929 0.1318429500937496 0 +103 0.02100734787201835 0.2697593571899669 0 +104 0.01851852992889785 0.1780648309383702 0 +105 0.01975373751155829 0.3656157417764139 0 +106 0.1380828495313891 0.2223431538302859 0 +107 0.1201496047558009 0.1104539843316892 0 +108 0.07139464245883492 0.04753717343090719 0 +109 0.07739538243701875 0.05816568478517997 0 +110 0.1661194384426598 0.3762780017001449 0 +111 0.06465229353223484 0.05786888604733686 0 +112 0.07105341349439694 0.07001721252635634 0 +113 0.05861522941355696 0.04581978963585301 0 +114 0.07478423465478534 0.02596004423216106 0 +115 0.08580130556869289 0.06980106762279989 0 +116 0.07875897670187208 0.08385889455539031 0 +117 0.05933727478719727 0.08313739516659871 0 +118 0.06963180909921085 0.09946597589125211 0 +119 0.04373604708462311 0.4247000682430588 0 +120 0.05018461125170211 0.0576651049545695 0 +121 0.04383575155072134 0.04364036529518101 0 +122 0.05308423977787338 0.1007298456935292 0 +123 0.06026768400938082 0.1169388697726063 0 +124 0.09260069942883825 0.04675732479854042 0 +125 0.05241895952038644 0.033738109555809 0 +126 0.09577236979441861 0.08273154919100237 0 +127 0.1548718872951699 0.4250875442533762 0 +128 0.0796029183670445 0.1163210376820139 0 +129 0.07012956282174107 0.1349412197986927 0 +130 0.09186616799631503 0.1350618897644589 0 +131 0.08082346379129768 0.154234038702757 0 +132 0.05925710972582848 0.1536904705606629 0 +133 0.06969202572558973 0.1736850544991498 0 +134 0.09217784786145963 0.1740364107965329 0 +135 0.080893532004515 0.1935945814895841 0 +136 0.05801560449289274 0.1933902874772769 0 +137 0.06841461701762572 0.213333465533934 0 +138 0.09221546164463755 0.2134423478136739 0 +139 0.0806130551852524 0.233313636768646 0 +140 0.05774680448704785 0.2332046019544428 0 +141 0.06918854930778283 0.2532284043039585 0 +142 0.09223404834922018 0.2532713347736151 0 +143 0.08074331869067393 0.2731916274119177 0 +144 0.1035079659528782 0.2733372045806039 0 +145 0.09220337174764649 0.2931489864915742 0 +146 0.1149866296046554 0.2932934507286892 0 +147 0.1037865415050659 0.3130077170919205 0 +148 0.08051862029139184 0.3131938840073829 0 +149 0.09259596112878916 0.3337350478698562 0 +150 0.06869743384346712 0.3336422455724328 0 +151 0.08064636407954338 0.3545080076344458 0 +152 0.1050843827736617 0.3545223152537291 0 +153 0.09191242720320381 0.3751836218795056 0 +154 0.06805882971657123 0.3757805980726077 0 +155 0.1166948391081345 0.3768779509611103 0 +156 0.05796944430945299 0.3127298276054668 0 +157 0.05853108145104588 0.2735688391025518 0 +158 0.03585744416528873 0.05598092173235855 0 +159 0.1274644465238469 0.312392718689852 0 +160 0.1132072131726445 0.2544324249454356 0 +161 0.02999226382504142 0.02171427769454251 0 +162 0.1274279603995819 0.3527821904836079 0 +163 0.04142653826792067 0.07251832493186286 0 +164 0.08096915988832193 0.3980525607799223 0 +165 0.1373005010596673 0.2913819489525884 0 +166 0.04264719858715715 0.3328127997756822 0 +167 0.1032119882094014 0.06656374025721309 0 +168 0.05607839892235084 0.4001681134522417 0 +169 0.02008903317053173 0.1054653408212298 0 +170 0.1283147148669697 0.4250159207139266 0 +171 0.1422810506355163 0.4032093311405927 0 +172 0.1184113488758795 0.4015098368905009 0 +173 0.07051260719329369 0.4239148592308987 0 +174 0.03042852431995918 0.4018417860674787 0 +175 0.02580658107115223 0.04131184182483952 0 +176 0.01905045087616503 0.1558849221666856 0 +177 0.01954133817194477 0.2016779308198104 0 +178 0.0427833172331801 0.213255482718974 0 +179 0.1429900234574767 0.370761137260077 0 +180 0.1027751069466466 0.2339277023115111 0 +181 0.1177780713483441 0.2140515899438229 0 +182 0.1041777846194941 0.1940338781088829 0 +183 0.1263184518506628 0.1916492823489387 0 +184 0.02080652720360878 0.24666495796176 0 +185 0.02357190702433035 0.2944365679888086 0 +186 0.04880626575841359 0.1340662403175293 0 +187 0.0434480928423636 0.2553121506520857 0 +188 0.1674435892403259 0.4015919206630524 0 +189 0.0494466757227576 0.171629827729881 0 +190 0.1627689860766738 0.354438542905298 0 +191 0.1525080501262763 0.3097545231692506 0 +192 0.1429897772131071 0.2437450547252676 0 +193 0.1186553922880458 0.1329358844598852 0 +194 0.01993361944230817 0.3434833591832917 0 +195 0.06947313044938996 0.2931215233813735 0 +196 0.05664364406398878 0.3542249686913794 0 +197 0.04076684809711739 0.3769841735191442 0 +198 0.1155998534490028 0.3328717175066215 0 +199 0.03962794669564088 0.1155973602567471 0 +200 0.04768387176158576 0.2932899930945981 0 +201 0.1150451017153435 0.1750258018428739 0 +202 0.1041333495601513 0.154707196500627 0 +203 0.06676521619099966 0.03726118409468218 0 +204 0.1245564375672111 0.2731887057929904 0 +205 0.08920612742700816 0.09779840498492323 0 +206 0.1115290495922837 0.09148682220718604 0 +207 0.1372398283630715 0.3307903156507633 0 +208 0.1111060315077069 0.04027977608364816 0 +209 0.1762224938033398 0.4220727350885489 0 +210 0.03675572492997919 0.1873564672312474 0 +211 0.01961391361922765 0.06271363624274287 0 +212 0.05751442962374353 0.06836637076828966 0 +213 0.02135115161915691 0.4226540679944191 0 +214 0.03723294704337082 0.09030058861964865 0 +215 0.03738072457138537 0.2347954568901953 0 +216 0.124345419127193 0.2375719584957789 0 +217 0.07378245709528866 0.01887698788822231 0 +218 0.06408821928287393 0.02466342664566726 0 +219 0.03949008565136098 0.1513865714227176 0 +220 0.03802227607859741 0.3104213676140947 0 +221 0.03594900954042596 0.3546242085891823 0 +222 0.1007260709059902 0.1137296832592864 0 +223 0.09961981749232507 0.03912868911460388 0 +224 0.103243109254024 0.0504705216264261 0 +225 0.03884846019026879 0.2772733816056022 0 +226 0.1324665038759431 0.1741985317363396 0 +227 0.0451101341239674 0.02034447130451936 0 +228 0.01818982198572698 0.3852041297462915 0 +229 0.154886954952322 0.2867884496832667 0 +230 0.1458487882854625 0.3483346332192421 0 +231 0.09019370105560537 0.05700003083431768 0 +232 0.07434183967469142 0.03313355449413863 0 +233 0.1326932695607706 0.3885869061519017 0 +234 0.1008419574095459 0.3941221326582853 0 +235 0.03856866126946384 0.03211288700412127 0 +236 0.01515976897923875 0.02731575021966587 0 +237 0.1537061168927705 0.3894316045750065 0 +238 0.1312394965510263 0.2569318252183505 0 +239 0.03244490923001331 0.1675204846187883 0 +240 0.1373679689793669 0.2047043313246152 0 +241 0.1441660498856124 0.06033681059305995 0 +242 0.1298646262102409 0.02642767730534307 0 +243 0.1538805615573401 0.1037219982998548 0 +244 0.1825384064299202 0.03539276453578749 0 +245 0.152165063509461 0.02357050807568878 0 +246 0.1802188712308234 0.1327945365414866 0 +247 0.1810614466244464 0.0815351683196967 0 +248 0.1645008605880404 0.170452199398905 0 +249 0.1794888720156318 0.2013157894736849 0 +250 0.1810685457101859 0.05801879394362827 0 +251 0.1767328497359234 0.1067715647465455 0 +252 0.1486594521300853 0.08185845891986593 0 +253 0.1617458546646725 0.1479983448334832 0 +254 0.1800579830917969 0.2240142707009324 0 +255 0.1582184517169784 0.1252836168320661 0 +256 0.1385435679692824 0.03687776649524302 0 +257 0.161604210688566 0.04314183433969498 0 +258 0.1814079862314785 0.1568483656744833 0 +259 0.1828289473684211 0.2486842105263164 0 +260 0.1703036687884799 0.02115535474830363 0 +261 0.1819564642227489 0.1798014391358503 0 +262 0.1631404858437857 0.068445525449567 0 +263 0.1850532992466589 0.2733959534404529 0 +264 0.1658208183126526 0.1886296750753721 0 +265 0.1643546216063115 0.08981986497485502 0 +266 0.1389282654315143 0.0212030012819428 0 +267 0.1293516588741512 0.01711282011379056 0 +268 0.1872233408384842 0.2959417214191695 0 +269 0.1855143321029322 0.01732702153865059 0 +$EndNodes +$Elements +575 +1 1 2 16 1 1 9 +2 1 2 16 1 9 10 +3 1 2 16 1 10 11 +4 1 2 16 1 11 12 +5 1 2 16 1 12 13 +6 1 2 16 1 13 14 +7 1 2 16 1 14 2 +8 1 2 17 2 5 15 +9 1 2 17 2 15 16 +10 1 2 17 2 16 17 +11 1 2 17 2 17 1 +12 1 2 17 3 5 18 +13 1 2 17 3 18 19 +14 1 2 17 3 19 6 +15 1 2 17 4 6 20 +16 1 2 17 4 20 21 +17 1 2 17 4 21 22 +18 1 2 17 4 22 2 +19 1 2 17 5 7 23 +20 1 2 17 5 23 24 +21 1 2 17 5 24 6 +22 1 2 17 6 8 25 +23 1 2 17 6 25 26 +24 1 2 17 6 26 7 +25 1 2 17 7 8 27 +26 1 2 17 7 27 28 +27 1 2 17 7 28 5 +28 1 2 14 8 1 29 +29 1 2 14 8 29 30 +30 1 2 14 8 30 31 +31 1 2 14 8 31 32 +32 1 2 14 8 32 33 +33 1 2 14 8 33 34 +34 1 2 14 8 34 35 +35 1 2 14 8 35 36 +36 1 2 14 8 36 37 +37 1 2 14 8 37 38 +38 1 2 14 8 38 39 +39 1 2 14 8 39 40 +40 1 2 14 8 40 41 +41 1 2 14 8 41 42 +42 1 2 14 8 42 43 +43 1 2 14 8 43 44 +44 1 2 14 8 44 45 +45 1 2 14 8 45 46 +46 1 2 14 8 46 4 +47 1 2 15 9 3 47 +48 1 2 15 9 47 48 +49 1 2 15 9 48 49 +50 1 2 15 9 49 50 +51 1 2 15 9 50 51 +52 1 2 15 9 51 52 +53 1 2 15 9 52 4 +54 1 2 14 10 3 53 +55 1 2 14 10 53 54 +56 1 2 14 10 54 55 +57 1 2 14 10 55 56 +58 1 2 14 10 56 57 +59 1 2 14 10 57 58 +60 1 2 14 10 58 59 +61 1 2 14 10 59 60 +62 1 2 14 10 60 61 +63 1 2 14 10 61 62 +64 1 2 14 10 62 63 +65 1 2 14 10 63 64 +66 1 2 14 10 64 65 +67 1 2 14 10 65 66 +68 1 2 14 10 66 67 +69 1 2 14 10 67 68 +70 1 2 14 10 68 69 +71 1 2 14 10 69 70 +72 1 2 14 10 70 2 +73 1 2 17 11 3 71 +74 1 2 17 11 71 72 +75 1 2 17 11 72 73 +76 1 2 17 11 73 74 +77 1 2 17 11 74 75 +78 1 2 17 11 75 76 +79 1 2 17 11 76 77 +80 1 2 17 11 77 78 +81 1 2 17 11 78 79 +82 1 2 17 11 79 80 +83 1 2 17 11 80 81 +84 1 2 17 11 81 82 +85 1 2 17 11 82 83 +86 1 2 17 11 83 84 +87 1 2 17 11 84 85 +88 1 2 17 11 85 86 +89 1 2 17 11 86 87 +90 1 2 17 11 87 88 +91 1 2 17 11 88 7 +92 2 2 13 1 19 12 6 +93 2 2 13 1 18 5 11 +94 2 2 13 1 17 1 9 +95 2 2 13 1 14 2 22 +96 2 2 13 1 11 5 10 +97 2 2 13 1 10 5 15 +98 2 2 13 1 6 12 13 +99 2 2 13 1 6 13 20 +100 2 2 13 1 10 16 9 +101 2 2 13 1 9 16 17 +102 2 2 13 1 18 11 19 +103 2 2 13 1 21 13 14 +104 2 2 13 1 21 14 22 +105 2 2 13 1 15 16 10 +106 2 2 13 1 11 12 19 +107 2 2 13 1 20 13 21 +108 2 2 12 2 26 91 90 +109 2 2 12 2 90 91 19 +110 2 2 12 2 18 91 89 +111 2 2 12 2 89 91 25 +112 2 2 12 2 25 91 26 +113 2 2 12 2 19 91 18 +114 2 2 12 2 28 89 27 +115 2 2 12 2 23 90 24 +116 2 2 12 2 18 89 5 +117 2 2 12 2 5 89 28 +118 2 2 12 2 6 90 19 +119 2 2 12 2 24 90 6 +120 2 2 12 2 7 90 23 +121 2 2 12 2 26 90 7 +122 2 2 12 2 8 89 25 +123 2 2 12 2 27 89 8 +124 2 2 13 3 193 222 130 +125 2 2 13 3 218 227 125 +126 2 2 13 3 167 206 87 +127 2 2 13 3 8 98 25 +128 2 2 13 3 17 236 161 +129 2 2 13 3 52 213 4 +130 2 2 13 3 103 187 184 +131 2 2 13 3 1 29 17 +132 2 2 13 3 103 225 187 +133 2 2 13 3 71 47 3 +134 2 2 13 3 29 236 17 +135 2 2 13 3 30 211 175 +136 2 2 13 3 163 214 117 +137 2 2 13 3 126 206 167 +138 2 2 13 3 15 227 218 +139 2 2 13 3 71 209 47 +140 2 2 13 3 205 222 206 +141 2 2 13 3 205 206 126 +142 2 2 13 3 88 167 87 +143 2 2 13 3 102 219 186 +144 2 2 13 3 130 202 193 +145 2 2 13 3 107 222 193 +146 2 2 13 3 8 99 98 +147 2 2 13 3 164 173 101 +148 2 2 13 3 119 173 168 +149 2 2 13 3 117 214 122 +150 2 2 13 3 163 211 96 +151 2 2 13 3 199 214 169 +152 2 2 13 3 136 178 137 +153 2 2 13 3 158 211 163 +154 2 2 13 3 137 178 140 +155 2 2 13 3 168 197 174 +156 2 2 13 3 186 199 102 +157 2 2 13 3 7 208 88 +158 2 2 13 3 101 234 164 +159 2 2 13 3 190 230 179 +160 2 2 13 3 15 218 217 +161 2 2 13 3 88 224 167 +162 2 2 13 3 168 173 164 +163 2 2 13 3 197 228 174 +164 2 2 13 3 110 188 73 +165 2 2 13 3 122 214 199 +166 2 2 13 3 73 188 72 +167 2 2 13 3 94 194 166 +168 2 2 13 3 93 204 165 +169 2 2 13 3 218 232 114 +170 2 2 13 3 168 174 119 +171 2 2 13 3 45 213 174 +172 2 2 13 3 159 207 191 +173 2 2 13 3 175 211 158 +174 2 2 13 3 110 190 179 +175 2 2 13 3 176 219 102 +176 2 2 13 3 165 229 93 +177 2 2 13 3 179 230 162 +178 2 2 13 3 32 169 96 +179 2 2 13 3 159 191 165 +180 2 2 13 3 193 202 97 +181 2 2 13 3 200 225 185 +182 2 2 13 3 175 236 30 +183 2 2 13 3 127 171 170 +184 2 2 13 3 46 213 45 +185 2 2 13 3 117 212 163 +186 2 2 13 3 206 222 107 +187 2 2 13 3 169 214 96 +188 2 2 13 3 172 234 101 +189 2 2 13 3 127 170 48 +190 2 2 13 3 203 232 218 +191 2 2 13 3 166 220 94 +192 2 2 13 3 47 209 127 +193 2 2 13 3 32 96 31 +194 2 2 13 3 51 173 119 +195 2 2 13 3 201 202 134 +196 2 2 13 3 102 199 169 +197 2 2 13 3 191 207 92 +198 2 2 13 3 170 172 101 +199 2 2 13 3 185 225 103 +200 2 2 13 3 102 169 33 +201 2 2 13 3 149 198 147 +202 2 2 13 3 47 127 48 +203 2 2 13 3 49 170 101 +204 2 2 13 3 171 172 170 +205 2 2 13 3 34 176 102 +206 2 2 13 3 33 169 32 +207 2 2 13 3 147 198 159 +208 2 2 13 3 96 214 163 +209 2 2 13 3 101 173 50 +210 2 2 13 3 174 228 45 +211 2 2 13 3 51 119 52 +212 2 2 13 3 127 188 171 +213 2 2 13 3 48 170 49 +214 2 2 13 3 105 228 197 +215 2 2 13 3 49 101 50 +216 2 2 13 3 155 162 152 +217 2 2 13 3 97 202 201 +218 2 2 13 3 36 177 104 +219 2 2 13 3 34 102 33 +220 2 2 13 3 50 173 51 +221 2 2 13 3 104 176 35 +222 2 2 13 3 17 161 16 +223 2 2 13 3 203 218 125 +224 2 2 13 3 36 104 35 +225 2 2 13 3 35 176 34 +226 2 2 13 3 159 165 146 +227 2 2 13 3 95 178 177 +228 2 2 13 3 180 181 138 +229 2 2 13 3 135 182 134 +230 2 2 13 3 152 198 149 +231 2 2 13 3 139 180 138 +232 2 2 13 3 38 184 95 +233 2 2 13 3 95 177 37 +234 2 2 13 3 155 179 162 +235 2 2 13 3 40 185 103 +236 2 2 13 3 136 210 178 +237 2 2 13 3 138 182 135 +238 2 2 13 3 144 160 142 +239 2 2 13 3 37 177 36 +240 2 2 13 3 181 183 182 +241 2 2 13 3 142 180 139 +242 2 2 13 3 103 184 39 +243 2 2 13 3 181 182 138 +244 2 2 13 3 148 156 150 +245 2 2 13 3 94 185 41 +246 2 2 13 3 38 95 37 +247 2 2 13 3 160 180 142 +248 2 2 13 3 143 144 142 +249 2 2 13 3 40 103 39 +250 2 2 13 3 42 94 41 +251 2 2 13 3 141 142 139 +252 2 2 13 3 145 146 144 +253 2 2 13 3 39 184 38 +254 2 2 13 3 41 185 40 +255 2 2 13 3 178 215 140 +256 2 2 13 3 137 138 135 +257 2 2 13 3 141 143 142 +258 2 2 13 3 143 145 144 +259 2 2 13 3 154 197 168 +260 2 2 13 3 86 107 85 +261 2 2 13 3 147 159 146 +262 2 2 13 3 137 139 138 +263 2 2 13 3 185 220 200 +264 2 2 13 3 85 193 84 +265 2 2 13 3 42 194 94 +266 2 2 13 3 76 92 75 +267 2 2 13 3 148 195 156 +268 2 2 13 3 79 93 78 +269 2 2 13 3 80 192 79 +270 2 2 13 3 81 106 80 +271 2 2 13 3 77 191 76 +272 2 2 13 3 84 97 83 +273 2 2 13 3 156 166 150 +274 2 2 13 3 140 141 139 +275 2 2 13 3 107 193 85 +276 2 2 13 3 43 194 42 +277 2 2 13 3 141 157 143 +278 2 2 13 3 182 201 134 +279 2 2 13 3 136 137 135 +280 2 2 13 3 133 135 134 +281 2 2 13 3 84 193 97 +282 2 2 13 3 140 187 141 +283 2 2 13 3 79 192 93 +284 2 2 13 3 137 140 139 +285 2 2 13 3 76 191 92 +286 2 2 13 3 106 192 80 +287 2 2 13 3 141 187 157 +288 2 2 13 3 133 136 135 +289 2 2 13 3 133 134 131 +290 2 2 13 3 153 155 152 +291 2 2 13 3 145 147 146 +292 2 2 13 3 44 105 43 +293 2 2 13 3 162 198 152 +294 2 2 13 3 132 189 133 +295 2 2 13 3 133 189 136 +296 2 2 13 3 143 195 145 +297 2 2 13 3 151 152 149 +298 2 2 13 3 87 206 86 +299 2 2 13 3 157 195 143 +300 2 2 13 3 105 194 43 +301 2 2 13 3 145 148 147 +302 2 2 13 3 75 190 74 +303 2 2 13 3 187 215 184 +304 2 2 13 3 145 195 148 +305 2 2 13 3 132 133 131 +306 2 2 13 3 92 190 75 +307 2 2 13 3 129 132 131 +308 2 2 13 3 148 150 149 +309 2 2 13 3 119 213 52 +310 2 2 13 3 166 196 150 +311 2 2 13 3 194 221 166 +312 2 2 13 3 151 154 153 +313 2 2 13 3 165 204 146 +314 2 2 13 3 154 168 164 +315 2 2 13 3 144 204 160 +316 2 2 13 3 31 211 30 +317 2 2 13 3 151 196 154 +318 2 2 13 3 183 240 82 +319 2 2 13 3 74 110 73 +320 2 2 13 3 148 149 147 +321 2 2 13 3 113 203 125 +322 2 2 13 3 151 153 152 +323 2 2 13 3 129 186 132 +324 2 2 13 3 196 197 154 +325 2 2 13 3 129 131 130 +326 2 2 13 3 154 164 153 +327 2 2 13 3 210 239 104 +328 2 2 13 3 150 151 149 +329 2 2 13 3 134 202 131 +330 2 2 13 3 150 196 151 +331 2 2 13 3 123 186 129 +332 2 2 13 3 129 130 128 +333 2 2 13 3 189 239 210 +334 2 2 13 3 124 231 224 +335 2 2 13 3 74 190 110 +336 2 2 13 3 16 227 15 +337 2 2 13 3 183 201 182 +338 2 2 13 3 98 223 25 +339 2 2 13 3 181 216 106 +340 2 2 13 3 195 200 156 +341 2 2 13 3 162 207 198 +342 2 2 13 3 180 216 181 +343 2 2 13 3 157 200 195 +344 2 2 13 3 123 199 186 +345 2 2 13 3 156 220 166 +346 2 2 13 3 131 202 130 +347 2 2 13 3 224 231 167 +348 2 2 13 3 123 128 118 +349 2 2 13 3 122 199 123 +350 2 2 13 3 123 129 128 +351 2 2 13 3 132 219 189 +352 2 2 13 3 146 204 144 +353 2 2 13 3 110 237 188 +354 2 2 13 3 130 222 128 +355 2 2 13 3 113 125 121 +356 2 2 13 3 155 234 172 +357 2 2 13 3 177 210 104 +358 2 2 13 3 122 123 118 +359 2 2 13 3 97 226 83 +360 2 2 13 3 106 240 181 +361 2 2 13 3 82 240 81 +362 2 2 13 3 93 229 78 +363 2 2 13 3 77 229 191 +364 2 2 13 3 181 240 183 +365 2 2 13 3 191 229 165 +366 2 2 13 3 187 225 157 +367 2 2 13 3 82 226 183 +368 2 2 13 3 45 228 44 +369 2 2 13 3 96 211 31 +370 2 2 13 3 86 206 107 +371 2 2 13 3 126 167 115 +372 2 2 13 3 117 122 118 +373 2 2 13 3 158 163 120 +374 2 2 13 3 178 210 177 +375 2 2 13 3 184 215 95 +376 2 2 13 3 124 223 98 +377 2 2 13 3 174 213 119 +378 2 2 13 3 188 209 72 +379 2 2 13 3 201 226 97 +380 2 2 13 3 112 117 116 +381 2 2 13 3 95 215 178 +382 2 2 13 3 128 205 118 +383 2 2 13 3 166 221 196 +384 2 2 13 3 30 236 29 +385 2 2 13 3 94 220 185 +386 2 2 13 3 153 234 155 +387 2 2 13 3 116 126 115 +388 2 2 13 3 99 203 108 +389 2 2 13 3 117 118 116 +390 2 2 13 3 121 175 158 +391 2 2 13 3 198 207 159 +392 2 2 13 3 108 111 109 +393 2 2 13 3 116 205 126 +394 2 2 13 3 108 113 111 +395 2 2 13 3 118 205 116 +396 2 2 13 3 121 158 120 +397 2 2 13 3 189 210 136 +398 2 2 13 3 172 233 155 +399 2 2 13 3 72 209 71 +400 2 2 13 3 160 216 180 +401 2 2 13 3 112 116 115 +402 2 2 13 3 167 231 115 +403 2 2 13 3 113 121 120 +404 2 2 13 3 99 108 100 +405 2 2 13 3 140 215 187 +406 2 2 13 3 106 216 192 +407 2 2 13 3 112 115 109 +408 2 2 13 3 163 212 120 +409 2 2 13 3 111 112 109 +410 2 2 13 3 162 230 207 +411 2 2 13 3 197 221 105 +412 2 2 13 3 26 208 7 +413 2 2 13 3 161 227 16 +414 2 2 13 3 186 219 132 +415 2 2 13 3 113 120 111 +416 2 2 13 3 196 221 197 +417 2 2 13 3 127 209 188 +418 2 2 13 3 108 109 100 +419 2 2 13 3 171 233 172 +420 2 2 13 3 105 221 194 +421 2 2 13 3 99 100 98 +422 2 2 13 3 125 235 121 +423 2 2 13 3 208 224 88 +424 2 2 13 3 25 223 26 +425 2 2 13 3 216 238 192 +426 2 2 13 3 83 226 82 +427 2 2 13 3 179 237 110 +428 2 2 13 3 8 232 99 +429 2 2 13 3 175 235 161 +430 2 2 13 3 108 203 113 +431 2 2 13 3 160 238 216 +432 2 2 13 3 233 237 179 +433 2 2 13 3 204 238 160 +434 2 2 13 3 78 229 77 +435 2 2 13 3 104 239 176 +436 2 2 13 3 219 239 189 +437 2 2 13 3 200 220 156 +438 2 2 13 3 112 212 117 +439 2 2 13 3 109 231 100 +440 2 2 13 3 44 228 105 +441 2 2 13 3 183 226 201 +442 2 2 13 3 188 237 171 +443 2 2 13 3 157 225 200 +444 2 2 13 3 161 236 175 +445 2 2 13 3 81 240 106 +446 2 2 13 3 92 230 190 +447 2 2 13 3 128 222 205 +448 2 2 13 3 192 238 93 +449 2 2 13 3 171 237 233 +450 2 2 13 3 155 233 179 +451 2 2 13 3 120 212 111 +452 2 2 13 3 227 235 125 +453 2 2 13 3 26 223 208 +454 2 2 13 3 27 114 8 +455 2 2 13 3 164 234 153 +456 2 2 13 3 111 212 112 +457 2 2 13 3 99 232 203 +458 2 2 13 3 100 124 98 +459 2 2 13 3 28 217 27 +460 2 2 13 3 93 238 204 +461 2 2 13 3 100 231 124 +462 2 2 13 3 121 235 175 +463 2 2 13 3 115 231 109 +464 2 2 13 3 207 230 92 +465 2 2 13 3 124 224 223 +466 2 2 13 3 223 224 208 +467 2 2 13 3 217 218 114 +468 2 2 13 3 176 239 219 +469 2 2 13 3 161 235 227 +470 2 2 13 3 27 217 114 +471 2 2 13 3 114 232 8 +472 2 2 13 3 4 213 46 +473 2 2 13 3 15 217 5 +474 2 2 13 3 5 217 28 +475 2 2 13 4 80 259 79 +476 2 2 13 4 249 254 81 +477 2 2 13 4 254 259 80 +478 2 2 13 4 81 254 80 +479 2 2 13 4 24 267 23 +480 2 2 13 4 56 74 75 +481 2 2 13 4 57 76 58 +482 2 2 13 4 56 75 57 +483 2 2 13 4 245 260 257 +484 2 2 13 4 251 255 243 +485 2 2 13 4 57 75 76 +486 2 2 13 4 79 259 78 +487 2 2 13 4 259 263 78 +488 2 2 13 4 246 255 251 +489 2 2 13 4 245 257 256 +490 2 2 13 4 6 267 24 +491 2 2 13 4 242 266 256 +492 2 2 13 4 257 260 244 +493 2 2 13 4 256 266 245 +494 2 2 13 4 78 263 77 +495 2 2 13 4 263 268 77 +496 2 2 13 4 76 268 58 +497 2 2 13 4 246 258 253 +498 2 2 13 4 64 261 258 +499 2 2 13 4 65 258 246 +500 2 2 13 4 21 245 20 +501 2 2 13 4 250 257 244 +502 2 2 13 4 247 250 68 +503 2 2 13 4 69 244 70 +504 2 2 13 4 69 250 244 +505 2 2 13 4 67 251 247 +506 2 2 13 4 62 249 63 +507 2 2 13 4 253 255 246 +508 2 2 13 4 246 251 66 +509 2 2 13 4 65 246 66 +510 2 2 13 4 67 247 68 +511 2 2 13 4 68 250 69 +512 2 2 13 4 81 264 249 +513 2 2 13 4 66 251 67 +514 2 2 13 4 61 259 254 +515 2 2 13 4 88 241 87 +516 2 2 13 4 87 252 86 +517 2 2 13 4 61 254 62 +518 2 2 13 4 86 243 85 +519 2 2 13 4 84 253 83 +520 2 2 13 4 258 261 248 +521 2 2 13 4 83 248 82 +522 2 2 13 4 241 252 87 +523 2 2 13 4 86 252 243 +524 2 2 13 4 62 254 249 +525 2 2 13 4 83 253 248 +526 2 2 13 4 250 262 257 +527 2 2 13 4 85 255 84 +528 2 2 13 4 243 255 85 +529 2 2 13 4 243 265 251 +530 2 2 13 4 84 255 253 +531 2 2 13 4 77 268 76 +532 2 2 13 4 70 269 2 +533 2 2 13 4 2 269 22 +534 2 2 13 4 247 262 250 +535 2 2 13 4 63 261 64 +536 2 2 13 4 64 258 65 +537 2 2 13 4 22 260 21 +538 2 2 13 4 88 256 241 +539 2 2 13 4 245 266 20 +540 2 2 13 4 21 260 245 +541 2 2 13 4 256 257 241 +542 2 2 13 4 60 259 61 +543 2 2 13 4 253 258 248 +544 2 2 13 4 257 262 241 +545 2 2 13 4 251 265 247 +546 2 2 13 4 23 267 242 +547 2 2 13 4 60 263 259 +548 2 2 13 4 247 265 262 +549 2 2 13 4 82 264 81 +550 2 2 13 4 249 261 63 +551 2 2 13 4 244 269 70 +552 2 2 13 4 241 262 252 +553 2 2 13 4 59 263 60 +554 2 2 13 4 59 268 263 +555 2 2 13 4 248 264 82 +556 2 2 13 4 262 265 252 +557 2 2 13 4 58 268 59 +558 2 2 13 4 252 265 243 +559 2 2 13 4 22 269 260 +560 2 2 13 4 261 264 248 +561 2 2 13 4 260 269 244 +562 2 2 13 4 249 264 261 +563 2 2 13 4 242 267 266 +564 2 2 13 4 266 267 20 +565 2 2 13 4 23 242 7 +566 2 2 13 4 20 267 6 +567 2 2 13 4 53 3 71 +568 2 2 13 4 7 256 88 +569 2 2 13 4 242 256 7 +570 2 2 13 4 53 72 54 +571 2 2 13 4 71 72 53 +572 2 2 13 4 54 73 55 +573 2 2 13 4 72 73 54 +574 2 2 13 4 55 74 56 +575 2 2 13 4 55 73 74 +$EndElements diff --git a/examples/TFM/bc_dir/bc_config b/examples/TFM/bc_dir/bc_config new file mode 100644 index 0000000..921516e --- /dev/null +++ b/examples/TFM/bc_dir/bc_config @@ -0,0 +1,13 @@ +[DIRICHLET] +u_c = wall -> [0.0, 0.0] + bottom -> [0.0, 0.0] + +[SLIP] +u_d = wall|bottom + +[ZERO_STRESS] +u_c = outlet +u_d = outlet + +[ZERO_BACKFLOW] +alpha_c = outlet|wall|bottom diff --git a/examples/TFM/config b/examples/TFM/config new file mode 100644 index 0000000..9e0a650 --- /dev/null +++ b/examples/TFM/config @@ -0,0 +1,51 @@ +[MESH] +filename = bc_2d.msh + +[FINITE ELEMENT SPACE] +elements = u_c -> HDiv + u_d -> HDiv + p -> L2 + alpha_c -> L2 +interpolant_order = 3 + +[DG] +DG = True +interior_penalty_coefficient = 10.0 + +[SOLVER] +linear_solver = direct +preconditioner = default +linearization_method = Picard +nonlinear_solver = default +nonlinear_tolerance = relative -> 1e-5 + absolute -> 1e-5 +nonlinear_max_iterations = 10 +relaxation_factors = 1.0, 1.0, 1.0, 1.0 + +[TRANSIENT] +transient = True +scheme = implicit euler +time_range = 0.0, 3.0 +dt = 2e-3 + +[VISUALIZATION] +save_to_file = True +save_type = .vtu +save_vtu_each_timestep = True +save_frequency = 5, numit + +[OTHER] +num_threads = 6 +model = TwoFluidModel +run_dir = . +slope_limiter = True +diffusion_switch = True + +[TFM] +canonical_form = C-TFM +IME = drag -> Tomiyama + lift -> Tomiyama + virtual_mass -> ConstantCoefficient + laminar_dispersion -> ConstantCoefficient +lift_wall_deactivation = True +lift_wall_boundaries = wall|bottom diff --git a/examples/TFM/ic_dir/ic_config b/examples/TFM/ic_dir/ic_config new file mode 100644 index 0000000..80b608f --- /dev/null +++ b/examples/TFM/ic_dir/ic_config @@ -0,0 +1,7 @@ +[TwoFluidModel] +# Initially quiescent liquid column. +u_c = all -> [0.0, 0.0] +u_d = all -> [0.0, 0.0] + +# Column initially contains only the continuous phase. +alpha_c = all -> 1.0 diff --git a/examples/TFM/model_dir/model_config b/examples/TFM/model_dir/model_config new file mode 100644 index 0000000..cd05652 --- /dev/null +++ b/examples/TFM/model_dir/model_config @@ -0,0 +1,24 @@ +[PARAMETERS] +# Phase densities (kg/m^3) +rho_c = all -> 998.2 +rho_d = all -> 1.204 + +# Kinematic viscosities (m^2/s) +nu_c = all -> 1.0038e-6 +nu_d = all -> 1.5158e-5 + +# Surface tension (N/m) +sigma_c = all -> 0.072 + +# Bubble diameter (m) +dp = all -> 3e-3 + +# Interphase force constants +c_vm = all -> 0.5 +cdis = all -> 4.545 +d_artificial = all -> 1e-5 + +[INJECTION] +region = injection +mass_flow_rate = 0.3 +velocity = 0.2 diff --git a/examples/TFM/ref_sol_dir/ref_sol_config b/examples/TFM/ref_sol_dir/ref_sol_config new file mode 100644 index 0000000..b0ece3e --- /dev/null +++ b/examples/TFM/ref_sol_dir/ref_sol_config @@ -0,0 +1,3 @@ +[REFERENCE SOLUTIONS] + +[METRICS] diff --git a/opencmp/config_functions/boundary_conditions.py b/opencmp/config_functions/boundary_conditions.py index 9d55ff4..3c886f5 100644 --- a/opencmp/config_functions/boundary_conditions.py +++ b/opencmp/config_functions/boundary_conditions.py @@ -33,10 +33,34 @@ def __init__(self, config_rel_path: str, import_dir: str, mesh: Mesh, bc_types: super().__init__(config_rel_path, import_dir, mesh, t_param) # Load the BC dict from the BC configfile. + # These types are fully determined by their name, so they are marker-only + # (``u_d = wall|bottom``) instead of the parser's ``marker -> value``. + marker_only_types = ({'SLIP', 'ZERO_STRESS', 'ZERO_GRADIENT', 'ZERO_BACKFLOW'} + & set(bc_types)) self.bc_dict, self.bc_re_parse_dict = self.config.get_three_level_dict(self.import_dir, None, self.t_param, new_variables, - white_list=bc_types, - ignore=['VERTICES', 'CENTROIDS']) + white_list=[bc_type for bc_type in bc_types + if bc_type not in marker_only_types], + ignore=['VERTICES', 'CENTROIDS'] + + list(marker_only_types)) + for bc_type in marker_only_types: + marker_dict = {} + section = self.config[bc_type] if self.config.has_section(bc_type) else {} + for var, marker_expression in section.items(): + markers = [marker.strip() for marker in marker_expression.split('|') if marker.strip()] + if not markers: + raise ValueError( + "Boundary condition '[{}] {}' must specify at least one mesh marker." + .format(bc_type, var) + ) + if len(markers) != len(set(markers)): + raise ValueError( + "Boundary condition '[{}] {}' contains duplicate mesh markers." + .format(bc_type, var) + ) + marker_dict[var] = {marker: [] for marker in markers} + self.bc_dict[bc_type.lower()] = marker_dict + self.bc_re_parse_dict[bc_type.lower()] = {var: {} for var in marker_dict} # Used to keep track of which unknown variables the user has already been warned about # A variable present in the config file that is not in the current model is not necessarily a bug. diff --git a/opencmp/helpers/dg.py b/opencmp/helpers/dg.py index 3ad41c6..8f7e9fa 100644 --- a/opencmp/helpers/dg.py +++ b/opencmp/helpers/dg.py @@ -102,6 +102,28 @@ def weighted_grad_avg(q: CoefficientFunction, c: CoefficientFunction) -> Coeffic return 0.5 * (c * Grad(q) + c.Other() * Grad(q).Other()) +def weighted_trans_grad_avg(q: CoefficientFunction, c: CoefficientFunction) -> CoefficientFunction: + """ + Returns the average of the transposed gradient of a field weighted by a (possibly discontinuous) coefficient. + + Note this is NOT weighted_grad_avg(q, c).trans in general: this weights on the right, + (Grad(q) * c)^T, and the two only coincide for scalar c. + + Args: + q: The field. + c: The coefficient weighting the gradient on each side of the facet. + + Returns: + The average of (Grad(q) * c)^T 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 * ((Grad(q) * c).trans + (Grad(q.Other()) * c.Other()).trans) + else: + return 0.5 * ((Grad(q) * c).trans + (Grad(q).Other() * c.Other()).trans) + + def weighted_div_avg(q: CoefficientFunction, c: CoefficientFunction) -> CoefficientFunction: """ Returns the average of the divergence of a field weighted by a (possibly discontinuous) coefficient. diff --git a/opencmp/helpers/ngsolve_.py b/opencmp/helpers/ngsolve_.py index 1838d14..de16f33 100644 --- a/opencmp/helpers/ngsolve_.py +++ b/opencmp/helpers/ngsolve_.py @@ -45,6 +45,21 @@ def construct_identity_mat(dim: int) -> CoefficientFunction: return ngs.CoefficientFunction(tuple(lst), dims=(dim, dim)) +def curl_3d(gfu: GridFunction) -> CoefficientFunction: + """ + Returns the curl of a 3D vector field computed from its gradient. + + Args: + gfu: The 3D vector field. + + Returns: + The curl of gfu. + """ + + g = ngs.grad(gfu) + return ngs.CoefficientFunction((g[5] - g[7], -(g[2] - g[6]), g[1] - g[3])) + + def get_special_functions(mesh: Mesh, nu: float) \ -> Tuple[CoefficientFunction, CoefficientFunction, CoefficientFunction, CoefficientFunction]: """ diff --git a/opencmp/helpers/wall_func.py b/opencmp/helpers/wall_func.py index afe2416..9d775b5 100644 --- a/opencmp/helpers/wall_func.py +++ b/opencmp/helpers/wall_func.py @@ -22,6 +22,48 @@ import ngsolve as ngs +def wall_distance(mesh: ngs.comp.Mesh, wall_boundary: str = 'wall', + order: int = 2, relax: float = 0.1) -> ngs.GridFunction: + """ + Distance to ``wall_boundary`` from a regularized Eikonal solve. + + Free of any turbulence modelling, so models that need only a wall distance + (e.g. the two-fluid lift wall-deactivation taper) can use it without + constructing a :class:`KEpsilonWallFunction`. + + Args: + mesh: The mesh used for the simulation. + wall_boundary: Boundary marker(s) to measure the distance from. + order: Order of the H1 space the distance field lives in. + relax: Regularization, as a multiple of the local mesh size. + + Returns: + The wall-distance field. + """ + + eps = relax * ngs.specialcf.mesh_size + fes = ngs.H1(mesh, order=order, dirichlet=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 + + class KEpsilonWallFunction: """Wall-layer eddy viscosity and dissipation for high-Re k-epsilon. @@ -180,27 +222,7 @@ def _find_physical_wall_facets(self): 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 + return wall_distance(self.mesh, self.wall_boundary, order, relax) # ------------------------------------------------------------------ # Per-iteration update diff --git a/opencmp/models/__init__.py b/opencmp/models/__init__.py index f3999b2..a67dba1 100644 --- a/opencmp/models/__init__.py +++ b/opencmp/models/__init__.py @@ -27,6 +27,7 @@ from .stokes_dim import StokesDIM from .multi_component_ins import MultiComponentINS from .k_epsilon import KEpsilonINS +from .tfm import TwoFluidModel models_dict = {"INS": INS, @@ -36,7 +37,8 @@ "Stokes": Stokes, "Stokes-DIM": StokesDIM, "MultiComponentINS": MultiComponentINS, - "KEpsilonINS": KEpsilonINS} + "KEpsilonINS": KEpsilonINS, + "TwoFluidModel": TwoFluidModel} # Helper functions from .misc import get_model_class diff --git a/opencmp/models/tfm.py b/opencmp/models/tfm.py new file mode 100644 index 0000000..1cd6204 --- /dev/null +++ b/opencmp/models/tfm.py @@ -0,0 +1,969 @@ +######################################################################################################################## +# 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 # +# . # +######################################################################################################################## + +import logging +from typing import Dict, List, Optional, Union + +import ngsolve as ngs +from ngsolve import GridFunction, FESpace, BilinearForm, LinearForm, Preconditioner, Parameter +from ngsolve.comp import ProxyFunction + +from . import Model +from ..helpers.dg import avg, jump, grad_avg, weighted_grad_avg, weighted_div_avg +from ..helpers.math import tanh, Max, Min +from ..helpers.ngsolve_ import get_special_functions, curl_3d +from ..helpers.limiter import Limiter +from ..helpers.error import norm, mean +from ..helpers.wall_func import wall_distance + + +class TwoFluidModel(Model): + """ + Laminar Euler-Euler two-fluid model for two-phase pipe flows. + + Implements four coupled equations: u_c (HDiv/BDM), u_d (HDiv/BDM), + p (L2), alpha_c (L2). The nonlinear system is resolved by Picard + (fixed-point) iteration each time step. Turbulence closure is left to a + subclass (cf. INS / KEpsilonINS). + """ + + # Interphase momentum exchange mechanisms and the closure models available to each. + IME_MODELS = {'drag': ('Tomiyama', 'SchillerNaumann'), + 'lift': ('Tomiyama', 'LegendreMagnaudet'), + 'virtual_mass': ('ConstantCoefficient',), + 'laminar_dispersion': ('ConstantCoefficient',)} + + # ------------------------------------------------------------------ + # Abstract-method overrides — bookkeeping + # ------------------------------------------------------------------ + + @staticmethod + def allows_explicit_schemes() -> bool: + return False + + def _define_model_components(self) -> Dict[str, Optional[int]]: + return {'u_c': 0, 'u_d': 1, 'p': 2, 'alpha_c': 3} + + def _define_model_local_error_components(self) -> Dict[str, bool]: + return {'u_c': True, 'u_d': True, 'p': False, 'alpha_c': True} + + def _define_time_derivative_components(self) -> List[Dict[str, bool]]: + return [{'u_c': True, 'u_d': True, 'p': False, 'alpha_c': True}] + + def _define_num_weak_forms(self) -> int: + return 1 + + def _define_bc_types(self) -> List[str]: + return ['dirichlet', 'zero_stress', 'zero_gradient', 'zero_backflow', 'slip'] + + # ------------------------------------------------------------------ + # Phase 1 — lifecycle hooks + # ------------------------------------------------------------------ + + def _tfm_option(self, key: str, val_type, default): + """Read an optional model-specific setting from [TFM].""" + if not self.config.has_option('TFM', key): + return default + try: + return self.config.get_item(['TFM', key], val_type, quiet=True) + except (TypeError, ValueError) as exc: + raise ValueError("Invalid [TFM] value for '{}'.".format(key)) from exc + + def _other_option(self, key: str, val_type, default): + """Read an optional generic numerical setting from [OTHER].""" + if not self.config.has_option('OTHER', key): + return default + return self.config.get_item(['OTHER', key], val_type, quiet=True) + + def _parse_ime(self) -> Dict[str, str]: + """Read [TFM] IME, which maps each active mechanism to its closure model.""" + if not self.config.has_option('TFM', 'ime'): + return {'drag': 'Tomiyama'} + if not self.config['TFM']['ime'].strip(): + return {} + + ime = {key.lower(): value for key, value in + self.config.get_dict(['TFM', 'ime'], '', None, all_str=True).items()} + + unknown_ime = set(ime) - set(self.IME_MODELS) + if unknown_ime: + raise ValueError('Unknown [TFM] IME mechanism(s): {}.'.format(', '.join(sorted(unknown_ime)))) + for mechanism, model in ime.items(): + if model not in self.IME_MODELS[mechanism]: + raise ValueError("[TFM] IME '{}' must use one of: {}.".format( + mechanism, ', '.join(self.IME_MODELS[mechanism]))) + if 'laminar_dispersion' in ime and 'drag' not in ime: + raise ValueError("[TFM] IME 'laminar_dispersion' requires 'drag'.") + + return ime + + def _pre_init(self) -> None: + allowed_keys = {'canonical_form', 'ime', 'lift_wall_deactivation', 'lift_wall_boundaries'} + unknown_keys = set(self.config['TFM']) - allowed_keys if self.config.has_section('TFM') else set() + if unknown_keys: + raise ValueError('Unknown [TFM] option(s): {}.'.format(', '.join(sorted(unknown_keys)))) + + self.canonical_form = self._tfm_option('canonical_form', str, 'B-TFM') + if self.canonical_form not in ('B-TFM', 'C-TFM'): + raise ValueError("[TFM] canonical_form must be 'B-TFM' or 'C-TFM'.") + + self.slope_limiter = self._other_option('slope_limiter', bool, True) + + ime = self._parse_ime() + self.drag_switch = 'drag' in ime + self.VM_switch = 'virtual_mass' in ime + self.Disp_switch = 'laminar_dispersion' in ime + self.Lift_switch = 'lift' in ime + + self.drag_model = ime.get('drag', 'Tomiyama') + self.lift_model = ime.get('lift', 'Tomiyama') + + # Artificial diffusion on the alpha_c transport. + self.diffusion_switch = self._other_option('diffusion_switch', bool, False) + self.mean_zero_pressure = self._other_option('mean_zero_pressure', bool, False) + + # Lift wall-deactivation taper (off by default; costs one wall-distance solve). + self.lift_wall_deactivation = self._tfm_option('lift_wall_deactivation', bool, False) + self.lift_wall_boundaries = self._tfm_option('lift_wall_boundaries', str, '') + if self.lift_wall_deactivation and not self.Lift_switch: + raise ValueError('lift_wall_deactivation requires lift in [TFM] IME.') + + # Loaded from model_dir/model_config after ModelFunctions is initialized. + self.injection_switch = False + self.injection_region = '' + self.inj_mass_flowrate = 0.0 + self.inj_velocity = 0.0 + + def _velocity_space_boundaries(self, var: str) -> str: + """Combine full Dirichlet and normal-only slip constraints for H(div).""" + dirichlet = list(self.BC.get('dirichlet', {}).get(var, {})) + slip = list(self.BC.get('slip', {}).get(var, {})) + overlap = set(dirichlet) & set(slip) + if overlap: + raise ValueError( + "Boundary marker(s) {} cannot be both DIRICHLET and SLIP for '{}'." + .format(', '.join(sorted(overlap)), var) + ) + return '|'.join(dirichlet + slip) + + def _construct_fes(self) -> FESpace: + # The weak form assumes H(div)-conforming velocities (the mixture mass + # conservation and the UDS momentum fluxes are written on u.n) and a + # discontinuous pressure, so those two are constrained rather than free. + for name in ('u_c', 'u_d'): + if self.element[name] != 'HDiv': + raise ValueError("TwoFluidModel requires an HDiv element for '{}', got '{}'." + .format(name, self.element[name])) + if self.element['p'] != 'L2': + raise ValueError("TwoFluidModel requires an L2 element for 'p', got '{}'." + .format(self.element['p'])) + + scalar_ord = max(self.interp_ord - 1, 0) + # HDiv strongly constrains only the normal trace. DIRICHLET markers also + # receive the weak full-vector terms below; SLIP markers do not, leaving + # zero tangential traction as their natural boundary condition. + uc_dirichlet = self._velocity_space_boundaries('u_c') + ud_dirichlet = self._velocity_space_boundaries('u_d') + fes_uc = ngs.HDiv(self.mesh, order=self.interp_ord, + dirichlet=uc_dirichlet, dgjumps=self.DG) + fes_ud = ngs.HDiv(self.mesh, order=self.interp_ord, + dirichlet=ud_dirichlet, dgjumps=self.DG) + fes_p = ngs.L2(self.mesh, order=scalar_ord, dgjumps=self.DG) + # alpha_c is an ordinary scalar transport variable -- any element the user asks for. + # No `dirichlet` here: its BCs are imposed weakly through the UDS flux terms. + fes_ac = getattr(ngs, self.element['alpha_c'])(self.mesh, order=scalar_ord, dgjumps=self.DG) + spaces = [fes_uc, fes_ud, fes_p, fes_ac] + if getattr(self, 'mean_zero_pressure', False): + spaces.append(ngs.NumberSpace(self.mesh)) + return FESpace(spaces, dgjumps=self.DG) + + def _set_model_parameters(self) -> None: + p = self.model_functions.model_parameters_dict + self.rho_c = p['rho_c']['all'] + self.rho_d = p['rho_d']['all'] + self.nu_c = p['nu_c']['all'] + self.nu_d = p['nu_d']['all'] + self.sigma_c = p['sigma_c']['all'] + self.dp = p['dp']['all'] + self.C_VM = p['c_vm']['all'] + self.Cdis = p['cdis']['all'] + self.D_art = p['d_artificial']['all'] + self.f = self.model_functions.model_functions_dict.get('source', {}) + + model_config = self.model_functions.config + self.injection_switch = model_config.has_section('INJECTION') + if self.injection_switch: + allowed = {'region', 'mass_flow_rate', 'velocity'} + unknown = set(model_config['INJECTION']) - allowed + if unknown: + raise ValueError('Unknown model [INJECTION] option(s): {}.' + .format(', '.join(sorted(unknown)))) + for required in allowed: + if not model_config.has_option('INJECTION', required): + raise ValueError("model_dir/model_config [INJECTION] requires '{}'." + .format(required)) + self.injection_region = model_config.get_item(['INJECTION', 'region'], str) + self.inj_mass_flowrate = model_config.get_item(['INJECTION', 'mass_flow_rate'], float) + self.inj_velocity = model_config.get_item(['INJECTION', 'velocity'], float) + + if self.mesh.dim == 2: + self.gravity = ngs.CoefficientFunction((0.0, -9.81)) + else: + self.gravity = ngs.CoefficientFunction((0.0, -9.81, 0.0)) + + def _post_init(self) -> None: + self._validate_bc_variables() + self.nonlinear = True + self.linearize = 'Picard' + + nonlinear_tolerance = self.config.get_dict( + ['SOLVER', 'nonlinear_tolerance'], self.run_dir, None) + self.abs_nonlinear_tolerance = nonlinear_tolerance['absolute'] + self.rel_nonlinear_tolerance = nonlinear_tolerance['relative'] + self.nonlinear_max_iters = self.config.get_item( + ['SOLVER', 'nonlinear_max_iterations'], int) + if self.nonlinear_max_iters < 1: + raise ValueError('nonlinear_max_iterations must be >= 1.') + + try: + relax = self.config.get_list(['SOLVER', 'relaxation_factors'], float) + except Exception: + relax = [] + n_comp = len(self.model_components) + self.relax_factors = relax if len(relax) == n_comp else [1.0] * n_comp + + # Picard iterate (updated at the start of every Picard sub-iteration) + self.UIter = ngs.GridFunction(self.fes) + self.UIter.vec.data = self.IC.vec + + # Piecewise-constant closure GFs (live outside self.fes) + _cl_fes = ngs.L2(self.mesh, order=0) + self.Cd_gfu = ngs.GridFunction(_cl_fes) + self.Cl_gfu = ngs.GridFunction(_cl_fes) + + # Special DG functions (self.nu = ipc * interp_ord^2 is the penalty parameter). + self.n, self.h, self.penalty_interior, self.IM = get_special_functions(self.mesh, self.nu) + + # Slope limiter instance + self._lim = Limiter(self.mesh) + + # Wall distance for the lift taper. Solved once (Eikonal Newton solve), and + # only when the taper is actually on. + if self.lift_wall_deactivation: + if not self.lift_wall_boundaries: + raise ValueError( + 'lift_wall_deactivation requires lift_wall_boundaries to be specified.' + ) + self._wall_dist = wall_distance(self.mesh, self.lift_wall_boundaries) + else: + self._wall_dist = None + + # ------------------------------------------------------------------ + # Internal helpers — closure coefficients + # ------------------------------------------------------------------ + + def _bc_regex(self, bc_type: str, varname: str) -> str: + return '|'.join(self.BC.get(bc_type, {}).get(varname, {}).keys()) + + # Each outflow condition applies to exactly one kind of variable. + BC_VARIABLES = {'zero_stress': ('u_c', 'u_d'), + 'zero_gradient': ('alpha_c',), + 'zero_backflow': ('alpha_c',)} + + def _validate_bc_variables(self) -> None: + for bc_type, allowed in self.BC_VARIABLES.items(): + for var in self.BC.get(bc_type, {}): + if var not in allowed: + raise ValueError( + "[{}] is only meaningful for {}, not '{}'." + .format(bc_type.upper(), ' and '.join(allowed), var)) + + overlap = set(self.BC.get('zero_gradient', {}).get('alpha_c', {})) \ + & set(self.BC.get('zero_backflow', {}).get('alpha_c', {})) + if overlap: + raise ValueError( + "Boundary marker(s) {} cannot be both ZERO_GRADIENT and ZERO_BACKFLOW " + "for 'alpha_c'.".format(', '.join(sorted(overlap)))) + + def time_derivative_terms(self, gfu_lst: List[List[GridFunction]], scheme: str, + step: int = 1): + """Add the virtual-mass cross derivatives to the standard time terms.""" + a, L = super().time_derivative_terms(gfu_lst, scheme, step) + if not self.VM_switch: + return a, L + + U, V = self.get_trial_and_test_functions() + comp = self.model_components + uc, ud = U[comp['u_c']], U[comp['u_d']] + vc, vd = V[comp['u_c']], V[comp['u_d']] + + Ac = self.UIter.components[comp['alpha_c']] + Ad = 1 - Ac + VM_d = self.rho_c[0] * self.C_VM[0] / self.rho_d[0] + VM_c = -Ad / Ac * self.C_VM[0] + + if scheme in ('explicit euler', 'implicit euler', 'crank nicolson', 'adaptive imex pred'): + current_coefficient = 1.0 + old = gfu_lst[1] + old_difference = old[comp['u_d']] - old[comp['u_c']] + elif scheme == 'CNLF': + current_coefficient = 1.0 + old = gfu_lst[2] + old_difference = old[comp['u_d']] - old[comp['u_c']] + elif scheme == 'SBDF': + current_coefficient = 11.0 / 6.0 + old_difference = ( + 3.0 * (gfu_lst[1][comp['u_d']] - gfu_lst[1][comp['u_c']]) + - 1.5 * (gfu_lst[2][comp['u_d']] - gfu_lst[2][comp['u_c']]) + + (1.0 / 3.0) * (gfu_lst[3][comp['u_d']] - gfu_lst[3][comp['u_c']]) + ) + elif scheme in ('RK 222', 'RK 232'): + current_coefficient = 1.0 + old = gfu_lst[step] + old_difference = old[comp['u_d']] - old[comp['u_c']] + else: + raise ValueError('Scheme "{}" is not implemented'.format(scheme)) + + current_difference = ud - uc + a[0] += current_coefficient * ( + VM_d * current_difference * vd + VM_c * current_difference * vc) + L[0] += VM_d * old_difference * vd + VM_c * old_difference * vc + return a, L + + def _get_drag_coeff(self, wd, wc, Ad, ts: int) -> ngs.CoefficientFunction: + rho_c = self.rho_c[ts]; rho_d = self.rho_d[ts] + nu_c = self.nu_c[ts]; dp = self.dp[ts] + sigma_c = self.sigma_c[ts] + g_mag = ngs.Norm(self.gravity) + Eo = g_mag * (rho_c - rho_d) * dp**2 / sigma_c + Re = ngs.Norm(wd - wc) * dp / nu_c + if self.drag_model == 'Tomiyama': + Cd = ngs.IfPos(Re, + Max(Min(24/Re*(1 + 0.15*Re**0.687), 72/Re), + 8/3*Eo/(Eo + 4)), + ngs.CoefficientFunction(0.0)) + else: + Cd = ngs.IfPos(Re - 1000, + ngs.CoefficientFunction(0.44), + 24/Re*(1 + 0.15*Re**0.687)) + return Cd + + def _get_lift_coeff(self, wd, wc, ts: int) -> ngs.CoefficientFunction: + rho_c = self.rho_c[ts]; rho_d = self.rho_d[ts] + nu_c = self.nu_c[ts]; dp = self.dp[ts] + sigma_c = self.sigma_c[ts] + g_mag = ngs.Norm(self.gravity) + Eo = g_mag * (rho_c - rho_d) * dp**2 / sigma_c + Re = ngs.Norm(wd - wc) * dp / nu_c + Sr = dp**2 / (Re * nu_c + 1e-30) * ngs.Norm(ngs.grad(wc)) + if self.lift_model == 'LegendreMagnaudet': + # Guarded denominator: at zero slip Re and Sr both vanish, giving 0/0. + ClLow = (6*2.255)**2 * Sr**2 / (ngs.pi**4 * Re * (Sr + 0.2*Re)**3 + 1e-10) + ClHigh = (0.5 * (Re + 16) / (Re + 29))**2 + Cl = (ClLow + ClHigh)**0.5 + else: # Tomiyama + fEo = 0.00105*Eo**3 - 0.0159*Eo**2 - 0.0204*Eo + 0.474 + Cl = (ngs.IfPos(Eo - 4, 0, 1) * Min(0.288*tanh(0.121*Re), fEo) + + ngs.IfPos(Eo - 4, 1, 0) * ngs.IfPos(10.7 - Eo, 1, 0) * fEo + + ngs.IfPos(Eo - 10.7, 1, 0) * (-0.288)) + + # Optional near-wall lift taper. + if self.lift_wall_deactivation: + dp_factor = 3 + xw = self._wall_dist + s = 2*xw/(dp_factor*dp) - 1 + activation = (ngs.IfPos(xw - dp_factor*dp, 1, 0) + + ngs.IfPos(dp_factor*dp - xw, 1, 0) + * ngs.IfPos(xw - 0.5*dp_factor*dp, 1, 0) + * (3*s**2 - 2*s**3)) + Cl = Cl * activation + return Cl + + # ------------------------------------------------------------------ + # Internal helpers — UDS numerical fluxes + # ------------------------------------------------------------------ + + def _NF_UDS_mass(self, u, w, facet: str = 'Interior', bl: bool = True): + """UDS numerical flux for scalar advection (alpha_c).""" + n = self.n + if facet == 'Interior': + return avg(u) * (w * n) + 0.5 * ngs.Norm(w * n) * jump(u) + if facet == 'Dirichlet': + sign = 1.0 if bl else -1.0 + return 0.5 * u * (w * n) + sign * 0.5 * u * ngs.Norm(w * n) + if facet == 'Neumann': + return u * Max(w * n, ngs.CoefficientFunction(0.0)) + #return u * w * n + + def _NF_UDS_mom(self, u, w, facet: str = 'Interior', bl: bool = True): + """UDS numerical flux for momentum (vector) advection.""" + n = self.n + if facet == 'Interior': + return (ngs.OuterProduct(avg(u), w) + + 0.5 * ngs.Norm(w * n) * ngs.OuterProduct(jump(u), n)) + if facet == 'Dirichlet': + sign = 1.0 if bl else -1.0 + return (0.5 * ngs.OuterProduct(u, w) + + sign * 0.5 * ngs.Norm(w * n) * ngs.OuterProduct(u, n)) + if facet == 'Neumann': + return ngs.OuterProduct(u, n) * Max(w * n, ngs.CoefficientFunction(0.0)) + + # ------------------------------------------------------------------ + # Phase 2 — bilinear form (spatial operator) + # ------------------------------------------------------------------ + + 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: + ts = time_step + n = self.n + IM = self.IM + pen = self.penalty_interior + comp = self.model_components + + # Trial / test functions + uc, ud = U[comp['u_c']], U[comp['u_d']] + p = U[comp['p']] + alpha_c = U[comp['alpha_c']] + vc, vd = V[comp['u_c']], V[comp['u_d']] + q, r = V[comp['p']], V[comp['alpha_c']] + + # Picard iterate fields + wc = self.UIter.components[comp['u_c']] + wd = self.UIter.components[comp['u_d']] + Ac = self.UIter.components[comp['alpha_c']] + Ad = 1 - Ac + + # Physical parameters + rho_c = self.rho_c[ts]; rho_d = self.rho_d[ts] + nu_c = self.nu_c[ts]; nu_d = self.nu_d[ts] + dp = self.dp[ts] + + # BC regex strings + ac_d_reg = self.dirichlet_names.get('alpha_c', '') + ac_zg_reg = self._bc_regex('zero_gradient', 'alpha_c') + ac_zb_reg = self._bc_regex('zero_backflow', 'alpha_c') + ac_n_reg = '|'.join(reg for reg in (ac_zg_reg, ac_zb_reg) if reg) + uc_d_reg = self.dirichlet_names.get('u_c', '') + uc_n_reg = self._bc_regex('zero_stress', 'u_c') + ud_d_reg = self.dirichlet_names.get('u_d', '') + ud_n_reg = self._bc_regex('zero_stress', 'u_d') + + a = ngs.CoefficientFunction(0.0) * ngs.dx + + # ============================================================ + # 1. Alpha_c (dispersed-phase mass conservation) — a terms only + # (a_dt = alpha_c*r*dx is handled by time_derivative_terms) + # ============================================================ + # div(u_d)*r is treated EXPLICITLY (lagged on the Picard wind w_d) and + # lives in the linear form as +dt*div(w_d)*r -- see construct_linear. + a += (dt * -alpha_c * (wd * ngs.grad(r))) * ngs.dx + a += (dt * jump(r) * self._NF_UDS_mass(alpha_c, wd)) * ngs.dx(skeleton=True) + if ac_d_reg: + a += (dt * r * self._NF_UDS_mass(alpha_c, wd, 'Dirichlet', True)) * self._ds(ac_d_reg) + if ac_zg_reg: + a += (dt * r * alpha_c * (wd * n)) * self._ds(ac_zg_reg) + if ac_zb_reg: + # Outflow only; the inflow half is prescribed in construct_linear. + a += (dt * r * alpha_c * Max(wd * n, ngs.CoefficientFunction(0.0))) * self._ds(ac_zb_reg) + + # Artificial diffusion (SIPG) on alpha_c -- stabilises sharp fronts. + if self.diffusion_switch: + art = self.D_art[ts] + a += (dt * art * ngs.grad(alpha_c) * ngs.grad(r)) * ngs.dx + a += (-dt * art * (n * grad_avg(r)) * jump(alpha_c)) * ngs.dx(skeleton=True) + a += (dt * art * (pen * jump(alpha_c)) * jump(r)) * ngs.dx(skeleton=True) + a += (-dt * art * (n * grad_avg(alpha_c)) * jump(r)) * ngs.dx(skeleton=True) + if ac_d_reg: + a += (-dt * art * alpha_c * (ngs.grad(r) * n)) * self._ds(ac_d_reg) + a += (dt * art * (pen * alpha_c - ngs.grad(alpha_c) * n) * r) * self._ds(ac_d_reg) + if ac_n_reg: + a += (-dt * art * (n * ngs.grad(alpha_c)) * r) * self._ds(ac_n_reg) + + # ============================================================ + # 2. Pressure / mixture mass conservation — NOT dt-scaled + # ============================================================ + a += (Ac * ngs.div(uc) * q) * ngs.dx + a += (ngs.grad(Ac) * uc * q) * ngs.dx + a += (ngs.div(ud) * q) * ngs.dx + a += (-Ac * ngs.div(ud) * q) * ngs.dx + a += (-ngs.grad(Ac) * ud * q) * ngs.dx + if self.mean_zero_pressure: + pressure_mean, pressure_mean_test = U[-1], V[-1] + a += (pressure_mean * q + p * pressure_mean_test) * ngs.dx + + # ============================================================ + # 3 & 4. Phase momentum conservation (c and d) + # ============================================================ + for (u_tr, w_pi, v_ts, A_pi, rho, nu_lam, d_reg, n_reg, phase) in [ + (uc, wc, vc, Ac, rho_c, nu_c, uc_d_reg, uc_n_reg, 'c'), + (ud, wd, vd, Ad, rho_d, nu_d, ud_d_reg, ud_n_reg, 'd'), + ]: + tau = ngs.grad(u_tr) + ngs.grad(u_tr).trans - 2.0/3.0 * ngs.div(u_tr) * IM + + # Advection — bulk (Picard-linearised convection) + a += (dt * -ngs.div(w_pi) * u_tr * v_ts) * ngs.dx + a += (dt * -ngs.InnerProduct(ngs.grad(v_ts), ngs.OuterProduct(u_tr, w_pi))) * ngs.dx + + # Pressure (shared mixture pressure p) + a += (dt * -p / rho * ngs.div(v_ts)) * ngs.dx + + # Advection — interior facets + a += (dt * ngs.InnerProduct(ngs.OuterProduct(jump(v_ts), n), + self._NF_UDS_mom(u_tr, w_pi))) * ngs.dx(skeleton=True) + # Advection — Dirichlet boundary (bilinear) + if d_reg: + a += (dt * ngs.InnerProduct(ngs.OuterProduct(v_ts, n), + self._NF_UDS_mom(u_tr, w_pi, 'Dirichlet', True))) \ + * self._ds(d_reg) + # Advection — Neumann boundary (bilinear) + if n_reg: + a += (dt * ngs.InnerProduct(ngs.OuterProduct(v_ts, n), + self._NF_UDS_mom(u_tr, w_pi, 'Neumann', True))) \ + * self._ds(n_reg) + + # Viscous: B-TFM (phase c only) or C-TFM (both phases) + do_viscous = (self.canonical_form == 'B-TFM' and phase == 'c') or \ + (self.canonical_form == 'C-TFM') + if do_viscous: + nu = ngs.CoefficientFunction(nu_lam) + + # grad of the phase fraction. A_d = 1 - A_c is a CF expression and + # ngs.grad() only accepts a GridFunction/proxy, so use grad(A_d) = + # -grad(A_c) with the actual A_c iterate (also tracks Picard updates). + grad_A = ngs.grad(Ac) if phase == 'c' else -ngs.grad(Ac) + + if self.canonical_form == 'B-TFM': + factor = nu / A_pi + factor_avg = avg(nu / A_pi) + a += (dt * -nu * ngs.InnerProduct( + tau, ngs.OuterProduct(v_ts, grad_A / (A_pi**2)))) * ngs.dx + else: # C-TFM + factor = ngs.CoefficientFunction(nu) + factor_avg = avg(ngs.CoefficientFunction(nu)) + if phase == 'c': + a += (dt * -nu * ngs.InnerProduct( + tau, ngs.OuterProduct(v_ts, grad_A / A_pi))) * ngs.dx + else: + a += (dt * -nu * ngs.InnerProduct( + tau, ngs.OuterProduct(v_ts, + grad_A / (A_pi + 1e-5)))) * ngs.dx + + # Bulk viscous + a += (dt * factor * ngs.InnerProduct(ngs.grad(v_ts), tau)) * ngs.dx + + # Interior IP-DG (SIPG for viscous stress) + stress_avg = (weighted_grad_avg(u_tr, factor) + weighted_grad_avg(u_tr, factor).trans + - 2.0/3.0 * weighted_div_avg(u_tr, factor) * IM) + a += (dt * -ngs.InnerProduct(stress_avg, + ngs.OuterProduct(jump(v_ts), n))) * ngs.dx(skeleton=True) + a += (dt * factor_avg * pen + * ngs.InnerProduct(ngs.OuterProduct(jump(u_tr), n), + ngs.OuterProduct(jump(v_ts), n))) * ngs.dx(skeleton=True) + # The reference form uses lower-case grad(v), which is Grad(v)^T + # for a vector. Transpose the outer product to express that term + # with OpenCMP's Grad-based helper. + a += (dt * -ngs.InnerProduct(ngs.OuterProduct(n, jump(u_tr)), + weighted_grad_avg(v_ts, factor))) * ngs.dx(skeleton=True) + a += (dt * -ngs.InnerProduct(ngs.OuterProduct(jump(u_tr), n), + weighted_grad_avg(v_ts, factor))) * ngs.dx(skeleton=True) + + # Dirichlet boundary — viscous bilinear (Nitsche) + if d_reg: + a += (dt * -factor * ngs.InnerProduct(tau, ngs.OuterProduct(v_ts, n))) \ + * self._ds(d_reg) + a += (dt * factor * pen + * ngs.InnerProduct(ngs.OuterProduct(u_tr, n), + ngs.OuterProduct(v_ts, n))) * self._ds(d_reg) + a += (dt * -factor * ngs.InnerProduct( + ngs.OuterProduct(u_tr, n), ngs.grad(v_ts))) * self._ds(d_reg) + a += (dt * -factor * ngs.InnerProduct( + ngs.OuterProduct(u_tr, n), ngs.grad(v_ts).trans)) * self._ds(d_reg) + + # ============================================================ + # 5. Interphase closures — bilinear contributions + # ============================================================ + + # — Drag (dt-scaled) + if self.drag_switch: + Cd = self._get_drag_coeff(wd, wc, Ad, ts) + drag_factor_d = 0.75 * Cd * rho_c / (rho_d * dp) * ngs.Norm(wd - wc) + drag_factor_c = -0.75 * Cd * Ad / (Ac * dp) * ngs.Norm(wd - wc) + a += (dt * drag_factor_d * (ud - uc) * vd) * ngs.dx + a += (dt * drag_factor_c * (ud - uc) * vc) * ngs.dx + + # — Virtual mass: a_dt (NOT dt-scaled) + a (dt-scaled) + if self.VM_switch: + C_VM = self.C_VM[ts] + VM_d = rho_c * C_VM / rho_d # phase-d coefficient (positive) + VM_c = -Ad / Ac * C_VM # phase-c coefficient (negative) + + # Spatial VM terms (dt-scaled), phase d + a += (dt * VM_d * (-ngs.InnerProduct(ngs.grad(vd), ngs.OuterProduct(ud, wd)) + + ngs.InnerProduct(ngs.grad(vd), ngs.OuterProduct(uc, wc)) + + (-ngs.div(wd)*ud*vd + ngs.div(wc)*uc*vd))) * ngs.dx + a += (dt * VM_d * ngs.InnerProduct(ngs.OuterProduct(jump(vd), n), + self._NF_UDS_mom(ud, wd))) * ngs.dx(skeleton=True) + a += (dt * -VM_d * ngs.InnerProduct(ngs.OuterProduct(jump(vd), n), + self._NF_UDS_mom(uc, wc))) * ngs.dx(skeleton=True) + if ud_d_reg: + a += (dt * VM_d * ngs.InnerProduct(ngs.OuterProduct(vd, n), + self._NF_UDS_mom(ud, wd, 'Dirichlet', True))) \ + * self._ds(ud_d_reg) + if uc_d_reg: + a += (dt * -VM_d * ngs.InnerProduct(ngs.OuterProduct(vd, n), + self._NF_UDS_mom(uc, wc, 'Dirichlet', True))) \ + * self._ds(uc_d_reg) + if ud_n_reg: + a += (dt * VM_d * ngs.InnerProduct(ngs.OuterProduct(vd, n), + self._NF_UDS_mom(ud, wd, 'Neumann', True))) \ + * self._ds(ud_n_reg) + if uc_n_reg: + a += (dt * -VM_d * ngs.InnerProduct(ngs.OuterProduct(vd, n), + self._NF_UDS_mom(uc, wc, 'Neumann', True))) \ + * self._ds(uc_n_reg) + + # Spatial VM terms (dt-scaled), phase c + a += (dt * VM_c * (-ngs.InnerProduct(ngs.grad(vc), ngs.OuterProduct(ud, wd)) + + ngs.InnerProduct(ngs.grad(vc), ngs.OuterProduct(uc, wc)) + + (-ngs.div(wd)*ud*vc + ngs.div(wc)*uc*vc))) * ngs.dx + grad_VM_c = C_VM * ngs.grad(Ac) / (Ac**2) + vm_flux_difference = (ngs.OuterProduct(ud, wd) + - ngs.OuterProduct(uc, wc)) + a += (dt * -(vm_flux_difference.trans * grad_VM_c) * vc) * ngs.dx + a += (dt * VM_c * ngs.InnerProduct(ngs.OuterProduct(jump(vc), n), + self._NF_UDS_mom(ud, wd))) * ngs.dx(skeleton=True) + a += (dt * -VM_c * ngs.InnerProduct(ngs.OuterProduct(jump(vc), n), + self._NF_UDS_mom(uc, wc))) * ngs.dx(skeleton=True) + if ud_d_reg: + a += (dt * VM_c * ngs.InnerProduct(ngs.OuterProduct(vc, n), + self._NF_UDS_mom(ud, wd, 'Dirichlet', True))) \ + * self._ds(ud_d_reg) + if uc_d_reg: + a += (dt * -VM_c * ngs.InnerProduct(ngs.OuterProduct(vc, n), + self._NF_UDS_mom(uc, wc, 'Dirichlet', True))) \ + * self._ds(uc_d_reg) + if ud_n_reg: + a += (dt * VM_c * ngs.InnerProduct(ngs.OuterProduct(vc, n), + self._NF_UDS_mom(ud, wd, 'Neumann', True))) \ + * self._ds(ud_n_reg) + if uc_n_reg: + a += (dt * -VM_c * ngs.InnerProduct(ngs.OuterProduct(vc, n), + self._NF_UDS_mom(uc, wc, 'Neumann', True))) \ + * self._ds(uc_n_reg) + + # — Lift (dt-scaled) + if self.Lift_switch: + Cl = self._get_lift_coeff(wd, wc, ts) + if self.mesh.dim == 2: + curl_wc = ngs.grad(wc)[1] - ngs.grad(wc)[2] + lift_mom_d = Cl * rho_c/rho_d * curl_wc * ( + (ud[1] - uc[1]) * vd[0] - (ud[0] - uc[0]) * vd[1]) + lift_mom_c = Cl * (-Ad)/(Ac + 1e-30) * curl_wc * ( + (ud[1] - uc[1]) * vc[0] - (ud[0] - uc[0]) * vc[1]) + else: + lift_mom_d = Cl * rho_c/rho_d * (ngs.Cross(ud - uc, curl_3d(wc)) * vd) + lift_mom_c = Cl * (-Ad)/(Ac + 1e-30) * (ngs.Cross(ud - uc, curl_3d(wc)) * vc) + a += (dt * lift_mom_d) * ngs.dx + a += (dt * lift_mom_c) * ngs.dx + + return [a] + + def construct_bilinear_time_coefficient(self, + U: List[ProxyFunction], + V: List[ProxyFunction], + dt: Parameter, + time_step: int) -> List: + return [ngs.CoefficientFunction(0.0) * ngs.dx] + + # ------------------------------------------------------------------ + # Phase 2 — linear form (source and BC terms) + # ------------------------------------------------------------------ + + def construct_linear(self, + V: List[ProxyFunction], + gfu_0: Optional[List[GridFunction]], + dt: Parameter, + time_step: int) -> List: + ts = time_step + n = self.n + pen = self.penalty_interior + comp = self.model_components + + vc, vd = V[comp['u_c']], V[comp['u_d']] + q, r = V[comp['p']], V[comp['alpha_c']] + + # Picard iterate (for nonlinear coefficient terms) + wc = self.UIter.components[comp['u_c']] + wd = self.UIter.components[comp['u_d']] + Ac = self.UIter.components[comp['alpha_c']] + Ad = 1 - Ac + + # Physical parameters + rho_c = self.rho_c[ts]; rho_d = self.rho_d[ts] + nu_c = self.nu_c[ts]; nu_d = self.nu_d[ts] + dp = self.dp[ts] + + # BC regex + ac_d_reg = self.dirichlet_names.get('alpha_c', '') + ac_zg_reg = self._bc_regex('zero_gradient', 'alpha_c') + ac_zb_reg = self._bc_regex('zero_backflow', 'alpha_c') + ac_n_reg = '|'.join(reg for reg in (ac_zg_reg, ac_zb_reg) if reg) + uc_d_reg = self.dirichlet_names.get('u_c', '') + uc_n_reg = self._bc_regex('zero_stress', 'u_c') + ud_d_reg = self.dirichlet_names.get('u_d', '') + ud_n_reg = self._bc_regex('zero_stress', 'u_d') + + L = ngs.CoefficientFunction(0.0) * ngs.dx + + # Optional manufactured/general volume sources. Their signs follow the + # four strong equations represented by (u_c, u_d, p, alpha_c). + if 'u_c' in self.f: + L += (dt * self.f['u_c'][ts] * vc) * ngs.dx + if 'u_d' in self.f: + L += (dt * self.f['u_d'][ts] * vd) * ngs.dx + if 'p' in self.f: + L += (self.f['p'][ts] * q) * ngs.dx + if 'alpha_c' in self.f: + L += (dt * self.f['alpha_c'][ts] * r) * ngs.dx + + # ============================================================ + # 1. Alpha_c — linear contributions + # ============================================================ + # Explicit (lagged) dispersed-phase compressibility term: the bilinear + # form drops -div(u_d)*r; it is carried here on the Picard wind w_d. + L += (dt * ngs.div(wd) * r) * ngs.dx + for marker, val_list in self.BC.get('dirichlet', {}).get('alpha_c', {}).items(): + val = val_list[ts] + L += (dt * -r * self._NF_UDS_mass(val, wd, 'Dirichlet', False)) * self._ds(marker) + if ac_zb_reg: + # Backflow enters as pure continuous phase (alpha_c = 1). + L += (dt * -r * Min(wd * n, ngs.CoefficientFunction(0.0))) * self._ds(ac_zb_reg) + # Artificial diffusion — Dirichlet boundary (linear part). + if self.diffusion_switch: + art = self.D_art[ts] + for marker, val_list in self.BC.get('dirichlet', {}).get('alpha_c', {}).items(): + val = val_list[ts] + L += (dt * art * pen * val * r) * self._ds(marker) + L += (-dt * art * val * (ngs.grad(r) * n)) * self._ds(marker) + + # ============================================================ + # 2 & 3. Phase momentum — linear terms (gravity + BCs) + # ============================================================ + for (v_ts, w_pi, A_pi, rho, nu_lam, d_var, phase) in [ + (vc, wc, Ac, rho_c, nu_c, 'u_c', 'c'), + (vd, wd, Ad, rho_d, nu_d, 'u_d', 'd'), + ]: + # Body force + L += (dt * self.gravity * v_ts) * ngs.dx + + # Viscous parameters (needed for Nitsche BC terms) + do_viscous = (self.canonical_form == 'B-TFM' and phase == 'c') or \ + (self.canonical_form == 'C-TFM') + if do_viscous: + nu_eff = ngs.CoefficientFunction(nu_lam) + factor = nu_eff / A_pi if self.canonical_form == 'B-TFM' else ngs.CoefficientFunction(nu_eff) + + # Convective Dirichlet BC (linear part: outflow from prescribed BC) + for marker, u_bc_list in self.BC.get('dirichlet', {}).get(d_var, {}).items(): + u_bc = u_bc_list[ts] + L += (dt * -ngs.InnerProduct(ngs.OuterProduct(v_ts, n), + self._NF_UDS_mom(u_bc, w_pi, 'Dirichlet', False))) \ + * self._ds(marker) + # Nitsche viscous Dirichlet (linear part) + if do_viscous: + L += (dt * factor * pen + * ngs.InnerProduct(ngs.OuterProduct(u_bc, n), + ngs.OuterProduct(v_ts, n))) * self._ds(marker) + L += (dt * -factor * ngs.InnerProduct( + ngs.OuterProduct(u_bc, n), ngs.grad(v_ts))) * self._ds(marker) + L += (dt * -factor * ngs.InnerProduct( + ngs.OuterProduct(u_bc, n), ngs.grad(v_ts).trans)) * self._ds(marker) + + # ZERO_STRESS contributes no data by definition. + + # ============================================================ + # 4. Virtual mass — l_dt (NOT dt-scaled, uses UOld) + l (dt-scaled) + # ============================================================ + if self.VM_switch: + C_VM = self.C_VM[ts] + VM_d = rho_c * C_VM / rho_d + VM_c = -Ad / Ac * C_VM + + # l (dt-scaled): convective Dirichlet BC contributions + for marker, u_bc_list in self.BC.get('dirichlet', {}).get('u_d', {}).items(): + u_bc = u_bc_list[ts] + L += (dt * -VM_d * ngs.InnerProduct(ngs.OuterProduct(vd, n), + self._NF_UDS_mom(u_bc, wd, 'Dirichlet', False))) \ + * self._ds(marker) + L += (dt * -VM_c * ngs.InnerProduct(ngs.OuterProduct(vc, n), + self._NF_UDS_mom(u_bc, wd, 'Dirichlet', False))) \ + * self._ds(marker) + for marker, u_bc_list in self.BC.get('dirichlet', {}).get('u_c', {}).items(): + u_bc = u_bc_list[ts] + L += (dt * VM_d * ngs.InnerProduct(ngs.OuterProduct(vd, n), + self._NF_UDS_mom(u_bc, wc, 'Dirichlet', False))) \ + * self._ds(marker) + L += (dt * VM_c * ngs.InnerProduct(ngs.OuterProduct(vc, n), + self._NF_UDS_mom(u_bc, wc, 'Dirichlet', False))) \ + * self._ds(marker) + + # ============================================================ + # 5. Laminar dispersion (linear source term) + # ============================================================ + if self.Disp_switch and self.drag_switch: + Cd = self._get_drag_coeff(wd, wc, Ad, ts) + H = 1 - 1.166*(1 - Ac) + 0.5*(1 - Ac)**2 + Cdis = self.Cdis[ts] + disp_d = 0.75*Cd*Cdis*rho_c/rho_d * H * ngs.Norm(wd - wc)**2 + disp_c = -0.75*Cd*Cdis*Ad/(Ac + 1e-30) * H * ngs.Norm(wd - wc)**2 + L += (dt * disp_d * ngs.grad(Ac) * vd) * ngs.dx + L += (dt * disp_c * ngs.grad(Ac) * vc) * ngs.dx + + # ============================================================ + # 6. Interior injection source (on the 'injection' material region) + # ============================================================ + if self.injection_switch and self.inj_mass_flowrate: + mdot = self.inj_mass_flowrate + dx_inj = ngs.dx(definedon=self.mesh.Materials(self.injection_region)) + L += (dt * -mdot * r) * dx_inj # dispersed-phase mass + L += (dt * mdot * q) * dx_inj # mixture mass + if self.mesh.dim == 2: + inj_vec = ngs.CoefficientFunction((0, self.inj_velocity)) + else: + inj_vec = ngs.CoefficientFunction((0, self.inj_velocity, 0)) + L += (dt * mdot * inj_vec * vc) * dx_inj # continuous-phase momentum + + return [L] + + def construct_imex_explicit(self, + V: List[ProxyFunction], + gfu_0: Optional[List[GridFunction]], + dt: Parameter, + time_step: int) -> List: + return [ngs.CoefficientFunction(0.0) * ngs.dx] + + # ------------------------------------------------------------------ + # Phase 3 — nonlinear (Picard) solver + # ------------------------------------------------------------------ + + def solve_single_step(self, + a_lst: List[BilinearForm], + L_lst: List[LinearForm], + precond_lst: List[Preconditioner], + gfu: GridFunction, + time_step: int = 0) -> bool: + # Cold-start guard: the adaptive solvers hand us fresh GridFunctions + # (gfu_long/gfu_short) on the first step with only the velocity Dirichlet + # BCs applied -- alpha_c is still zero. A zero alpha_c is a degenerate + # linearization wind for this model (drag /Ac -> /0), so seed it from the + # initial condition before it is used as UIter. (Keys on alpha_c since + # velocity BCs leave it untouched.) + comp = self.model_components + if gfu.components[comp['alpha_c']].vec.Norm() == 0.0: + gfu.vec.data = self.IC.vec + + gfu_prev = ngs.GridFunction(self.fes) + + for _it in range(self.nonlinear_max_iters): + gfu_prev.vec.data = gfu.vec + + # Update Picard wind so that Assemble() picks up the new values + self.UIter.vec.data = gfu.vec + + 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) + + # Per-component relaxation + for i, rf in enumerate(self.relax_factors): + if rf < 1.0: + gfu.components[i].vec.data = ( + rf * gfu.components[i].vec + + (1.0 - rf) * gfu_prev.components[i].vec) + + # Bound-preserving Bezier/Bernstein limiter: GUARANTEES the per-element + # polynomial stays in bounds everywhere. Order is the + # scalar L2 order used for alpha_c in _construct_fes. + if self.slope_limiter: + ord_s = max(self.interp_ord - 1, 0) + ac = gfu.components[comp['alpha_c']] + self._lim.bezier_bound(ac, ac.space, ord_s, (0.0, 1.0)) + + # Match INS convergence testing: evaluate each local-error component + # independently and exclude pressure. + converged = True + for name, include in self.model_local_error_components.items(): + if not include: + continue + i = comp[name] + fes_component = self.fes.components[i] + err = norm('l2_norm', gfu_prev.components[i], gfu.components[i], + self.mesh, fes_component, average=False) + solution_norm = mean(gfu.components[i], self.mesh) + tolerance = (self.abs_nonlinear_tolerance + + self.rel_nonlinear_tolerance * solution_norm) + if err >= tolerance: + converged = False + break + + if converged: + logging.info(f'TFM Picard converged in {_it + 1} iteration(s).') + return True + + logging.warning('TFM Picard did NOT converge within ' + f'{self.nonlinear_max_iters} iterations.') + return False + + def update_linearization(self, gfu: GridFunction) -> None: + self.UIter.vec.data = gfu.vec + + def linearized_solve(self, a_assembled: BilinearForm, L_assembled: LinearForm, + precond: Preconditioner, gfu: GridFunction): + """Perform one stationary Picard iteration and report its change.""" + previous = ngs.GridFunction(self.fes) + previous.vec.data = gfu.vec + self.linear_solve(a_assembled, L_assembled, precond, gfu) + + # The stationary solver performs one Picard solve per outer iteration, + # so apply the same component-wise under-relaxation used by + # solve_single_step for transient TFM solves. + for i, relaxation in enumerate(self.relax_factors): + if relaxation < 1.0: + gfu.components[i].vec.data = ( + relaxation * gfu.components[i].vec + + (1.0 - relaxation) * previous.components[i].vec) + + error = 0.0 + solution_norm = 0.0 + for name, include in self.model_local_error_components.items(): + if not include: + continue + component = self.model_components[name] + space = self.fes.components[component] + error = max(error, norm('l2_norm', previous.components[component], + gfu.components[component], self.mesh, space, + average=False)) + solution_norm = max(solution_norm, mean(gfu.components[component], self.mesh)) + return error, solution_norm diff --git a/pytests/full_system/tfm/bc_dir/bc_config b/pytests/full_system/tfm/bc_dir/bc_config new file mode 100644 index 0000000..c2b1104 --- /dev/null +++ b/pytests/full_system/tfm/bc_dir/bc_config @@ -0,0 +1,15 @@ +[DIRICHLET] +u_c = bottom -> [1-0.01*cos(pi*x)*sin(pi*y), 0.01*sin(pi*x)*cos(pi*y)] + right -> [1-0.01*cos(pi*x)*sin(pi*y), 0.01*sin(pi*x)*cos(pi*y)] + top -> [1-0.01*cos(pi*x)*sin(pi*y), 0.01*sin(pi*x)*cos(pi*y)] + left -> [1-0.01*cos(pi*x)*sin(pi*y), 0.01*sin(pi*x)*cos(pi*y)] +u_d = bottom -> [1+0.01*sin(pi*x)*sin(pi*y), 0.01*cos(pi*x)*cos(pi*y)] + right -> [1+0.01*sin(pi*x)*sin(pi*y), 0.01*cos(pi*x)*cos(pi*y)] + top -> [1+0.01*sin(pi*x)*sin(pi*y), 0.01*cos(pi*x)*cos(pi*y)] + left -> [1+0.01*sin(pi*x)*sin(pi*y), 0.01*cos(pi*x)*cos(pi*y)] +alpha_c = bottom -> 0.9+0.002*sin(pi*x+pi/4)*cos(pi*y) + right -> 0.9+0.002*sin(pi*x+pi/4)*cos(pi*y) + top -> 0.9+0.002*sin(pi*x+pi/4)*cos(pi*y) + left -> 0.9+0.002*sin(pi*x+pi/4)*cos(pi*y) + +[SLIP] diff --git a/pytests/full_system/tfm/config b/pytests/full_system/tfm/config new file mode 100644 index 0000000..25fedd1 --- /dev/null +++ b/pytests/full_system/tfm/config @@ -0,0 +1,47 @@ +[MESH] +filename = pytests/mesh_files/unit_square_coarse.vol + +[DIM] +diffuse_interface_method = False + +[FINITE ELEMENT SPACE] +elements = u_c -> HDiv + u_d -> HDiv + p -> L2 + alpha_c -> L2 +interpolant_order = 2 + +[DG] +DG = True +interior_penalty_coefficient = 20.0 + +[SOLVER] +linear_solver = direct +preconditioner = default +nonlinear_solver = NoMixing +nonlinear_tolerance = relative -> 1e-5 + absolute -> 1e-8 +nonlinear_max_iterations = 100 +relaxation_factors = 0.2, 0.2, 0.2, 0.2 + +[TRANSIENT] +transient = False + +[VISUALIZATION] +save_to_file = False + +[ERROR ANALYSIS] +check_error = False + +[OTHER] +num_threads = 1 +model = TwoFluidModel +run_dir = pytests/full_system/tfm +slope_limiter = False +diffusion_switch = False +mean_zero_pressure = True + +[TFM] +canonical_form = C-TFM +IME = +lift_wall_deactivation = False diff --git a/pytests/full_system/tfm/ic_dir/ic_config b/pytests/full_system/tfm/ic_dir/ic_config new file mode 100644 index 0000000..480f08f --- /dev/null +++ b/pytests/full_system/tfm/ic_dir/ic_config @@ -0,0 +1,5 @@ +[TwoFluidModel] +u_c = all -> [1.0, 0.0] +u_d = all -> [1.0, 0.0] +p = all -> 0.0 +alpha_c = all -> 0.9 diff --git a/pytests/full_system/tfm/import_functions.py b/pytests/full_system/tfm/import_functions.py new file mode 100644 index 0000000..5581458 --- /dev/null +++ b/pytests/full_system/tfm/import_functions.py @@ -0,0 +1,159 @@ +import ngsolve as ngs +from opencmp.helpers.math import Max, Min, tanh +from opencmp.models.tfm import TwoFluidModel + + +ACTIVE_IME = frozenset() + + +def set_active_ime(ime): + global ACTIVE_IME + ACTIVE_IME = frozenset(ime) + + +def ime_config(ime): + """Build a [TFM] IME value from mechanism names, using each one's first closure model.""" + return '\n'.join('{} -> {}'.format(mechanism, TwoFluidModel.IME_MODELS[mechanism][0]) + for mechanism in ime) + + +def exact_solution(t=None): + k = ngs.pi + uc = ngs.CoefficientFunction(( + 1.0 - 0.01 * ngs.cos(k * ngs.x) * ngs.sin(k * ngs.y), + 0.01 * ngs.sin(k * ngs.x) * ngs.cos(k * ngs.y), + )) + ud = ngs.CoefficientFunction(( + 1.0 + 0.01 * ngs.sin(k * ngs.x) * ngs.sin(k * ngs.y), + 0.01 * ngs.cos(k * ngs.x) * ngs.cos(k * ngs.y), + )) + p = -0.0025 * (ngs.cos(2 * k * ngs.x) + ngs.cos(2 * k * ngs.y)) + alpha_c = 0.9 + 0.002 * ngs.sin(k * ngs.x + ngs.pi / 4) * ngs.cos(k * ngs.y) + return uc, ud, p, alpha_c + + +def _value(t_param, time_step, component): + return exact_solution(t_param[time_step])[component] + + +def exact_uc(t_param, model_variables, mesh, time_step): + return _value(t_param, time_step, 0) + + +def exact_ud(t_param, model_variables, mesh, time_step): + return _value(t_param, time_step, 1) + + +def exact_p(t_param, model_variables, mesh, time_step): + return _value(t_param, time_step, 2) + + +def exact_alpha_c(t_param, model_variables, mesh, time_step): + return _value(t_param, time_step, 3) + + +def _div_vector(v): + return v[0].Diff(ngs.x) + v[1].Diff(ngs.y) + + +def _grad_scalar(s): + return ngs.CoefficientFunction((s.Diff(ngs.x), s.Diff(ngs.y))) + + +def _grad_vector(v): + return ngs.CoefficientFunction(( + (v[0].Diff(ngs.x), v[1].Diff(ngs.x)), + (v[0].Diff(ngs.y), v[1].Diff(ngs.y)), + ), dims=(2, 2)) + + +def _div_tensor(T): + return ngs.CoefficientFunction(( + T[0].Diff(ngs.x) + T[2].Diff(ngs.y), + T[1].Diff(ngs.x) + T[3].Diff(ngs.y), + )) + + +def source_alpha_c(t_param, model_variables, mesh, time_step): + _, ud, _, alpha_c = exact_solution() + return _div_vector(alpha_c * ud) - _div_vector(ud) + + +def source_mixture_mass(t_param, model_variables, mesh, time_step): + uc, ud, _, alpha_c = exact_solution() + return _div_vector(alpha_c * uc) + _div_vector((1 - alpha_c) * ud) + + +def _momentum_source(phase): + uc, ud, p, alpha_c = exact_solution() + if phase == 'c': + u, rho, nu, alpha = uc, 1.0, 0.01, alpha_c + else: + u, rho, nu, alpha = ud, 1.0, 0.01, 1 - alpha_c + convection = _div_tensor(ngs.OuterProduct(u, u)) - u * _div_vector(u) + source = convection + _grad_scalar(p) / rho + grad_u = _grad_vector(u) + identity = ngs.CoefficientFunction(((1, 0), (0, 1)), dims=(2, 2)) + stress = grad_u + grad_u.trans - (2.0 / 3.0) * _div_vector(u) * identity + alpha_denominator = alpha if phase == 'c' else alpha + 1e-5 + source -= nu * (_div_tensor(stress) + + stress.trans * _grad_scalar(alpha) / alpha_denominator) + gravity = ngs.CoefficientFunction((0.0, -9.81)) + source -= gravity + + relative_velocity = ud - uc + relative_speed = ngs.Norm(relative_velocity) + rho_c = rho_d = 1.0 + dp = 1.0 + c_vm = 0.2 + cdis = 0.1 + ad = 1 - alpha_c + + re = relative_speed * dp / 0.01 + cd = ngs.IfPos( + re, + Max(Min(24 / re * (1 + 0.15 * re**0.687), 72 / re), 0.0), + ngs.CoefficientFunction(0.0), + ) + + if 'drag' in ACTIVE_IME: + if phase == 'd': + source += 0.75 * cd * rho_c / (rho_d * dp) * relative_speed * relative_velocity + else: + source += -0.75 * cd * ad / (alpha_c * dp) * relative_speed * relative_velocity + + acceleration_difference = ( + _div_tensor(ngs.OuterProduct(ud, ud)) - ud * _div_vector(ud) + - _div_tensor(ngs.OuterProduct(uc, uc)) + uc * _div_vector(uc) + ) + if 'virtual_mass' in ACTIVE_IME: + coefficient = (rho_c * c_vm / rho_d if phase == 'd' + else -ad * c_vm / alpha_c) + source += coefficient * acceleration_difference + + if 'lift' in ACTIVE_IME: + # Tomiyama coefficient with Eo=0 for this equal-density MMS case. + cl = Min(0.288 * tanh(0.121 * re), 0.474) + curl_uc = _grad_vector(uc)[1] - _grad_vector(uc)[2] + lift_vector = ngs.CoefficientFunction((relative_velocity[1], + -relative_velocity[0])) + coefficient = (rho_c / rho_d if phase == 'd' else -ad / alpha_c) + source += cl * coefficient * curl_uc * lift_vector + + if 'laminar_dispersion' in ACTIVE_IME: + hindered = 1 - 1.166 * ad + 0.5 * ad**2 + coefficient = (0.75 * cd * cdis * rho_c / rho_d * hindered * relative_speed**2 + if phase == 'd' + else -0.75 * cd * cdis * ad / alpha_c * hindered * relative_speed**2) + # Dispersion is assembled on the linear-form side in TwoFluidModel. + source -= coefficient * _grad_scalar(alpha_c) + + return source + + +def source_uc(t_param, model_variables, mesh, time_step): + return _momentum_source('c') + + +def source_ud(t_param, model_variables, mesh, time_step): + return _momentum_source('d') diff --git a/pytests/full_system/tfm/model_dir/model_config b/pytests/full_system/tfm/model_dir/model_config new file mode 100644 index 0000000..0139901 --- /dev/null +++ b/pytests/full_system/tfm/model_dir/model_config @@ -0,0 +1,16 @@ +[PARAMETERS] +rho_c = all -> 1.0 +rho_d = all -> 1.0 +nu_c = all -> 0.01 +nu_d = all -> 0.01 +sigma_c = all -> 0.072 +dp = all -> 1.0 +c_vm = all -> 0.2 +cdis = all -> 0.1 +d_artificial = all -> 0.0 + +[FUNCTIONS] +source = u_c -> IMPORT(source_uc) + u_d -> IMPORT(source_ud) + p -> IMPORT(source_mixture_mass) + alpha_c -> IMPORT(source_alpha_c) diff --git a/pytests/full_system/tfm/ref_sol_dir/ref_sol_config b/pytests/full_system/tfm/ref_sol_dir/ref_sol_config new file mode 100644 index 0000000..73df1bf --- /dev/null +++ b/pytests/full_system/tfm/ref_sol_dir/ref_sol_config @@ -0,0 +1,8 @@ +[REFERENCE SOLUTIONS] +u_c = [-cos(2*pi*x)*sin(2*pi*y), sin(2*pi*x)*cos(2*pi*y)] +u_d = [sin(2*pi*x)*sin(2*pi*y), cos(2*pi*x)*cos(2*pi*y)] +p = -0.25*(cos(4*pi*x)+cos(4*pi*y)) +alpha_c = 0.9+0.02*sin(2*pi*x+pi/4)*cos(2*pi*y) + +[METRICS] +l2_norm = u_c, u_d, p, alpha_c diff --git a/pytests/full_system/tfm/test_tfm_conservation.py b/pytests/full_system/tfm/test_tfm_conservation.py new file mode 100644 index 0000000..ea96ae4 --- /dev/null +++ b/pytests/full_system/tfm/test_tfm_conservation.py @@ -0,0 +1,188 @@ +import math +from pathlib import Path + +import ngsolve as ngs +import pytest + +from opencmp.helpers.math import Max, Min, tanh +from opencmp.config_functions import ConfigParser +from opencmp.models import get_model_class +from opencmp.solvers import get_solver_class +from import_functions import (_div_tensor, _div_vector, _grad_scalar, + exact_solution, ime_config, set_active_ime) + + +def _integrate_vector(field, mesh): + return tuple(ngs.Integrate(field[i], mesh) for i in range(mesh.dim)) + + +def _physical_ime_forces(mechanism): + uc, ud, _, ac = exact_solution() + ad = 1 - ac + rho_c = rho_d = dp = 1.0 + nu_c = 0.01 + c_vm = 0.2 + cdis = 0.1 + relative = ud - uc + speed = ngs.Norm(relative) + re = speed * dp / nu_c + cd = ngs.IfPos( + re, + Max(Min(24 / re * (1 + 0.15 * re**0.687), 72 / re), 0.0), + ngs.CoefficientFunction(0.0), + ) + + per_mass_d = ngs.CoefficientFunction((0.0, 0.0)) + per_mass_c = ngs.CoefficientFunction((0.0, 0.0)) + + if mechanism == 'drag': + per_mass_d = 0.75 * cd * rho_c / (rho_d * dp) * speed * relative + per_mass_c = -0.75 * cd * ad / (ac * dp) * speed * relative + elif mechanism == 'laminar_dispersion': + hindered = 1 - 1.166 * ad + 0.5 * ad**2 + per_mass_d = (0.75 * cd * cdis * rho_c / rho_d + * hindered * speed**2 * _grad_scalar(ac)) + per_mass_c = (-0.75 * cd * cdis * ad / ac + * hindered * speed**2 * _grad_scalar(ac)) + elif mechanism == 'virtual_mass': + acceleration_difference = ( + _div_tensor(ngs.OuterProduct(ud, ud)) - ud * _div_vector(ud) + - _div_tensor(ngs.OuterProduct(uc, uc)) + uc * _div_vector(uc) + ) + per_mass_d = rho_c * c_vm / rho_d * acceleration_difference + per_mass_c = -ad * c_vm / ac * acceleration_difference + elif mechanism == 'lift': + cl = Min(0.288 * tanh(0.121 * re), 0.474) + grad_uc = ngs.CoefficientFunction(( + (uc[0].Diff(ngs.x), uc[1].Diff(ngs.x)), + (uc[0].Diff(ngs.y), uc[1].Diff(ngs.y)), + ), dims=(2, 2)) + curl_uc = grad_uc[1] - grad_uc[2] + lift_direction = ngs.CoefficientFunction((relative[1], -relative[0])) + per_mass_d = cl * rho_c / rho_d * curl_uc * lift_direction + per_mass_c = -cl * ad / ac * curl_uc * lift_direction + else: + raise ValueError(mechanism) + + return ad * rho_d * per_mass_d, ac * rho_c * per_mass_c + + +def test_global_mass_conservation(capsys): + mesh = ngs.Mesh('pytests/mesh_files/unit_square_coarse.vol') + for _ in range(3): + mesh.Refine() + + uc, ud, _, ac = exact_solution() + ad = 1 - ac + normal = ngs.specialcf.normal(mesh.dim) + + dispersed_source = _div_vector(ad * ud) + mixture_source = _div_vector(ac * uc) + _div_vector(ad * ud) + + dispersed_volume = ngs.Integrate(dispersed_source, mesh) + dispersed_boundary = ngs.Integrate(ad * ud * normal, mesh, ngs.BND) + mixture_volume = ngs.Integrate(mixture_source, mesh) + mixture_boundary = ngs.Integrate((ac * uc + ad * ud) * normal, + mesh, ngs.BND) + + dispersed_residual = dispersed_volume - dispersed_boundary + mixture_residual = mixture_volume - mixture_boundary + print('TFM global dispersed mass residual: {:.6e}'.format(dispersed_residual)) + print('TFM global mixture mass residual: {:.6e}'.format(mixture_residual)) + + assert dispersed_residual == pytest.approx(0.0, abs=1e-11) + assert mixture_residual == pytest.approx(0.0, abs=1e-11) + + +def test_discrete_all_ime_global_mass_conservation(tmp_path: Path): + ime = ('drag', 'laminar_dispersion', 'virtual_mass', 'lift') + set_active_ime(ime) + mesh = ngs.Mesh('pytests/mesh_files/unit_square_coarse.vol') + for _ in range(2): + mesh.Refine() + mesh_file = tmp_path / 'unit_square_refined.vol' + mesh.ngmesh.Save(str(mesh_file)) + + config = ConfigParser('pytests/full_system/tfm/config') + config.set('MESH', 'filename', str(mesh_file)) + config.set('TFM', 'IME', ime_config(ime)) + config.set('TFM', 'lift_wall_deactivation', 'False') + solver = get_solver_class(config)(get_model_class('TwoFluidModel', False), + config) + + exact = exact_solution(0.0) + for name, component in solver.model.model_components.items(): + solver.model.UIter.components[component].Set(exact[component]) + solver._create_linear_and_bilinear_forms() + solver._create_preconditioners() + solver._assemble() + solver.gfu.vec.data = solver.model.IC.vec + solver._apply_boundary_conditions() + solver.model.linear_solve(solver.a[0], solver.L[0], + solver.preconditioners[0], solver.gfu) + + algebraic_residual = solver.L[0].vec.CreateVector() + algebraic_residual.data = (solver.L[0].vec + - solver.a[0].mat * solver.gfu.vec) + pressure_component = solver.model.model_components['p'] + pressure_dofs = solver.model.fes.Range(pressure_component) + mixture_weak_residual = algebraic_residual[pressure_dofs].Norm() + + uc = solver.gfu.components[solver.model.model_components['u_c']] + ud = solver.gfu.components[solver.model.model_components['u_d']] + ac = solver.gfu.components[solver.model.model_components['alpha_c']] + normal = ngs.specialcf.normal(solver.model.mesh.dim) + + mixture_source = solver.model.f['p'][0] + alpha_source = solver.model.f['alpha_c'][0] + mixture_strong_defect = ( + exact[3] * ngs.div(uc) + _grad_scalar(exact[3]) * uc + + ngs.div(ud) - exact[3] * ngs.div(ud) + - _grad_scalar(exact[3]) * ud - mixture_source + ) + mixture_strong_l2 = math.sqrt(ngs.Integrate( + mixture_strong_defect**2, solver.model.mesh)) + mixture_residual = ( + ngs.Integrate(mixture_source, solver.model.mesh) + - ngs.Integrate((exact[3] * uc + (1 - exact[3]) * ud) * normal, + solver.model.mesh, ngs.BND) + ) + dispersed_residual = ( + ngs.Integrate(-alpha_source, solver.model.mesh) + - ngs.Integrate((1 - ac) * exact[1] * normal, + solver.model.mesh, ngs.BND) + ) + print('TFM discrete all-IME mixture mass residual: {:.6e}'.format( + mixture_residual)) + print('TFM discrete all-IME dispersed mass residual: {:.6e}'.format( + dispersed_residual)) + print('TFM discrete all-IME mixture weak residual: {:.6e}'.format( + mixture_weak_residual)) + print('TFM discrete all-IME mixture strong L2 defect: {:.6e}'.format( + mixture_strong_l2)) + assert mixture_residual == pytest.approx(0.0, abs=1e-10) + assert dispersed_residual == pytest.approx(0.0, abs=1e-10) + assert mixture_weak_residual < 1e-10 + + +@pytest.mark.parametrize('mechanism', + ['drag', 'laminar_dispersion', 'virtual_mass', 'lift']) +def test_ime_action_reaction_conservation(mechanism): + mesh = ngs.Mesh('pytests/mesh_files/unit_square_coarse.vol') + force_d, force_c = _physical_ime_forces(mechanism) + residual = _integrate_vector(force_d + force_c, mesh) + print('TFM {} action-reaction residual: ({:.6e}, {:.6e})'.format( + mechanism, residual[0], residual[1])) + assert math.hypot(*residual) < 1e-12 + + +def test_combined_ime_action_reaction_conservation(): + mesh = ngs.Mesh('pytests/mesh_files/unit_square_coarse.vol') + total = ngs.CoefficientFunction((0.0, 0.0)) + for mechanism in ('drag', 'laminar_dispersion', 'virtual_mass', 'lift'): + force_d, force_c = _physical_ime_forces(mechanism) + total += force_d + force_c + residual = _integrate_vector(total, mesh) + print('TFM combined IME action-reaction residual: ({:.6e}, {:.6e})'.format( + residual[0], residual[1])) + assert math.hypot(*residual) < 1e-12 diff --git a/pytests/full_system/tfm/test_tfm_convergence.py b/pytests/full_system/tfm/test_tfm_convergence.py new file mode 100644 index 0000000..61fc22c --- /dev/null +++ b/pytests/full_system/tfm/test_tfm_convergence.py @@ -0,0 +1,98 @@ +import math +from pathlib import Path + +import ngsolve as ngs +import pytest + +from opencmp.config_functions import ConfigParser +from opencmp.models import get_model_class +from opencmp.solvers import get_solver_class +from import_functions import exact_solution, ime_config, set_active_ime + + +def _l2_errors(solution, model, time): + exact = exact_solution(time) + errors = {} + for name, component in model.model_components.items(): + difference = solution.components[component] - exact[component] + if name == 'p': + difference -= ngs.Integrate(difference, model.mesh) / ngs.Integrate(1, model.mesh) + errors[name] = math.sqrt(ngs.Integrate(difference * difference, model.mesh)) + uc = solution.components[model.model_components['u_c']] + ud = solution.components[model.model_components['u_d']] + ac = exact[3] + grad_ac = ngs.CoefficientFunction((ac.Diff(ngs.x), ac.Diff(ngs.y))) + mixture_defect = ( + ac * ngs.div(uc) + grad_ac * uc + + ngs.div(ud) - ac * ngs.div(ud) - grad_ac * ud + - model.f['p'][0] + ) + errors['mixture_mass'] = math.sqrt(ngs.Integrate( + mixture_defect**2, model.mesh)) + return errors + + +def _solve_about_exact_picard_state(solver): + exact = exact_solution(0.0) + for name, component in solver.model.model_components.items(): + solver.model.UIter.components[component].Set(exact[component]) + + solver._create_linear_and_bilinear_forms() + solver._create_preconditioners() + solver._assemble() + solver.gfu.vec.data = solver.gfu_0_list[0].vec + solver._apply_boundary_conditions() + solver.model.linear_solve(solver.a[0], solver.L[0], + solver.preconditioners[0], solver.gfu) + return solver.gfu + + +@pytest.mark.parametrize('ime', [ + (), + ('drag',), + ('drag', 'laminar_dispersion'), + ('lift',), + ('virtual_mass',), + ('drag', 'laminar_dispersion', 'virtual_mass', 'lift'), +], ids=['no-ime', 'drag', 'drag-dispersion', 'lift', 'virtual-mass', 'all-ime']) +def test_tfm_manufactured_solution_h_convergence(tmp_path: Path, ime) -> None: + set_active_ime(ime) + case_name = '+'.join(ime) if ime else 'no-ime' + errors = [] + element_counts = [] + base_mesh = 'pytests/mesh_files/unit_square_coarse.vol' + + for level in range(1, 5): + mesh = ngs.Mesh(base_mesh) + for _ in range(level): + mesh.Refine() + mesh_file = tmp_path / ('unit_square_refined_{}.vol'.format(level)) + mesh.ngmesh.Save(str(mesh_file)) + + config = ConfigParser('pytests/full_system/tfm/config') + config.set('MESH', 'filename', str(mesh_file)) + config.set('TFM', 'IME', ime_config(ime)) + config.set('TFM', 'lift_wall_deactivation', 'False') + solver = get_solver_class(config)(get_model_class('TwoFluidModel', False), config) + solver.gfu_0_list[0].vec.data = solver.model.IC.vec + solution = _solve_about_exact_picard_state(solver) + errors.append(_l2_errors(solution, solver.model, 0.0)) + element_counts.append(solver.model.mesh.ne) + + for variable in ('u_c', 'u_d', 'p', 'alpha_c', 'mixture_mass'): + assert errors[-1][variable] < errors[0][variable] + rates = [ + math.log(errors[level][variable] / errors[level + 1][variable]) + / math.log(math.sqrt(element_counts[level + 1] / element_counts[level])) + for level in range(len(errors) - 1) + ] + print('TFM MMS [{}] {}: errors={}, orders={}'.format( + case_name, variable, + ['{:.6e}'.format(level[variable]) for level in errors], + ['{:.3f}'.format(rate) for rate in rates])) + # The strong mixture residual contains derivatives of the HDiv + # velocity error and therefore converges one order below the velocity + # L2 error. All solution variables should exceed first order; the + # strong mixture defect should approach its expected first order. + minimum_order = 0.9 if variable == 'mixture_mass' else 1.0 + assert min(rates[-2:]) > minimum_order, (variable, errors, rates) diff --git a/pytests/helpers/test_ngsolve_.py b/pytests/helpers/test_ngsolve_.py new file mode 100644 index 0000000..d86cac5 --- /dev/null +++ b/pytests/helpers/test_ngsolve_.py @@ -0,0 +1,73 @@ +######################################################################################################################## +# 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 # +# . # +######################################################################################################################## + +import ngsolve as ngs +from netgen.csg import unit_cube +from numpy import allclose +from pytest import fixture + +from opencmp.helpers.ngsolve_ import curl_3d + + +@fixture +def cube_mesh() -> ngs.Mesh: + """A coarse unit-cube mesh.""" + + return ngs.Mesh(unit_cube.GenerateMesh(maxh=0.4)) + + +def test_curl_3d(cube_mesh) -> None: + """ + curl_3d against a field whose curl is known analytically. + + F = (x*y + 2*z^2, 3*x^2 - y*z, x*z + 5*y^2) -> curl(F) = (11*y, 3*z, 5*x). + + Every component of the expected curl carries a different coefficient AND a + different variable, so a permuted or sign-flipped gradient index cannot pass + by coincidence. F is quadratic, so an order-2 space represents it exactly and + the comparison holds to machine precision. + """ + + fes = ngs.H1(cube_mesh, order=2, dim=3) + gfu = ngs.GridFunction(fes) + gfu.Set(ngs.CoefficientFunction((ngs.x * ngs.y + 2 * ngs.z ** 2, + 3 * ngs.x ** 2 - ngs.y * ngs.z, + ngs.x * ngs.z + 5 * ngs.y ** 2))) + + expected = ngs.CoefficientFunction((11 * ngs.y, 3 * ngs.z, 5 * ngs.x)) + + for point in [(0.3, 0.4, 0.6), (0.15, 0.85, 0.25), (0.5, 0.5, 0.5)]: + mip = cube_mesh(*point) + assert allclose(curl_3d(gfu)(mip), expected(mip), atol=1e-10) + + +def test_curl_3d_of_a_gradient_is_zero(cube_mesh) -> None: + """curl(grad(phi)) == 0 for any smooth phi -- an identity, not a hand-computed value.""" + + phi = ngs.x ** 2 * ngs.y + ngs.y * ngs.z ** 2 + 3 * ngs.x * ngs.z + + fes = ngs.H1(cube_mesh, order=2, dim=3) + gfu = ngs.GridFunction(fes) + # grad(phi) written out componentwise so the test does not depend on the same + # gradient-flattening convention curl_3d is being tested for. + gfu.Set(ngs.CoefficientFunction((2 * ngs.x * ngs.y + 3 * ngs.z, + ngs.x ** 2 + ngs.z ** 2, + 2 * ngs.y * ngs.z + 3 * ngs.x))) + + for point in [(0.3, 0.4, 0.6), (0.7, 0.2, 0.1)]: + mip = cube_mesh(*point) + assert allclose(curl_3d(gfu)(mip), (0.0, 0.0, 0.0), atol=1e-10) diff --git a/pytests/models/test_tfm.py b/pytests/models/test_tfm.py new file mode 100644 index 0000000..5bf3d16 --- /dev/null +++ b/pytests/models/test_tfm.py @@ -0,0 +1,189 @@ +from pathlib import Path +import math +import re + +import ngsolve as ngs +import pytest + +from opencmp.config_functions import ConfigParser +from opencmp.config_functions.boundary_conditions import BCFunctions +from opencmp.models.tfm import TwoFluidModel +import opencmp.models.tfm as tfm_module + + +def _tfm_from_config(tmp_path: Path, tfm: dict[str, str], other=None) -> TwoFluidModel: + config_path = tmp_path / 'config' + config_path.write_text('') + config = ConfigParser(str(config_path)) + config['TFM'] = tfm + config['OTHER'] = other or {} + model = TwoFluidModel.__new__(TwoFluidModel) + model.config = config + model._pre_init() + return model + + +_ALL_IME = ('drag -> Tomiyama\n' + 'lift -> LegendreMagnaudet\n' + 'virtual_mass -> ConstantCoefficient\n' + 'laminar_dispersion -> ConstantCoefficient') + + +def test_tfm_ime_enables_requested_mechanisms(tmp_path: Path) -> None: + model = _tfm_from_config(tmp_path, { + 'canonical_form': 'C-TFM', + 'IME': _ALL_IME, + 'lift_wall_deactivation': 'True', + 'lift_wall_boundaries': 'wall|bottom', + }) + + assert model.canonical_form == 'C-TFM' + assert model.drag_switch + assert model.VM_switch + assert model.Disp_switch + assert model.Lift_switch + assert model.drag_model == 'Tomiyama' + assert model.lift_model == 'LegendreMagnaudet' + assert model.lift_wall_deactivation + assert model.lift_wall_boundaries == 'wall|bottom' + + +@pytest.mark.parametrize(('tfm', 'message'), [ + ({'IME': 'drag -> Tomiyama\nbuoyancy -> Tomiyama'}, 'Unknown [TFM] IME mechanism'), + ({'IME': 'laminar_dispersion -> ConstantCoefficient'}, + "'laminar_dispersion' requires 'drag'"), + ({'canonical_form': 'invalid'}, 'canonical_form'), + ({'IME': 'drag -> invalid'}, "IME 'drag' must use one of"), + ({'IME': 'lift -> invalid'}, "IME 'lift' must use one of"), + ({'IME': 'virtual_mass -> Tomiyama'}, "IME 'virtual_mass' must use one of"), + ({'drag_model': 'Tomiyama'}, 'Unknown [TFM] option'), + ({'unexpected': 'value'}, 'Unknown [TFM] option'), + ({'IME': 'drag -> Tomiyama', 'lift_wall_deactivation': 'True'}, + 'lift_wall_deactivation requires lift'), +]) +def test_tfm_configuration_rejects_invalid_values( + tmp_path: Path, tfm: dict[str, str], message: str) -> None: + with pytest.raises(ValueError, match=re.escape(message)): + _tfm_from_config(tmp_path, tfm) + + +def test_slip_accepts_marker_only_syntax(tmp_path: Path) -> None: + bc_path = tmp_path / 'bc_config' + bc_path.write_text('[SLIP]\nu_d = wall|bottom\n') + + functions = BCFunctions(str(bc_path), str(tmp_path), None, + ['dirichlet', 'neumann', 'slip']) + boundary_conditions, dirichlet_names = functions.set_boundary_conditions( + ['dirichlet', 'neumann', 'slip']) + + assert boundary_conditions['slip']['u_d'] == {'wall': [], 'bottom': []} + assert 'u_d' not in dirichlet_names + + +def test_hdiv_space_combines_dirichlet_and_slip_markers(monkeypatch) -> None: + model = TwoFluidModel.__new__(TwoFluidModel) + model.element = {'u_c': 'HDiv', 'u_d': 'HDiv', 'p': 'L2', 'alpha_c': 'L2'} + model.interp_ord = 2 + model.mesh = object() + model.DG = True + model.BC = { + 'dirichlet': {'u_c': {'inlet': []}, 'u_d': {'outlet': []}}, + 'slip': {'u_c': {}, 'u_d': {'wall': [], 'bottom': []}}, + } + + hdiv_calls = [] + + def fake_hdiv(mesh, **kwargs): + hdiv_calls.append(kwargs) + return object() + + monkeypatch.setattr(tfm_module.ngs, 'HDiv', fake_hdiv) + monkeypatch.setattr(tfm_module.ngs, 'L2', lambda mesh, **kwargs: object()) + monkeypatch.setattr(tfm_module, 'FESpace', lambda spaces, **kwargs: spaces) + + model._construct_fes() + + assert hdiv_calls[0]['dirichlet'] == 'inlet' + assert hdiv_calls[1]['dirichlet'] == 'outlet|wall|bottom' + + +def test_dirichlet_and_slip_cannot_overlap() -> None: + model = TwoFluidModel.__new__(TwoFluidModel) + model.BC = { + 'dirichlet': {'u_d': {'wall': []}}, + 'slip': {'u_d': {'wall': []}}, + } + + with pytest.raises(ValueError, match='both DIRICHLET and SLIP'): + model._velocity_space_boundaries('u_d') + + +def test_hdiv_slip_constrains_normal_but_preserves_tangential_velocity() -> None: + mesh = ngs.Mesh('pytests/mesh_files/unit_square_coarse.vol') + for _ in range(2): + mesh.Refine() + + # In HDiv, marking the bottom as Dirichlet constrains only u.n. Project a + # horizontal field with a nonzero tangential component toward that wall. + fes = ngs.HDiv(mesh, order=2, dirichlet='bottom') + u, v = fes.TnT() + a = ngs.BilinearForm(fes) + a += (u * v + 1e-3 * ngs.div(u) * ngs.div(v)) * ngs.dx + L = ngs.LinearForm(fes) + target = ngs.CoefficientFunction((1.0, 0.0)) + L += target * v * ngs.dx + a.Assemble() + L.Assemble() + + velocity = ngs.GridFunction(fes) + velocity.vec.data = a.mat.Inverse( + freedofs=fes.FreeDofs(), inverse='umfpack') * L.vec + + normal = ngs.specialcf.normal(mesh.dim) + normal_l2_squared = ngs.Integrate( + (velocity * normal)**2, mesh, ngs.BND, + definedon=mesh.Boundaries('bottom')) + near_wall_tangential_energy = ngs.Integrate( + ngs.IfPos(0.125 - ngs.y, 1.0, 0.0) * velocity[0]**2, mesh) + projection_error = ngs.Integrate((velocity - target)**2, mesh) + + print('TFM SLIP bottom normal L2^2: {:.6e}'.format(normal_l2_squared)) + print('TFM SLIP near-wall tangential energy: {:.6e}'.format( + near_wall_tangential_energy)) + assert normal_l2_squared == pytest.approx(0.0, abs=1e-14) + assert near_wall_tangential_energy > 0.1 + assert math.sqrt(projection_error) < 1e-12 + + +def test_ime_closures_are_finite_and_zero_without_a_driver() -> None: + config = ConfigParser('pytests/full_system/tfm/config') + config.set('TFM', 'IME', _ALL_IME) + config.set('TFM', 'lift_wall_deactivation', 'False') + model = TwoFluidModel(config, [ngs.Parameter(0.0)]) + + velocity = ngs.CoefficientFunction((1.0, 0.0)) + for alpha_c in (0.1, 0.5, 0.9): + model.UIter.components[model.model_components['u_c']].Set(velocity) + model.UIter.components[model.model_components['u_d']].Set(velocity) + model.UIter.components[model.model_components['alpha_c']].Set(alpha_c) + + wc = model.UIter.components[model.model_components['u_c']] + wd = model.UIter.components[model.model_components['u_d']] + ac = model.UIter.components[model.model_components['alpha_c']] + ad = 1 - ac + cd = model._get_drag_coeff(wd, wc, ad, 0) + cl = model._get_lift_coeff(wd, wc, 0) + + # With no slip, drag and dispersion have zero magnitude; lift also has + # zero slip and zero carrier vorticity. Virtual mass is covered by its + # MMS convergence and action-reaction regression tests. + drag = cd * ngs.Norm(wd - wc) * (wd - wc) + dispersion = cd * ngs.Norm(wd - wc)**2 * ngs.grad(ac) + curl_wc = ngs.grad(wc)[1] - ngs.grad(wc)[2] + lift = cl * curl_wc * ngs.CoefficientFunction( + (wd[1] - wc[1], -(wd[0] - wc[0]))) + + for closure in (drag, dispersion, lift): + magnitude_squared = ngs.Integrate(closure * closure, model.mesh) + assert math.isfinite(magnitude_squared) + assert magnitude_squared == pytest.approx(0.0, abs=1e-24)