{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies\n", " , UnicodeSyntax, ScopedTypeVariables, NoMonomorphismRestriction #-}\n", "import Prelude ()\n", "import Data.Manifold.TreeCover\n", "import Data.Manifold.Shade\n", "import Data.Random\n", "import Data.Random.Manifold\n", "import Data.Manifold\n", "import Data.Manifold.Web\n", "import Data.Manifold.DifferentialEquation\n", "import Math.LinearMap.Category\n", "import Math.LinearMap.Category.Derivatives\n", "import Data.VectorSpace\n", "import Data.Function.Affine\n", "import Data.Basis (Basis)\n", "import Linear(V2(..), ex, ey, _x, _y)\n", "import Data.Semigroup\n", "import Data.List.NonEmpty (NonEmpty(..))\n", "import qualified Data.Foldable as Hask\n", "import qualified Control.Monad as Hask\n", "import qualified Control.Comonad.Cofree as Hask\n", "import Control.Lens\n", ":opt no-lint\n", "import Control.Category.Constrained.Prelude\n", "import Control.Arrow.Constrained\n", "test x = x `seq` putStrLn \"ok\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From [diagrams](http://projects.haskell.org/diagrams):" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import Diagrams.Prelude (p2, circle, (&), (^&), moveTo, opacity, orange, fromVertices, Point(P))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "type X = ℝ\n", "type T = ℝ\n", "type U = ℝ\n", "type Ðx'U = ℝ\n", "type Ðt'U = ℝ\n", "type Ðx²'U = ℝ\n", "type x × y = ℝ²\n", "type HeatFlow = (U, Ðx'U)\n", "et = ey :: Basis ℝ²\n", "_t = _y :: Lens' (X × T) T\n", "_u = _1 :: Lens' (U, Ðx'U) U\n", "_ðx'u = _2 :: Lens' (U, Ðx'U) Ðx'U" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plotting\n", "---\n", "\n", "From [dynamic-plot](http://hackage.haskell.org/package/dynamic-plot):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "ok" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "ok" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import Graphics.Dynamic.Plot.R2\n", "import Data.Colour.Names\n", "import Data.Colour.Manifold\n", "\n", "type Colourscheme y = Shade' y -> Shade (Colour ℝ)\n", "\n", "colourscheme :: Colourscheme ℝ\n", "colourscheme = f `seq` \\(Shade' y e) -> f (Shade y $ dualNorm e)\n", " where Just (f :: Shade ℝ -> Shade (Colour ℝ))\n", " = rangeWithinVertices (0, neutral)\n", " [ (2, red) ]\n", " Just neutral = toInterior (dimgrey :: Colour ℝ)\n", "test colourscheme\n", "prettyWebPlot :: PointsWeb ℝ² y -> DynamicPlottable\n", "prettyWebPlot w = plot [ diagramPlot . opacity 0.1 $ fromVertices [P r₁, P r₂]\n", " | ((r₁@(V2 x₁ y₁),_),(r₂@(V2 x₂ y₂),_)) <- edg ]\n", " where edg = webEdges w\n", " \n", "colourscheme_heat :: Colourscheme U\n", "colourscheme_heat = f `seq` \\(Shade' y e) -> f (Shade y $ dualNorm e)\n", " where Just (f :: Shade U -> Shade (Colour ℝ))\n", " = rangeWithinVertices (0, neutral)\n", " [ (1, orange)\n", " {-, ((0,2), seagreen)-} ]\n", " Just neutral = toInterior (dimgrey :: Colour ℝ)\n", "test colourscheme_heat" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [], "source": [ "uncertainFieldSeqPlot :: (Hask.Functor m, Hask.Foldable m)\n", " => Colourscheme y -> (m () -> DynamicPlottable)\n", " -> Hask.Cofree m (PointsWeb ℝ² (Shade' y)) -> DynamicPlottable\n", "uncertainFieldSeqPlot cscm errDisp = plotLatest . listd 0\n", " where listd i ~(a Hask.:< q) = step : case Hask.toList q of\n", " [] -> [errDisp (const () <$> q) <> step]\n", " l -> foldr1 (>>) $ listd (i+1)<$>l\n", " where step = plot (fmap cscm a) & legendName (\"i = \"++show i)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import Control.Monad.Trans.Except\n", "pdeApproachPlot :: Colourscheme y\n", " -> Hask.Cofree (Except (PropagationInconsistency (X × T) (Shade' y)))\n", " (PointsWeb ℝ² (Shade' y))\n", " -> DynamicPlottable\n", "pdeApproachPlot cscm = uncertainFieldSeqPlot cscm $ \\e -> case runExcept e of\n", " Left (PropagationInconsistency pro _)\n", " -> shapePlot (Hask.fold [circle 0.02 & moveTo (P p) | (p,_) <- pro]) & tint orange" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "trimDecimals :: ℝ -> ℝ\n", "trimDecimals = (/q) . fromInteger . round . (*q)\n", " where q = 1000\n", "\n", "import Control.Applicative ((<|>))\n", "honestStrategy, cheatingStrategy :: ∀ x y. (Refinable y, y ~ Interior y)\n", " => InformationMergeStrategy [] (Except (PropagationInconsistency x (Shade' y))) (x, Shade' y) (Shade' y)\n", "honestStrategy = indicateInconsistencies intersectShade's\n", "cheatingStrategy = indicateInconsistencies (\\l -> intersectShade's l <|> mixShade's l)\n", "\n", "import Control.Monad.Trans.Writer\n", "stubbornStrategy :: ∀ x y. (Refinable y, y ~ Interior y)\n", " => InformationMergeStrategy [] (WriterT [PropagationInconsistency x (Shade' y)] Maybe) (x, Shade' y) (Shade' y)\n", "stubbornStrategy = postponeInconsistencies intersectShade's\n", "\n", "veryVague :: (LSpace (Needle x), Fractional' (Scalar (Needle x)))\n", " => Shade' x -> Shade' x\n", "veryVague (Shade' x ex) = Shade' x $ scaleNorm 0.01 ex" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Grid generation\n", "---" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "rectangular, perturbedHexagonal :: (x~ℝ, y~ℝ)\n", " => ((x,x),Needle x) -> ((y,y),Needle y) -> [ℝ²]\n", "\n", "rectangular ((xl,xr),δx) ((yb,yt),δy)\n", " = [V2 x y | x<-[xl-δx/2, xl+δx/2 .. xr], y<-[yb-δy/2, yb+δy/2 .. yt]]\n", "\n", "perturbedHexagonal ((xl,xr),δx) ((yb,yt),δy)\n", " = [ V2 x y\n", " | y'<-trimDecimals<$>[yb - δy/2, yb+δy/2 .. yt]\n", " , let xlq = xl + skew*(y' - yt)\n", " , x<-trimDecimals<$>[xlq-δx, xlq .. xr+δx]\n", " , x >= xl-δx\n", " , let y = trimDecimals $ y'\n", " + 0.005 * sin (37245*sin(2334*x)+y') -- pseudorandom pertubation\n", " ]\n", " where skew = 8/14" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Advection\n", "===\n", "\n", "Setup and analytic solution\n", "---" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "β = 0.5\n", "(xl,xr) = (0,4)\n", "\n", "exact_advect :: X -> T -> U\n", "exact_advect x t = u 2\n", " where u n = sin (n*pi*(x - xl + β*t)/l)\n", " l = xr - xl\n", "\n", "--plotWindow [ continFnPlot (`exact_advect`t) & legendName (\"t = \"++show t)\n", " -- | t<-[0,0.2 .. 1] ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 1D linear advection equation\n", "---\n", "$$\n", " \\frac{\\partial u}{\\partial t} = \\beta \\cdot \\frac{\\partial u}{\\partial x}\n", "$$\n", "with Dirichlet right boundary condition\n", "$$\n", " u(x_r) = 0\n", "$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "advectEqn :: DifferentialEqn (X × T) U\n", "advectEqn (Shade (V2 x t, _) _) = LocalDifferentialEqn {\n", " _rescanDifferentialEqn = \\d⁰ d¹ d²\n", " -> let dx¹ = projectShade (lensEmbedding (1*∂id/∂_x)) d¹ :: Shade' Ðx'U\n", " dt¹ = projectShade (lensEmbedding (1*∂id/∂_t)) d¹ :: Shade' Ðx'U\n", " in ( return d⁰\n", " , mixShade's $ embedShade (lensEmbedding (recip β*∂id/∂_t)) dx¹\n", " :|[veryVague $ embedShade (lensEmbedding (β*∂id/∂_x)) dt¹] )\n", " }" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "startSt_advect = fromWebNodes euclideanMetric\n", " [ (V2 x t, if t<0\n", " then exact_advect x t|±|[1e-2]\n", " else 0|±|[10])\n", " | V2 x t <- rectangular ((xl,xr),δx₀) ((tStart,tEnd),δt₀) ]\n", " where (tStart, tEnd) = (0, 4)\n", " δx₀ = 1/3\n", " δt₀ = δx₀\n", "\n", "tfs_advect = iterateFilterDEqn_static_selective honestStrategy id (euclideanVolGoal 0.1) advectEqn startSt_advect\n", "--putStrLn \"Precalculate...\"\n", "--print . length . take 30 $ Hask.toList tfs_advect" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "GraphWindowSpecR2{lBound=-0.8079137577969195, rBound=4.41670473388506, bBound=-1.0778231696181226, tBound=4.336401908145449, xResolution=640, yResolution=480}" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plotWindow [ prettyWebPlot $ head (Hask.toList tfs_advect)\n", " , pdeApproachPlot colourscheme tfs_advect\n", " , dynamicAxes ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![1D linear advection equation solved](https://raw.githubusercontent.com/leftaroundabout/manifolds/master/manifolds/images/examples/PDE-solution-filter/Advection-simple.gif)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotWindow [plotLatest $\n", " [ plotMultiple [ plot (fromWebNodes euclideanMetric . map (first (^._x))\n", " $ sliceWeb_lin tf (normalPlane (V2 0 t) (V2 0 1)))\n", " & legendName (\"t = \"++show t)\n", " | t <- [0,0.5..2.5] ] & legendName (\"i = \"++show i)\n", " | (i,tf) <- zip [0..60] $ Hask.toList tfs_advect ] ] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Waves\n", "===\n", "\n", "Setup and analytic solution\n", "---" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c = 2\n", "(xl,xr) = (0,2)\n", "\n", "exact_wave :: X -> T -> (U, Ðx'U)\n", "exact_wave x t = u 1\n", " where u n = ( sin (n*pi*(x - xl - c*t)/l), n*pi/l*cos (n*pi*(x - xl - c*t)/l) )\n", " ^+^ ( sin (n*pi*(x - xl + c*t)/l), n*pi/l*cos (n*pi*(x - xl + c*t)/l) )\n", " l = xr - xl\n", "\n", "plotWindow [plotLatest [ continFnPlot (fst . (`exact_wave`t)) & legendName (\"t = \"++show t)\n", " | t<-[0,0.05 .. ] ] ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The 1D linear wave equation\n", "---\n", "$$\n", " \\frac{\\partial^2 u}{\\partial t^2} = c^2 \\cdot \\frac{\\partial^2 u}{\\partial x^2}\n", "$$\n", "with Dirichlet boundary condition\n", "$$\n", " u|_{\\partial\\Omega} = 0\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "waveEqn :: DifferentialEqn (X × T) ((Ðx'U,Ðt'U),Ðx²'U) (U, Ðx'U)\n", "waveEqn (Shade (V2 x t, _) _) = LocalDifferentialEqn {\n", " _predictDerivatives = factoriseShade >>> first factoriseShade >>> \\((dx¹,dt¹),dx²)\n", " -> mixShade's $ embedShade (lensEmbedding (recip (c^2)*∂id/∂_t)) dx²\n", " :|[veryVague $ embedShade (lensEmbedding (1*∂id/∂_x)) dx¹]\n", " , _rescanDerivatives = \\d⁰ d¹ d²\n", " -> let dx¹ = projectShade (lensEmbedding (1*∂id/∂_x)) d¹\n", " in if t<0 || x<0 || x>xr\n", " then ( return $ exact_wave x 0|±|[1e-6]\n", " , return $ dx¹ ✠ projectShade (lensEmbedding (recip α*∂id/∂_t)) d¹ )\n", " else ( return $ d⁰\n", " , return $ dx¹ ✠ projectShade (lensEmbedding (1*∂id/∂_x.∂_x)) d²\n", " )\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Heat equation\n", "===\n", "\n", "Setup and analytic solution\n", "---" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "α = 0.1\n", "(xl,xr) = (0,2)\n", "\n", "exact_heat :: X -> T -> U\n", "exact_heat x t = u 1 - u 3/3\n", " where u n = sin (n*pi*(x-xl)/l) * exp(-α * (n*pi/l)^2 * t)\n", " l = xr - xl\n", "\n", "--plotWindow [ continFnPlot (`exact_heat`t) & legendName (\"t = \"++show t)\n", " -- | t<-[0,0.2 .. 1] ]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The heat equation\n", "---\n", "$$\n", " \\frac{\\partial u}{\\partial t} = \\alpha \\cdot \\frac{\\partial^2u}{\\partial x^2}\n", "$$\n", "With homogeneous Dirichlet boundary conditions\n", "$$\n", " u |_{\\partial\\Omega} = 0\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "heatEqn :: DifferentialEqn (X × T) (Ðx'U,Ðx²'U) U\n", "heatEqn (Shade (V2 x t, _) _) = LocalDifferentialEqn {\n", " _predictDerivatives = factoriseShade >>> \\(dx¹,dx²)\n", " -> mixShade's $ embedShade (lensEmbedding (recip α*∂id/∂_t)) dx²\n", " :|[veryVague $ embedShade (lensEmbedding (1*∂id/∂_x)) dx¹]\n", " , _rescanDerivatives = \\d⁰ d¹ d²\n", " -> let dx¹ = projectShade (lensEmbedding (1*∂id/∂_x)) d¹\n", " in if t<0 || x<0 || x>xr\n", " then ( return $ exact_heat x 0|±|[1e-6]\n", " , return $ dx¹ ✠ projectShade (lensEmbedding (recip α*∂id/∂_t)) d¹ )\n", " else ( return $ d⁰\n", " , return $ dx¹ ✠ projectShade (lensEmbedding (1*∂id/∂_x.∂_x)) d²\n", " )\n", " }" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "startSt_heat = fromWebNodes euclideanMetric\n", " [ (V2 x t, 0|±|[2])\n", " | V2 x t <- cubic ((xl,xr),δx₀) ((tStart,tEnd),δt₀) ]\n", " where (tStart, tEnd) = (0, 1)\n", " δx₀ = 1/8\n", " δt₀ = δx₀\n", "\n", "tfs_heat = iterateFilterDEqn_static honestStrategy id heatEqn startSt_heat\n", "forM_ [Hask.toList tfs_heat ]\n", " $ \\tfs ->\n", " plotWindow\n", " [ prettyWebPlot $ head tfs\n", " , plotLatest [ plot (fmap colourscheme_heat tfi)\n", " & legendName (\"i = \"++show i)\n", " | (i,tfi) <- zip [0..] tfs ]\n", " , dynamicAxes ]\n", "Hask.mapM_ (\\st -> mapM_ (\\((p,y),(q,z)) -> putStrLn $\n", " show p ++ \":\\t\" ++ prettyShowShade' y\n", " ++ \"\\n\" ++ show q ++ \":\\t\" ++ prettyShowShade' z) (webEdges st) >> print())\n", " $ take 0 (Hask.toList tfs_heat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Examination of inconsistencies\n", "===" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pertubation :: X -> T -> U\n", "pertubation x t = 0.4 * exp (-((x-0.7)^2 + (t-0.6)^2)*3)\n", "\n", "startSt_perturbed = fromWebNodes euclideanMetric\n", " [ (V2 x t, (exact_heat x t + pertubation x t, 0)\n", " |±|[(0.01,0), (0,0.10)])\n", " | x<-[xl,xl+δx₀ .. xr], t<-[tStart, tStart+δt₀ .. tEnd]]\n", " where δx₀ = 0.0625\n", " δt₀ = 0.0625\n", "\n", "tfs_heat = iterateFilterDEqn_static strategy id heatEqn startSt_perturbed\n", "forM_ [Hask.toList tfs_heat ]\n", " $ \\tfs ->\n", " plotWindow\n", " [ prettyWebPlot $ head tfs\n", " , plotLatest [ plot (fmap colourscheme_heat tfi)\n", " & legendName (\"i = \"++show i)\n", " | (i,tfi) <- zip [0..] tfs ] ]\n", "Hask.mapM_ (\\st -> mapM_ (\\((p,y),(q,z)) -> putStrLn $\n", " show p ++ \":\\t\" ++ prettyShowShade' y\n", " ++ \"\\n\" ++ show q ++ \":\\t\" ++ prettyShowShade' z) (webEdges st) >> print())\n", " $ take 0 (Hask.toList tfs_heat)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tf_ð²x'U :: PointsWeb ℝ² (Shade' ℝ)\n", "Just tf_ð²x'U = fmap (snd . factoriseShade . snd)\n", " <$> rescanPDEOnWeb AbortOnInconsistency heatEqn startSt_heat\n", "length $ Hask.toList tf_ð²x'U\n", "-- Hask.toList $ tf_ð²x'U\n", "plotWindow [ prettyWebPlot $ tf_ð²x'U\n", " , plot $ colourscheme <$> tf_ð²x'U\n", " , dynamicAxes ]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "forM_ [Hask.toList tfs_heat] $ \\(tf:_) -> do\n", " let Just tfð = rescanPDEOnWeb AbortOnInconsistency heatEqn tf\n", " print . length $ Hask.toList tfð\n", " plotWindow [ prettyWebPlot $ tf\n", " , plot $ colourscheme . snd . factoriseShade . snd <$> tfð\n", " , dynamicAxes ]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "findErr :: (Hask.Monad m, Hask.Foldable m) => Hask.Cofree m a -> m ()\n", "findErr (a Hask.:< q) = case Hask.toList q of\n", " [] -> const () <$> q\n", " l -> foldr1 (>>) $ findErr<$>l\n", "\n", "case runExcept $ findErr tfs_advect of\n", " Left (PropagationInconsistency pro apr) -> do\n", " putStrLn $ \"apriori: \"++show apr\n", " forM_ pro $ \\(o,m) -> putStrLn $ show o ++ \"\\t-> \" ++ show m\n", " --plotWindow $ (plot apr & legendName \"apriori\")\n", " -- : [plot m & legendName (\"from \"++show o) | (o,m)<-pro]\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ ":t tfs_advect" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Haskell", "language": "haskell", "name": "haskell" }, "language_info": { "codemirror_mode": "ihaskell", "file_extension": ".hs", "name": "haskell", "version": "8.2.1" } }, "nbformat": 4, "nbformat_minor": 1 }