sundials-burgers.hs

{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedLists     #-}
{-# LANGUAGE NumDecimals         #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE DataKinds           #-}

{-# OPTIONS_GHC -Wall            #-}

module Main (main) where

import           Numeric.Sundials
import           Numeric.LinearAlgebra
import           Data.Csv
import           Data.Char
import           Data.ByteString.Lazy (putStr, writeFile)
import           Prelude hiding (putStr, writeFile)
import           Control.Exception
import           Data.Coerce
import           Katip.Monadic
import           GHC.Int


bigN :: Int
bigN = 101 -- 201

deltaX :: Double
deltaX = 1.0 / (fromIntegral bigN - 1)

bigC :: Matrix Double
bigC = assoc (bigN, bigN) 0.0 [ ((i, j), (1/(2*deltaX)) * f (i, j)) | i <- [0 .. bigN - 1]
                                                                    , j <- [0 .. bigN - 1]
                        ]
 where
   f (i, j)
            | i     == 0        =  0
            | i     == bigN - 1 =  0
            | i     == j        =  0
            | i - 1 == j        = -1
            | i + 1 == j        =  1
            -- | i     == 0 &&
            --   j     == bigN - 1 = -1
            -- | i     == bigN - 1  &&
            --   j     == 0        =  1
            | otherwise         =  0

bigD :: Matrix Double
bigD = assoc (bigN, bigN) 0.0 [ ((i, j), f (i, j)) | i <- [0 .. bigN - 1]
                                                   , j <- [0 .. bigN - 1]
                        ]
 where
   f (i, j) | i     == j        =  1
            | otherwise         =  0

bigV :: Vector Double
bigV = assoc bigN 0.0 [(i, f i) | i <- [0 .. bigN - 1]]
  where
  f i = (sin (2 * pi * fromIntegral i * deltaX))

sol' :: IO (Matrix Double)
sol' = do
  w <- runNoLoggingT $ solve (defaultOpts' HEUN_EULER_2_1_2) burgers
  case w of
    Left e  -> error $ show e
    Right y -> return (solutionMatrix y)

burgers :: OdeProblem
burgers = emptyOdeProblem
  { odeRhs = odeRhsPure $ \_t x -> coerce (cmap negate ((bigD #> (coerce x)) * (bigC #> (coerce x))))
  , odeJacobian = Nothing
  , odeEvents = mempty
  , odeEventHandler = nilEventHandler
  , odeMaxEvents = 0
  , odeInitCond = bigV
  , odeSolTimes = vector $ map (0.025 *) [0 .. 10]
  , odeTolerances = defaultTolerances
  }

myOptions :: EncodeOptions
myOptions = defaultEncodeOptions {
      encDelimiter = fromIntegral (ord ' ')
    }

main :: IO ()
main = do
  y <- sol'
  writeFile "burgers.txt" $ encodeWith myOptions $ map toList $ toRows y


defaultOpts' :: method -> ODEOpts method
defaultOpts' method = ODEOpts
  { maxNumSteps = 1e5
  , minStep     = 1.0e-14
  , fixedStep   = 0
  , maxFail     = 10
  , odeMethod   = method
  , initStep    = Nothing
  , jacobianRepr = DenseJacobian
  }

instance Element Int8

emptyOdeProblem :: OdeProblem
emptyOdeProblem = OdeProblem
      { odeRhs = error "emptyOdeProblem: no odeRhs provided"
      , odeJacobian = Nothing
      , odeInitCond = error "emptyOdeProblem: no odeInitCond provided"
      , odeEvents = mempty
      , odeTimeBasedEvents = TimeEventSpec $ return $ 1.0 / 0.0
      , odeEventHandler = nilEventHandler
      , odeMaxEvents = 100
      , odeSolTimes = error "emptyOdeProblem: no odeSolTimes provided"
      , odeTolerances = defaultTolerances
      }

nilEventHandler :: EventHandler
nilEventHandler _ _ _ = throwIO $ ErrorCall "nilEventHandler"

defaultTolerances :: Tolerances
defaultTolerances = Tolerances
  { absTolerances = Left 1.0e-6
  , relTolerance = 1.0e-10
  }