module Simulation.Aivika.Internal.Specs
(Specs(..),
Method(..),
Run(..),
Point(..),
EventPriority(..),
EventQueue(..),
newEventQueue,
basicTime,
integIterationBnds,
integIterationHiBnd,
integIterationLoBnd,
integPhaseBnds,
integPhaseHiBnd,
integPhaseLoBnd,
integTimes,
integPoints,
integPointsStartingFrom,
integStartPoint,
integStopPoint,
simulationStopPoint,
timeGrid,
pointAt,
delayPoint) where
import Data.IORef
import Simulation.Aivika.Generator
import qualified Simulation.Aivika.PriorityQueue.EventQueue as PQ
data Specs = Specs { Specs -> Double
spcStartTime :: Double,
Specs -> Double
spcStopTime :: Double,
Specs -> Double
spcDT :: Double,
Specs -> Method
spcMethod :: Method,
Specs -> GeneratorType
spcGeneratorType :: GeneratorType
}
data Method = Euler
| RungeKutta2
| RungeKutta4
| RungeKutta4b
deriving (Method -> Method -> Bool
(Method -> Method -> Bool)
-> (Method -> Method -> Bool) -> Eq Method
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Method -> Method -> Bool
== :: Method -> Method -> Bool
$c/= :: Method -> Method -> Bool
/= :: Method -> Method -> Bool
Eq, Eq Method
Eq Method =>
(Method -> Method -> Ordering)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Bool)
-> (Method -> Method -> Method)
-> (Method -> Method -> Method)
-> Ord Method
Method -> Method -> Bool
Method -> Method -> Ordering
Method -> Method -> Method
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Method -> Method -> Ordering
compare :: Method -> Method -> Ordering
$c< :: Method -> Method -> Bool
< :: Method -> Method -> Bool
$c<= :: Method -> Method -> Bool
<= :: Method -> Method -> Bool
$c> :: Method -> Method -> Bool
> :: Method -> Method -> Bool
$c>= :: Method -> Method -> Bool
>= :: Method -> Method -> Bool
$cmax :: Method -> Method -> Method
max :: Method -> Method -> Method
$cmin :: Method -> Method -> Method
min :: Method -> Method -> Method
Ord, Int -> Method -> ShowS
[Method] -> ShowS
Method -> String
(Int -> Method -> ShowS)
-> (Method -> String) -> ([Method] -> ShowS) -> Show Method
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Method -> ShowS
showsPrec :: Int -> Method -> ShowS
$cshow :: Method -> String
show :: Method -> String
$cshowList :: [Method] -> ShowS
showList :: [Method] -> ShowS
Show)
data Run = Run { Run -> Specs
runSpecs :: Specs,
Run -> Int
runIndex :: Int,
Run -> Int
runCount :: Int,
Run -> EventQueue
runEventQueue :: EventQueue,
Run -> Generator
runGenerator :: Generator
}
data Point = Point { Point -> Specs
pointSpecs :: Specs,
Point -> Run
pointRun :: Run,
Point -> Double
pointTime :: Double,
Point -> Int
pointPriority :: EventPriority,
Point -> Int
pointIteration :: Int,
Point -> Int
pointPhase :: Int
}
type EventPriority = Int
data EventQueue = EventQueue { EventQueue -> PriorityQueue (Point -> IO ())
queuePQ :: PQ.PriorityQueue (Point -> IO ()),
EventQueue -> IORef Bool
queueBusy :: IORef Bool,
EventQueue -> IORef Double
queueTime :: IORef Double
}
newEventQueue :: Specs -> IO EventQueue
newEventQueue :: Specs -> IO EventQueue
newEventQueue Specs
specs =
do IORef Bool
f <- Bool -> IO (IORef Bool)
forall a. a -> IO (IORef a)
newIORef Bool
False
IORef Double
t <- Double -> IO (IORef Double)
forall a. a -> IO (IORef a)
newIORef (Double -> IO (IORef Double)) -> Double -> IO (IORef Double)
forall a b. (a -> b) -> a -> b
$ Specs -> Double
spcStartTime Specs
specs
PriorityQueue (Point -> IO ())
pq <- IO (PriorityQueue (Point -> IO ()))
forall a. IO (PriorityQueue a)
PQ.newQueue
EventQueue -> IO EventQueue
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return EventQueue { queuePQ :: PriorityQueue (Point -> IO ())
queuePQ = PriorityQueue (Point -> IO ())
pq,
queueBusy :: IORef Bool
queueBusy = IORef Bool
f,
queueTime :: IORef Double
queueTime = IORef Double
t }
integIterations :: Specs -> [Int]
integIterations :: Specs -> [Int]
integIterations Specs
sc = [Int
i1 .. Int
i2] where
i1 :: Int
i1 = Specs -> Int
integIterationLoBnd Specs
sc
i2 :: Int
i2 = Specs -> Int
integIterationHiBnd Specs
sc
integIterationBnds :: Specs -> (Int, Int)
integIterationBnds :: Specs -> (Int, Int)
integIterationBnds Specs
sc = (Int
i1, Int
i2) where
i1 :: Int
i1 = Specs -> Int
integIterationLoBnd Specs
sc
i2 :: Int
i2 = Specs -> Int
integIterationHiBnd Specs
sc
integIterationLoBnd :: Specs -> Int
integIterationLoBnd :: Specs -> Int
integIterationLoBnd Specs
sc = Int
0
integIterationHiBnd :: Specs -> Int
integIterationHiBnd :: Specs -> Int
integIterationHiBnd Specs
sc =
let n :: Int
n = Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
round ((Specs -> Double
spcStopTime Specs
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
-
Specs -> Double
spcStartTime Specs
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc)
in if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
then
String -> Int
forall a. HasCallStack => String -> a
error (String -> Int) -> String -> Int
forall a b. (a -> b) -> a -> b
$
String
"The iteration number in the stop time has a negative value. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"Either the simulation specs are incorrect, " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"or a floating point overflow occurred, " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"for example, when using a too small integration time step. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"You have to define this time step regardless of " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"whether you actually use it or not, " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"for Aivika allows combining the ordinary differential equations " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"with the discrete event simulation within one model. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"So, if you are still using the 32-bit architecture and " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"you do need a small integration time step " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"for integrating the equations " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"then you might think of using the 64-bit architecture. " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"Although you could probably just forget " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"to increase the time step " String -> ShowS
forall a. [a] -> [a] -> [a]
++
String
"after increasing the stop time: integIterationHiBnd"
else Int
n
integPhases :: Specs -> [Int]
integPhases :: Specs -> [Int]
integPhases Specs
sc =
case Specs -> Method
spcMethod Specs
sc of
Method
Euler -> [Int
0]
Method
RungeKutta2 -> [Int
0, Int
1]
Method
RungeKutta4 -> [Int
0, Int
1, Int
2, Int
3]
Method
RungeKutta4b -> [Int
0, Int
1, Int
2, Int
3]
integPhaseBnds :: Specs -> (Int, Int)
integPhaseBnds :: Specs -> (Int, Int)
integPhaseBnds Specs
sc =
case Specs -> Method
spcMethod Specs
sc of
Method
Euler -> (Int
0, Int
0)
Method
RungeKutta2 -> (Int
0, Int
1)
Method
RungeKutta4 -> (Int
0, Int
3)
Method
RungeKutta4b -> (Int
0, Int
3)
integPhaseLoBnd :: Specs -> Int
integPhaseLoBnd :: Specs -> Int
integPhaseLoBnd Specs
sc = Int
0
integPhaseHiBnd :: Specs -> Int
integPhaseHiBnd :: Specs -> Int
integPhaseHiBnd Specs
sc =
case Specs -> Method
spcMethod Specs
sc of
Method
Euler -> Int
0
Method
RungeKutta2 -> Int
1
Method
RungeKutta4 -> Int
3
Method
RungeKutta4b -> Int
3
basicTime :: Specs -> Int -> Int -> Double
basicTime :: Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
ph =
if Int
ph Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then
String -> Double
forall a. HasCallStack => String -> a
error String
"Incorrect phase: basicTime"
else
Specs -> Double
spcStartTime Specs
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
n' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Method -> Int -> Double
forall {a}. (Eq a, Num a) => Method -> a -> Double
delta (Specs -> Method
spcMethod Specs
sc) Int
ph
where n' :: Double
n' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
delta :: Method -> a -> Double
delta Method
Euler a
0 = Double
0
delta Method
RungeKutta2 a
0 = Double
0
delta Method
RungeKutta2 a
1 = Specs -> Double
spcDT Specs
sc
delta Method
RungeKutta4 a
0 = Double
0
delta Method
RungeKutta4 a
1 = Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
delta Method
RungeKutta4 a
2 = Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
delta Method
RungeKutta4 a
3 = Specs -> Double
spcDT Specs
sc
delta Method
RungeKutta4b a
0 = Double
0
delta Method
RungeKutta4b a
1 = Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
delta Method
RungeKutta4b a
2 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
delta Method
RungeKutta4b a
3 = Specs -> Double
spcDT Specs
sc
integTimes :: Specs -> [Double]
integTimes :: Specs -> [Double]
integTimes Specs
sc = (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Double
t [Int
nl .. Int
nu]
where (Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
t :: Int -> Double
t Int
n = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0
integPoints :: Run -> [Point]
integPoints :: Run -> [Point]
integPoints Run
r = [Point]
points
where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
(Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
points :: [Point]
points = (Int -> Point) -> [Int] -> [Point]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Point
point [Int
nl .. Int
nu]
point :: Int -> Point
point Int
n = Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
pointRun :: Run
pointRun = Run
r,
pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
pointPriority :: Int
pointPriority = Int
0,
pointIteration :: Int
pointIteration = Int
n,
pointPhase :: Int
pointPhase = Int
0 }
integStartPoint :: Run -> Point
integStartPoint :: Run -> Point
integStartPoint Run
r = Int -> Point
point Int
nl
where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
(Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
point :: Int -> Point
point Int
n = Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
pointRun :: Run
pointRun = Run
r,
pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
pointPriority :: Int
pointPriority = Int
0,
pointIteration :: Int
pointIteration = Int
n,
pointPhase :: Int
pointPhase = Int
0 }
integStopPoint :: Run -> Point
integStopPoint :: Run -> Point
integStopPoint Run
r = Int -> Point
point Int
nu
where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
(Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
point :: Int -> Point
point Int
n = Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
pointRun :: Run
pointRun = Run
r,
pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
pointPriority :: Int
pointPriority = Int
0,
pointIteration :: Int
pointIteration = Int
n,
pointPhase :: Int
pointPhase = Int
0 }
simulationStopPoint :: Run -> Point
simulationStopPoint :: Run -> Point
simulationStopPoint Run
r = Run -> Double -> Int -> Point
pointAt Run
r (Specs -> Double
spcStopTime (Specs -> Double) -> Specs -> Double
forall a b. (a -> b) -> a -> b
$ Run -> Specs
runSpecs Run
r) Int
0
pointAt :: Run -> Double -> EventPriority -> Point
pointAt :: Run -> Double -> Int -> Point
pointAt Run
r Double
t Int
priority = Point
p
where sc :: Specs
sc = Run -> Specs
runSpecs Run
r
t0 :: Double
t0 = Specs -> Double
spcStartTime Specs
sc
dt :: Double
dt = Specs -> Double
spcDT Specs
sc
n :: Int
n = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor ((Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
dt)
p :: Point
p = Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
pointRun :: Run
pointRun = Run
r,
pointTime :: Double
pointTime = Double
t,
pointPriority :: Int
pointPriority = Int
priority,
pointIteration :: Int
pointIteration = Int
n,
pointPhase :: Int
pointPhase = -Int
1 }
integPointsStartingFrom :: Point -> [Point]
integPointsStartingFrom :: Point -> [Point]
integPointsStartingFrom Point
p = [Point]
points
where r :: Run
r = Point -> Run
pointRun Point
p
sc :: Specs
sc = Run -> Specs
runSpecs Run
r
(Int
nl, Int
nu) = Specs -> (Int, Int)
integIterationBnds Specs
sc
n0 :: Int
n0 = if Point -> Int
pointPhase Point
p Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
then Point -> Int
pointIteration Point
p
else Point -> Int
pointIteration Point
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
points :: [Point]
points = (Int -> Point) -> [Int] -> [Point]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Point
point [Int
n0 .. Int
nu]
point :: Int -> Point
point Int
n = Point { pointSpecs :: Specs
pointSpecs = Specs
sc,
pointRun :: Run
pointRun = Run
r,
pointTime :: Double
pointTime = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n Int
0,
pointPriority :: Int
pointPriority = Int
0,
pointIteration :: Int
pointIteration = Int
n,
pointPhase :: Int
pointPhase = Int
0 }
timeGrid :: Specs -> Int -> [(Int, Double)]
timeGrid :: Specs -> Int -> [(Int, Double)]
timeGrid Specs
sc Int
n =
let t0 :: Double
t0 = Specs -> Double
spcStartTime Specs
sc
t2 :: Double
t2 = Specs -> Double
spcStopTime Specs
sc
n' :: Int
n' = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Int
1
dt :: Double
dt = (Double
t2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n')
f :: Int -> (Int, Double)
f Int
i | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = (Int
i, Double
t0)
| Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
n' = (Int
i, Double
t2)
| Bool
otherwise = (Int
i, Double
t0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
dt)
in (Int -> (Int, Double)) -> [Int] -> [(Int, Double)]
forall a b. (a -> b) -> [a] -> [b]
map Int -> (Int, Double)
f [Int
0 .. Int
n']
delayPoint :: Point -> Int -> Point
delayPoint :: Point -> Int -> Point
delayPoint Point
p Int
dn
| Int
dn Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0 = String -> Point
forall a. HasCallStack => String -> a
error String
"Expected the positive number of iterations: delayPoint"
| Bool
otherwise =
let sc :: Specs
sc = Point -> Specs
pointSpecs Point
p
n :: Int
n = Point -> Int
pointIteration Point
p
ph :: Int
ph = Point -> Int
pointPhase Point
p
in if Int
ph Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0
then let t' :: Double
t' = Point -> Double
pointTime Point
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dn Double -> Double -> Double
forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc
n' :: Int
n' = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Integer) -> Double -> Integer
forall a b. (a -> b) -> a -> b
$ (Double
t' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Specs -> Double
spcStartTime Specs
sc) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
in Point
p { pointTime = t',
pointIteration = n',
pointPhase = -1 }
else let n' :: Int
n' = Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
dn
t' :: Double
t' = Specs -> Int -> Int -> Double
basicTime Specs
sc Int
n' Int
ph
in Point
p { pointTime = t',
pointIteration = n',
pointPhase = ph }