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