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
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: Method -> Method -> Bool
$c/= :: Method -> Method -> Bool
== :: Method -> Method -> Bool
$c== :: Method -> Method -> Bool
Eq, Eq 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
min :: Method -> Method -> Method
$cmin :: Method -> Method -> Method
max :: Method -> Method -> Method
$cmax :: Method -> Method -> Method
>= :: Method -> Method -> Bool
$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
compare :: Method -> Method -> Ordering
$ccompare :: Method -> Method -> Ordering
Ord, Int -> Method -> ShowS
[Method] -> ShowS
Method -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Method] -> ShowS
$cshowList :: [Method] -> ShowS
show :: Method -> String
$cshow :: Method -> String
showsPrec :: Int -> Method -> ShowS
$cshowsPrec :: Int -> 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 <- forall a. a -> IO (IORef a)
newIORef Bool
False
IORef Double
t <- forall a. a -> IO (IORef a)
newIORef forall a b. (a -> b) -> a -> b
$ Specs -> Double
spcStartTime Specs
specs
PriorityQueue (Point -> IO ())
pq <- forall a. IO (PriorityQueue a)
PQ.newQueue
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 = forall a b. (RealFrac a, Integral b) => a -> b
round ((Specs -> Double
spcStopTime Specs
sc forall a. Num a => a -> a -> a
-
Specs -> Double
spcStartTime Specs
sc) forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc)
in if Int
n forall a. Ord a => a -> a -> Bool
< Int
0
then
forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$
String
"The iteration number in the stop time has a negative value. " forall a. [a] -> [a] -> [a]
++
String
"Either the simulation specs are incorrect, " forall a. [a] -> [a] -> [a]
++
String
"or a floating point overflow occurred, " forall a. [a] -> [a] -> [a]
++
String
"for example, when using a too small integration time step. " forall a. [a] -> [a] -> [a]
++
String
"You have to define this time step regardless of " forall a. [a] -> [a] -> [a]
++
String
"whether you actually use it or not, " forall a. [a] -> [a] -> [a]
++
String
"for Aivika allows combining the ordinary differential equations " forall a. [a] -> [a] -> [a]
++
String
"with the discrete event simulation within one model. " forall a. [a] -> [a] -> [a]
++
String
"So, if you are still using the 32-bit architecture and " forall a. [a] -> [a] -> [a]
++
String
"you do need a small integration time step " forall a. [a] -> [a] -> [a]
++
String
"for integrating the equations " forall a. [a] -> [a] -> [a]
++
String
"then you might think of using the 64-bit architecture. " forall a. [a] -> [a] -> [a]
++
String
"Although you could probably just forget " forall a. [a] -> [a] -> [a]
++
String
"to increase the time step " 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 forall a. Ord a => a -> a -> Bool
< Int
0 then
forall a. HasCallStack => String -> a
error String
"Incorrect phase: basicTime"
else
Specs -> Double
spcStartTime Specs
sc forall a. Num a => a -> a -> a
+ Double
n' forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc forall a. Num a => a -> a -> a
+ forall {a}. (Eq a, Num a) => Method -> a -> Double
delta (Specs -> Method
spcMethod Specs
sc) Int
ph
where n' :: Double
n' = 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 forall a. Fractional a => a -> a -> a
/ Double
2
delta Method
RungeKutta4 a
2 = Specs -> Double
spcDT Specs
sc 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 forall a. Fractional a => a -> a -> a
/ Double
3
delta Method
RungeKutta4b a
2 = Double
2 forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc 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 = 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 = 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 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 = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor ((Double
t forall a. Num a => a -> a -> a
- Double
t0) 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 forall a. Eq a => a -> a -> Bool
== Int
0
then Point -> Int
pointIteration Point
p
else Point -> Int
pointIteration Point
p forall a. Num a => a -> a -> a
+ Int
1
points :: [Point]
points = 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' = forall a. Ord a => a -> a -> a
max (Int
n forall a. Num a => a -> a -> a
- Int
1) Int
1
dt :: Double
dt = (Double
t2 forall a. Num a => a -> a -> a
- Double
t0) forall a. Fractional a => a -> a -> a
/ (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n')
f :: Int -> (Int, Double)
f Int
i | Int
i forall a. Eq a => a -> a -> Bool
== Int
0 = (Int
i, Double
t0)
| Int
i forall a. Eq a => a -> a -> Bool
== Int
n' = (Int
i, Double
t2)
| Bool
otherwise = (Int
i, Double
t0 forall a. Num a => a -> a -> a
+ (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i) forall a. Num a => a -> a -> a
* Double
dt)
in 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 forall a. Ord a => a -> a -> Bool
<= Int
0 = 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 forall a. Ord a => a -> a -> Bool
< Int
0
then let t' :: Double
t' = Point -> Double
pointTime Point
p forall a. Num a => a -> a -> a
- forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
dn forall a. Num a => a -> a -> a
* Specs -> Double
spcDT Specs
sc
n' :: Int
n' = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a b. (RealFrac a, Integral b) => a -> b
floor forall a b. (a -> b) -> a -> b
$ (Double
t' forall a. Num a => a -> a -> a
- Specs -> Double
spcStartTime Specs
sc) forall a. Fractional a => a -> a -> a
/ Specs -> Double
spcDT Specs
sc
in Point
p { pointTime :: Double
pointTime = Double
t',
pointIteration :: Int
pointIteration = Int
n',
pointPhase :: Int
pointPhase = -Int
1 }
else let n' :: Int
n' = Int
n 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 :: Double
pointTime = Double
t',
pointIteration :: Int
pointIteration = Int
n',
pointPhase :: Int
pointPhase = Int
ph }