1
+ """
2
+ module TestLinearEquations - Test LinearEquations(..) in ModiaBase/src/EquationAndStateInfo.jl
3
+ """
4
+ module TestLinearEquations
5
+
6
+ using Test
7
+ using Modia
8
+ using Modia. TimerOutputs
9
+
10
+ println (" ... Test TestLinearEquations" )
11
+
12
+ # Solve A*x = b as residual equation
13
+ const A = [1.0 2.0 3.0 4.0 5.0 ;
14
+ 0.0 2.0 0.0 0.0 0.0 ;
15
+ 0.0 0.0 3.0 0.0 0.0 ;
16
+ 0.0 0.0 0.0 4.0 0.0 ;
17
+ 0.0 0.0 0.0 0.0 5.0 ]
18
+ const b = [- 1.0 , - 2.0 , - 3.0 , - 4.0 , - 5.0 ]
19
+ const isInitial = true # the first call must be performed with isInitial=true
20
+ const time = 0.0 # only for warning/error messages
21
+ const timer = TimerOutput () # timer to measure the time to solve the linear equation system
22
+ const leq = Modia. LinearEquations ([" x2" , " x1" ], [:x1 ], [1 ,5 ], 1 , false )
23
+
24
+ computeResidual1 (x1) = A* x1- b
25
+ computeResidual2 (x1,x2) = x1[1 ] + 2 * x1[2 ] + 3 * x1[3 ] + 4 * x2
26
+
27
+ # Solution
28
+ const x1_sol = A\ b
29
+ const x2_sol = - (x1_sol[1 ] + 2 * x1_sol[2 ] + 3 * x1_sol[3 ])/ 4.0
30
+
31
+ # Test function for ODE mode
32
+ function solveLinearEquation (leq,isInitial,time,timer)
33
+ local x1:: Vector{Float64} , x2
34
+ leq. mode = - 3
35
+ while Modia. LinearEquationsIteration! (leq, isInitial, time, timer)
36
+ x1 = leq. x_vec[1 ]
37
+ x2 = leq. x[1 ]
38
+ append! (leq. residuals, computeResidual1 (x1))
39
+ append! (leq. residuals, computeResidual2 (x1,x2))
40
+ end
41
+ return (x1,x2)
42
+ end
43
+
44
+ # Test function for DAE mode
45
+ function solveLinearEquation (leq,isInitial,solve_leq,isStoreResult,time,timer)
46
+ local x1, x2
47
+ leq. mode = - 3
48
+ while Modia. LinearEquationsIteration! (leq, isInitial, solve_leq, isStoreResult, time, timer)
49
+ x1 = leq. x_vec[1 ]
50
+ x2 = leq. x[1 ]
51
+ append! (leq. residuals, computeResidual1 (x1))
52
+ append! (leq. residuals, computeResidual2 (x1,x2))
53
+ end
54
+ return (x1,x2)
55
+ end
56
+
57
+ @testset " Test TestLinearEquations" begin
58
+ # Test ODE mode (x is always explicitly computed by solving the linear equation system)
59
+ leq. odeMode = true
60
+ (x1,x2) = solveLinearEquation (leq,isInitial,time,timer)
61
+
62
+ d1 = x1_sol- x1
63
+ err1 = sqrt (d1' * d1) + abs (x2_sol - x2)
64
+ @test isapprox (err1, 0.0 , atol= 1e-13 )
65
+
66
+
67
+ # Test DAE mode
68
+ leq. odeMode = false
69
+
70
+ # x is computed explicitly at initialization or at an event
71
+ (x1,x2) = solveLinearEquation (leq,false ,true ,false ,time,timer)
72
+ d2 = x1_sol- x1
73
+ err2 = sqrt (d2' * d2) + abs (x2_sol - x2)
74
+ @test isapprox (err2, 0.0 , atol= 1e-13 )
75
+
76
+ # x is computed by DAE solver
77
+ leq. x = [6.0 , 1.0 , 2.0 , 3.0 , 4.0 , 5.0 ] # Provide x from DAE solver
78
+ (x1,x2) = solveLinearEquation (leq,false ,false ,false ,time,timer)
79
+ residuals = leq. residuals # Store residuals in DAE solver
80
+ d3 = leq. x - vcat (x2,x1)
81
+ d4 = vcat (A* leq. x[2 : 6 ] - b, leq. x[2 ] + 2 * leq. x[3 ] + 3 * leq. x[4 ] + 4 * leq. x[1 ]) - residuals
82
+ err3 = sqrt (d3' * d3)
83
+ err4 = sqrt (d4' * d4)
84
+ @test isapprox (err3, 0.0 , atol= 1e-15 )
85
+ @test isapprox (err4, 0.0 , atol= 1e-15 )
86
+
87
+ # x is computed by DAE solver at a communication point (leq.residuals is not provided)
88
+ leq. x = [6.0 , 1.0 , 2.0 , 3.0 , 4.0 , 5.0 ] # Provide x from DAE solver
89
+ (x1,x2) = solveLinearEquation (leq,false ,false ,true ,time,timer)
90
+ d5 = leq. x - vcat (x2,x1)
91
+ err5 = sqrt (d5' * d5)
92
+ @test isapprox (err5, 0.0 , atol= 1e-15 )
93
+ end
94
+
95
+ end
0 commit comments