Mid-term evaluation has now arrived so it's time to summarize the work I've done and check which goals I have achieved. During this period I enjoied working with the community and the advices given by the mentors and by the other members have been really helpful. The most important change to the project is that, discussing with mentors, we decided to start implementing ode15i because it's more general than ode15s and to build ode15s later around it.
Here you can find the code I've written so far:
https://bitbucket.org/Francesco_Faccio/octave
The most difficult task in the first part of this project was to have Octave compiled with link to Sundials. After having accomplished this, I checked the presence and usability of the library nvector_serial which contains the implementation of IDADENSE and IDABAND modules. I aggregated its build flags with the flags of sundials_ida and included the header nvector_serial.h in the dld-function.
I checked the licenses of Sundials and SUPERLUMT (a package which will be used as a sparse direct solver, independent from Sundials) and they have 3-Clause license, so they are compatible with GNU-license and can be used.
After configuring the correct flags, I began writing a minimal wrapper of ode15i of the form:
[t, y] = ode15i (odefun, tspan, y0, yp0, options)
The first problem was to deal with Sundials' own types. Sundials uses realtype and N_Vector. An N_Vector is a vector of realtype, while a realtype can be both a float, a double or a long double, depending on how Sundials has been built (default type is double). I assumed to use the default double realtype and wrote functions N_Vector ColToNVec (ColumnVector data, long int n) and ColumnVector NVecToCol (N_Vector v, long int n), which convert an Octave ColumnVector to an N_Vector and viceversa.
I checked some minimal input conditions, wrote a few input validation tests, set AbsTol, Reltol, tspan, y0 and yp0 of type realtype or N_Vector.
Once preprocessed data, the moment of glory of Sundials' functions arrived. The first call was to IDACreate(), which creates an IDA memory block and returns a pointer which is then passed as the first argument to all subsequent ida function calls.
Sundials then asks to provide a function which computes the residual function in the DAE. This function must have the form:
int flag = resfun (realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user data)
As a temporary workaround I wrote a function which converts yy, yp and tt in ColumnVector(s), uses feval to evaluate the DAE (passed through a global pointer of type octave_function) and put the output in rr.
Then a call to IDAInit, IDASVtolerances (or IDASStolerances if AbsTol is a scalar), IDADense and IDADlsSetDenseJacFn (if supplied) sets up the linear solver.
Sundials accepts only a Jacobian Function of the form J = DF/DY + cj*DF/DYP, where cj is a scalar proportional to the inverse of the step size (cj is computed by IDA's solver). I used the same workaround of the residual to evaluate J when it's required.
Finally in the main loop a call to IDASolve solves the DAE system and gives the solution in output.
What this wrapper can solve:
this wrapper of ode15i can solve a system of differential equations of the form f(t, y, y') = 0, integrating from t0 to tf with initial conditions y0 and yp0. The output of the function is the solution of the DAE evaluated ONLY in the points supplied in tspan.
It accepts as option a scalar RelTol and a scalar or vector AbsTol. If the user wants to supply a Jacobian, it must be of the form J = DF/DY + cj DF/DYP. Both the system of DAE and the Jacobian must be a function handle.
I tested this wrapper using the 2 benchmark problems described in the previous post.
In Robertson chemical kinetics problem I found the right solution both passing a Jacobian and letting Sundials approximate it. The script I used is the following:
function res = robertsidae(t, y, yp)
res = [-(yp(1) + 0.04*y(1) - 1e4*y(2)*y(3));
-(yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2);
y(1) + y(2) + y(3) - 1];
endfunction
function jacc = jacobian(t, y, yp, c)
jacc = [-0.04-c, 1e4*y(3), 1e4*y(2);
0.04, -c-6*1e7*y(2)-1e4*y(3), -1e4*y(2);
1, 1, 1];
endfunction
options = odeset('RelTol', 1e-3, 'AbsTol', 1e-6, ...
'Jacobian', @jacobian);
y0 = [1; 0; 0];
yp0 = [-1e-4; 1e-4; 0];
tspan = [0 4*logspace(-6, 6)];
[t, y] = ode15i(@robertsidae, tspan, y0, yp0, options);
y(:, 2) = 1e4*y(:, 2);
figure(1);
semilogx(t, y)
xlabel('time');
ylabel('species concentration');
title('Robertson DAE problem with a Conservation Law');
legend ('y1', 'y2', 'y3');
The second problem was to find the solution of a 2-D heat equation semidiscretized to a DAE on the unit square, as described in the previous post.
After discretizing the domain, the system of 100 differential algebraic equations was passed to ode15i (as sparse methods are in progress, I used the dense solver of Sundials also for this problem, without passing the Jacobian).
Here you can find the script I wrote to solve the problem. A 3-D time dependent plot shows how the solution evolves in time and space.
uu0 = zeros(100, 1);
%Initialize uu in all grid points
for j = 1:10
yfact = (1 / 9).*(j - 1);
offset = 10.*(j - 1);
for i = 1:10
xfact= (1 / 9).*(i - 1);
loc = offset + (i - 1);
uu0(loc + 1) = 16 * xfact * (1 - xfact) * yfact * (1 - yfact);
endfor
endfor
up0 = zeros(100, 1);
%Set values of uu and up at boundary points
for j = 1:10
offset = 10 * (j - 1);
for i = 1:10
loc = offset + (i - 1);
if (j == 1 || j == 10 || i == 1 || i == 10 )
uu0 (loc + 1) = 0; up0 (loc + 1) = 0;
endif
endfor
endfor
function res = klu (t, uu, up)
res = [ uu(1);
uu(2);
uu(3);
uu(4);
uu(5);
uu(6);
uu(7);
uu(8);
uu(9);
uu(10);
uu(11);
up(12) - 81 * (uu(11) + uu(13) + uu(2) + uu(22) - 4*uu(12));
up(13) - 81 * (uu(12) + uu(14) + uu(3) + uu(23) - 4*uu(13));
up(14) - 81 * (uu(13) + uu(15) + uu(4) + uu(24) - 4*uu(14));
up(15) - 81 * (uu(14) + uu(16) + uu(5) + uu(25) - 4*uu(15));
up(16) - 81 * (uu(15) + uu(17) + uu(6) + uu(26) - 4*uu(16));
up(17) - 81 * (uu(16) + uu(18) + uu(7) + uu(27) - 4*uu(17));
up(18) - 81 * (uu(17) + uu(19) + uu(8) + uu(28) - 4*uu(18));
up(19) - 81 * (uu(18) + uu(20) + uu(9) + uu(29) - 4*uu(19));
uu(20);
uu(21);
up(22) - 81 * (uu(21) + uu(23) + uu(12) + uu(32) - 4*uu(22));
up(23) - 81 * (uu(22) + uu(24) + uu(13) + uu(33) - 4*uu(23));
up(24) - 81 * (uu(23) + uu(25) + uu(14) + uu(34) - 4*uu(24));
up(25) - 81 * (uu(24) + uu(26) + uu(15) + uu(35) - 4*uu(25));
up(26) - 81 * (uu(25) + uu(27) + uu(16) + uu(36) - 4*uu(26));
up(27) - 81 * (uu(26) + uu(28) + uu(17) + uu(37) - 4*uu(27));
up(28) - 81 * (uu(27) + uu(29) + uu(18) + uu(38) - 4*uu(28));
up(29) - 81 * (uu(28) + uu(30) + uu(19) + uu(39) - 4*uu(29));
uu(30);
uu(31);
up(32) - 81 * (uu(31) + uu(33) + uu(22) + uu(42) - 4*uu(32));
up(33) - 81 * (uu(32) + uu(34) + uu(23) + uu(43) - 4*uu(33));
up(34) - 81 * (uu(33) + uu(35) + uu(24) + uu(44) - 4*uu(34));
up(35) - 81 * (uu(34) + uu(36) + uu(25) + uu(45) - 4*uu(35));
up(36) - 81 * (uu(35) + uu(37) + uu(26) + uu(46) - 4*uu(36));
up(37) - 81 * (uu(36) + uu(38) + uu(27) + uu(47) - 4*uu(37));
up(38) - 81 * (uu(37) + uu(39) + uu(28) + uu(48) - 4*uu(38));
up(39) - 81 * (uu(38) + uu(40) + uu(29) + uu(49) - 4*uu(39));
uu(40);
uu(41);
up(42) - 81 * (uu(41) + uu(43) + uu(32) + uu(52) - 4*uu(42));
up(43) - 81 * (uu(42) + uu(44) + uu(33) + uu(53) - 4*uu(43));
up(44) - 81 * (uu(43) + uu(45) + uu(34) + uu(54) - 4*uu(44));
up(45) - 81 * (uu(44) + uu(46) + uu(35) + uu(55) - 4*uu(45));
up(46) - 81 * (uu(45) + uu(47) + uu(36) + uu(56) - 4*uu(46));
up(47) - 81 * (uu(46) + uu(48) + uu(37) + uu(57) - 4*uu(47));
up(48) - 81 * (uu(47) + uu(49) + uu(38) + uu(58) - 4*uu(48));
up(49) - 81 * (uu(48) + uu(50) + uu(39) + uu(59) - 4*uu(49));
uu(50);
uu(51);
up(52) - 81 * (uu(51) + uu(53) + uu(42) + uu(62) - 4*uu(52));
up(53) - 81 * (uu(52) + uu(54) + uu(43) + uu(63) - 4*uu(53));
up(54) - 81 * (uu(53) + uu(55) + uu(44) + uu(64) - 4*uu(54));
up(55) - 81 * (uu(54) + uu(56) + uu(45) + uu(65) - 4*uu(55));
up(56) - 81 * (uu(55) + uu(57) + uu(46) + uu(66) - 4*uu(56));
up(57) - 81 * (uu(56) + uu(58) + uu(47) + uu(67) - 4*uu(57));
up(58) - 81 * (uu(57) + uu(59) + uu(48) + uu(68) - 4*uu(58));
up(59) - 81 * (uu(58) + uu(50) + uu(49) + uu(69) - 4*uu(59));
uu(60);
uu(61);
up(62) - 81 * (uu(61) + uu(63) + uu(52) + uu(72) - 4*uu(62));
up(63) - 81 * (uu(62) + uu(64) + uu(53) + uu(73) - 4*uu(63));
up(64) - 81 * (uu(63) + uu(65) + uu(54) + uu(74) - 4*uu(64));
up(65) - 81 * (uu(64) + uu(66) + uu(55) + uu(75) - 4*uu(65));
up(66) - 81 * (uu(65) + uu(67) + uu(56) + uu(76) - 4*uu(66));
up(67) - 81 * (uu(66) + uu(68) + uu(57) + uu(77) - 4*uu(67));
up(68) - 81 * (uu(67) + uu(69) + uu(58) + uu(78) - 4*uu(68));
up(69) - 81 * (uu(68) + uu(60) + uu(59) + uu(79) - 4*uu(69));
uu(70);
uu(71);
up(72) - 81 * (uu(71) + uu(73) + uu(62) + uu(82) - 4*uu(72));
up(73) - 81 * (uu(72) + uu(74) + uu(63) + uu(83) - 4*uu(73));
up(74) - 81 * (uu(73) + uu(75) + uu(64) + uu(84) - 4*uu(74));
up(75) - 81 * (uu(74) + uu(76) + uu(65) + uu(85) - 4*uu(75));
up(76) - 81 * (uu(75) + uu(77) + uu(66) + uu(86) - 4*uu(76));
up(77) - 81 * (uu(76) + uu(78) + uu(67) + uu(87) - 4*uu(77));
up(78) - 81 * (uu(77) + uu(79) + uu(68) + uu(88) - 4*uu(78));
up(79) - 81 * (uu(78) + uu(70) + uu(69) + uu(89) - 4*uu(79));
uu(80);
uu(81);
up(82) - 81 * (uu(81) + uu(83) + uu(72) + uu(92) - 4*uu(82));
up(83) - 81 * (uu(82) + uu(84) + uu(73) + uu(93) - 4*uu(83));
up(84) - 81 * (uu(83) + uu(85) + uu(74) + uu(94) - 4*uu(84));
up(85) - 81 * (uu(84) + uu(86) + uu(75) + uu(95) - 4*uu(85));
up(86) - 81 * (uu(85) + uu(87) + uu(76) + uu(96) - 4*uu(86));
up(87) - 81 * (uu(86) + uu(88) + uu(77) + uu(97) - 4*uu(87));
up(88) - 81 * (uu(87) + uu(89) + uu(78) + uu(98) - 4*uu(88));
up(89) - 81 * (uu(88) + uu(90) + uu(79) + uu(99) - 4*uu(89));
uu(90);
uu(91);
uu(92);
uu(93);
uu(94);
uu(95);
uu(96);
uu(97);
uu(98);
uu(99);
uu(100)];
endfunction
tspan = linspace(0, 0.3, 100);
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, y] = ode15i(@klu, tspan, uu0, up0, options);
sol = zeros(10, 10, 100);
for z = 1:100
for i = 1:10
sol(:, i, z) = y(z, (((i - 1) * 10) + 1):(i * 10));
endfor
endfor
figure(1);
for k = 1:100
surf(sol(:, :, k))
axis([0 10 0 10 0 1]);
xlabel('X');
ylabel('Y');
zlabel('Heat');
title('2-D heat equation semidiscretized to a DAE on the unit square');
pause(0.03)
endfor
This is the solution at time t = 0.0091:
Larger problems will be used in order to test the efficiency of the code, because this 2 were solved almost immediately.
Nessun commento:
Posta un commento