Планшет — не роскошь

«Вы назовете это извращением. Но кто сказал, что извращение — это плохо?»
Один доцент нашей кафедры

Идея описанного ниже эксперимента возникла после серии услышанных и подслушанных высказываний о том, что современные планшеты не могут выступать в роли инструментов серьезной научно-исследовательской деятельности. Действительно, для многих пользователей работа с планшетником сводится к веб-серфингу, переписке по электронной почте и разного протокола мессенджерах, чтению книг, просмотру видео и иным преимущественно развлекательным целям. Также, как следует ожидать и как показывает недавний пост, планшетник является прекрасным мобильным подспорьем при работе с «офисными» приложениями. Однако аппаратные характеристики существующей техники позволяют задуматься — а насколько же эффективным станет планшет на изначальном для компьютеров поприще.

Техника

Подопытным стал TeXeT TM-7025, с 1 ГГц процессором, 512 Мб оперативной памяти и Android 4.0 на борту. Мой первый компьютер, приобретенный в 2004 году, был чуточку мощнее (Athlon XP 1.66 ГГц, 256 Мб памяти до и 1 Гб — после апгрейда), однако заменен был только в марте 2012, намолотив Фортраном за свою жизнь расчётов на курсовую, а затем в паре с ноутбуком HP 550 — на два диплома, обеспечив успешное окончание школы, бакалавриата и магистратуры по направлению теоретической физики со специализацией на вычислительной гидродинамике.

Стоит сказать, что подопытный планшет имеет 7-дюймовый экран, что резко ограничивает удобство редактирования исходного кода, однако внешний компьютер для этих целей применялся по минимуму. Сторонними средствами написана уже данная статья, хотя, чего греха таить, введение и этот раздел набраны на смартфоне.

Софт

В силу рода занятий, возиться с портированием доступного софта не собирался даже в предельном случае разбора задачи. Потому, если бы такого не нашлось, то и эксперимента бы не было. Однако, обзор Google Play показал, что в нём наличествует порт Octave на Android. Помимо него, конечно же, имеется и весьма обширный спектр более простых программ, но с их вычислительными возможностями ясности нет, зато Octave — программа известная и небольшой опыт возни с ней уже имеется. Потому выбирать было не из чего, но возможность заодно разобраться и с этой системой более детально стала дополнительным стимулом в работе. Для рисования графиков имеющийся пакет требует установки графического пакета droidplot — портированной тем же автором версии gnuplot.

Естественно, код программы проще писать отдельным текстом и просто загружать его в консоли Octave. В качестве редактора был взял первый подвернувшийся под руку, то бишь бесплатная версия DroidEdit.

Задача

Пробные запуски были сделаны на системах дифференциальных уравнений. Просто загоняем систему и рисуем график на выходе, всё встроенными функциями Octave. Неинтересно, быстро. Но работает, что воодушевляет.

Тестовой стала небезызвестная система Лоренца, берущая своё начало, как известно, из уравнений конвекции в приближении Буссинеска.

function xdot=f(x,t)  r=28.; s=10.; b=8./3.;  # система уравнений Лоренца xdot(1)=s*(x(2)-x(1)); xdot(2)=r*x(1)-x(2)-x(1)*x(3); xdot(3)=x(1)*x(2)-b*x(3);  endfunction  # временные границы, начальные условия t=linspace(0.,10.,250); x0=[7.;10.;5.]; lsode_options("integration method","non-stiff"); y=lsode("f",x0,t);  # график аттрактора в плоскости x-z plot(y(:,1),y(:,3)); 

В зависимости от управляющего параметра (нормированного числа Рэлея), эта система проходит ряд устойчивых состояний от полного равновесия в нуле до хаотического режима. Сгенерированные на планшете картинки наглядно их демонстрируют. При r = 0.5:

r = 4:

r = 16:

r = 25:

r = 28:

Это всё замечательно, но какая вычислительная гидродинамика без конечно-разностных методов? Следующим этапом доказательства того, что планшет может не только потреблять контент, стало моделирование совершенно стандартной задачи конвекции жидкости в замкнутой полости. В рамках одного из курсов её обязан реализовать каждый студент наших кафедр. Формулировка следующая: рассчитать поле скорости и температуры несжимаемой вязкой жидкости в подогреваемой сбоку квадратной полости. По уму, конечно, разумно бы сперва посмотреть чего попроще, например, задачи теплопроводности, но решено было не мелочиться.

Уравнения сразу даны в безразмерном виде, сокращающем всё многообразие физических параметров задачи до двух — числа Грасгофа Gr и числа Прандтля Pr. На всех границах скорость, и значит — функция тока — равна нулю (условие прилипания для вязкой жидкости), температура слева — 0, справа — 1, на верхней и нижней границе линейно растёт.

Для упрощения, задача рассматривается как двумерная — полость представляется бесконечно протяженным каналом квадратного сечения. Трёхмерные задачи — это гораздо более широкое богатство содержания, но в то же время уже совершенно иной уровень по требованиям к ресурсам, потому их численное моделирование широко развиваться стало только в последние пару-тройку десятилетий. Двумерный же подход позволяет выявить в первую очередь основные, всеобщие свойства гидродинамических систем, с гораздо меньшими усилиями. Прежде всего, связано это с особенностями двумерного течения жидкости — вместо векторных уравнений Навье-Стокса можно перейти к описанию скорости через две скалярные величины (т.н. двухполевой метод) — функцию тока (она же z-компонента векторного потенциала скорости) и завихренность (z-компонента ротора скорости). Примечательны они тем, что, во-первых, при переходе к функции тока автоматически и точно удовлетворяется условие несжимаемости (что и является основной проблемой моделирования в естественных переменных скорость-давление), а, во-вторых, с помощью функции тока крайне удобно изображать поле скорости — её изолинии соответствуют линиям тока, главное, разобраться с направлением течения.

Уравнения двухполевого метода:

Программа

Программа незамысловата, и являет собой простейшую явную конечно-разностную схему. Потому останавливаться детально на ней незачем. Просто приведу исходник. Даже не особо оптимизированный, хотя с учётом производительности планшета вопросы оптимальности кода становятся крайне существенными.

a = 8.0; ap = 8.0; maxp = 100;  Nx = 10; Ny = 10;  tmax = 10.;  Gr = 20000.; Pr = 1.; eps = 1.0e-2;  # выходные файлы o1id = fopen("/sdcard/Octave/Convection/psi(t)_Gr=20000.txt","w"); o2id = fopen("/sdcard/Octave/Convection/field_Gr=20000.txt","w");  # массивы/матрицы # завихренность phi0 = zeros(Nx+1,Ny+1,"single"); phi1 = zeros(Nx+1,Ny+1,"single");  # функция тока psi0 = zeros(Nx+1,Ny+1,"single"); psi1 = zeros(Nx+1,Ny+1,"single");  # температура T0 = zeros(Nx+1,Ny+1,"single"); T1 = zeros(Nx+1,Ny+1,"single");  # вспомогательно-оптимизационные параметры hx = 1.0 / Nx; hy = 1.0 / Ny; hxi = 1.0 / hx; hyi = 1.0 / hy;  ht = hy**2/ a; htp = hy**2 / ap;  htx2 = ht*hxi**2; hty2 = ht*hyi**2; htxy = 0.25*ht*hxi*hyi;  htpx2 = htp*hxi**2; htpy2 = htp*hyi**2; htpxy = 0.25*htp*hxi*hyi;  Pri = 1.0 / Pr;  # метки осей x = linspace(0.,1.,Nx+1); y = linspace(0.,1.,Ny+1); axis("xy");  # начальное распределение for i = 1:Nx+1  for j = 1:Ny+1   psi0(i,j) = 1.0e-1*(1.0d0 - (i-1)*hx)*(i-1)*hx*(1.0d0 - (j-1)*hy)*(j-1)*hy;   T0(i,j) = (i-1)*hx;   phi0(i,j) = 0.0;  endfor endfor  ct = 0.; q = 0;  # цикл по времени while(ct <= tmax)  # уравнения для phi, T  for i = 2:Nx   for j = 2:Ny    dpsidx = psi0(i+1,j) - psi0(i-1,j);    dpsidy = psi0(i,j+1) - psi0(i,j-1);     phi1(i,j) = phi0(i,j) + ...    (phi0(i+1,j) - 2.*phi0(i,j) + phi0(i-1,j))*htx2 + ...     (phi0(i,j+1) - 2.*phi0(i,j) + phi0(i,j-1))*hty2 + ...     htxy*( dpsidx*(phi0(i,j+1) - phi0(i,j-1)) - dpsidy*(phi0(i+1,j) - phi0(i-1,j)) ) + ...     0.5*ht*hxi*Gr*(T0(i+1,j) - T0(i-1,j));     T1(i,j) = T0(i,j) + Pri*( (T0(i+1,j) - 2.*T0(i,j) + T0(i-1,j))*htx2 + ...     (T0(i,j+1) - 2.*T0(i,j) + T0(i,j-1))*hty2 ) + ...     htxy*( dpsidx*(T0(i,j+1) - T0(i,j-1)) - dpsidy*(T0(i+1,j) - T0(i-1,j)) );   endfor  endfor    # граничные условия  for i = 1:Nx+1   T1(i,1) = (i-1)*hx;   T1(i,Ny+1) = (i-1)*hx;    phi1(i,1) = -2.*Nx*Nx*psi1(i,2);   phi1(i,Ny+1) = -2.*Nx*Nx*psi1(i,Ny);  endfor   for j = 1:Ny+1   T1(1,j) = 0.;   T1(Nx+1,j) = 1.;    phi1(1,j) = -2.*Ny*Ny*psi1(2,j);   phi1(Nx+1,j) = -2.*Ny*Ny*psi1(Nx,j);  endfor  # уравнение Пуассона для функции тока  p = 0;  ppp0 = 0;  ppp1 = 0;    do     for i = 1:Nx+1    for j = 1:Ny+1     psi1(i,j) = 0.;    endfor   endfor     for i = 2:Nx    for j = 2:Ny     psi1(i,j) = psi0(i,j) + htp*phi1(i,j) + ...      (psi0(i+1,j) - 2.*psi0(i,j) + psi0(i-1,j))*htpx2 + ...      (psi0(i,j+1) - 2.*psi0(i,j) + psi0(i,j-1))*htpy2;       ppp0=ppp0 + abs(psi0(i,j));       ppp1=ppp1 + abs(psi1(i,j));    endfor   endfor    for i = 1:Nx+1    psi1(i,1) = 0.;    psi1(i,Ny+1) = 0.;   endfor 		   for j = 1:Ny+1    psi1(1,j) = 0.;    psi1(Ny+1,j) = 0.;   endfor    for i = 1:Nx+1    for j = 1:Ny+1     psi0(i,j) = psi1(i,j);    endfor   endfor    p = p++;  until(p > maxp || abs(ppp1 - ppp0)/(ppp0+ppp1) < eps)    for i = 1:Nx+1    for j = 1:Ny+1     phi0(i,j) = phi1(i,j);     psi0(i,j) = psi1(i,j);     T0(i,j) = T1(i,j);    endfor   endfor  # максимальное значение функции тока   if(q == 1000)    fprintf(o1id,"%f %f\n",ct,max(max(psi0)));    q = 0;   endif    ct = ct + ht   q++; endwhile  # поля переменных  for i = 1:Nx+1  for j = 1:Ny+1   fprintf(o2id,"%f %f %f %f %f\n",(i-1)*hx,(j-1)*hy,psi0(i,j),T0(i,j),phi0(i,j));  endfor endfor  fclose(o1id); fclose(o2id);  # рисование и сохранение графиков # в силу невыясненных причин потребовалось транспонировать выводимые матрицы # как видно, генерируются картинки в eps-формате # похоже, порт Octave пока попросту других делать не умеет # можно строить в gnuplot / droidplot figure(1); contourf(x,y,T0'); title("T, Gr = 2000"); xlabel("x"); ylabel("y"); saveas(1,"/sdcard/Octave/Convection/T_20000.eps");  figure(1); contourf(x,y,psi0'); title("Psi, Gr = 2000"); xlabel("x"); ylabel("y"); saveas(1,"/sdcard/Octave/Convection/psi_20000.eps"); 

Результаты

Ну, и самое интересное. Скорость и температура жидкости при различной интенсивности подогрева, то бишь — разных Gr. Решения получены на сетке размером 10х10 узлов. Конечно, это тоже немного, однако, когда компьютеры были большими, нормой бывали расчёты и на сетках 4х4, 5х5.

Со значениями основных параметров Gr ~ 10 000, Pr = 1, eps = 0.01 и кодом, приведённым выше, планшет обрабатывает одну единицу безразмерного времени в течение примерно трёх-пяти минут, что вовсе не так уж и медленно, с учётом недостатка памяти, необходимости Android’у выполнять иные системные функции и отвоёвывать на то ресурсы системы, а также и на потенциальную медлительность самой системы Octave как интерпретатора.

При больших Gr расчёт двухвихревого режима течения занимает гораздо больше времени из-за существенного ухудшения сходимости итерационного процесса. Вызвано это ухудшение тем, что двухвихревое течение является не стационарным, а установившимся колебательным процессом — вихри то усиливаются, то ослабевают, попутно изменяя свои размеры в довольно широком диапазоне. Например, Gr = 120 000 на планшете считается примерно часа полтора-два, но и на обычной машине процесс идёт не очень радостно по сравнению со слабым подогревом. Да и улететь в NaN на таком режиме уже легче лёгкого, требуется значительно уменьшать шаг времени. Ускорение алгоритма в данной ситуации — вопрос особый, и требует пересмотра метода решения разностных уравнений.

Итак, полученные картинки. Просто и наглядно — видно, что в полости существует устойчивое вихревое течение и выраженный перенос тепла в первую очередь за счёт движения жидкости — всплывания у нагретой границы и погружения у холодной. При повышении интенсивности подогрева течение сосредотачивается в пограничном слое и переходит в двухвихревой режим. Цветовая шкала не рисовалась, но она стандартная — рост от фиолетового к коричневому через зелёный (как на географических картах).

Функция тока и температура при Gr = 2000:

при Gr = 20000:

и при Gr = 120000:

Заключение

Планшет — действительно полноценная вычислительная машина. Слабенькая, но упорная. И производить на нём несложные вычисления, скажем, чтобы проверить свежую идею или прикинуть основные направления дальнейшего поиска, абсолютно реалистично.

Не стоит забывать, что любой компьютер всегда остаётся компьютером в изначальном смысле этого слова.