PIC(Левченко)

Материал из theor
Перейти к: навигация, поиск

Программный комплекс CFhall

CFhall - это программный комплекс для моделирования плазмы методом "частица-в-ячейке". CFhall создан при подготовке к защите магистерской диссертации Перепёлкиной А.Ю., научный руководитель: Левченко В.Д.

Вопросы и комментарии о работе программы можно отправить на адрес .

Установка и запуск

В начале следует при помощи 'ssh' зайти на сервер bladerunner (194.67.66.81) (85.143.113.31). После этого Вы, скорее всего, окажетесь в своей домашней директории. Для того, чтобы получить версию программы, можно попробовать набрать

> bzr checkout /mpihome/NASTYA/CFHall_MEPhI

В появившейся директории CF_mephi нужно создать директорию для сброса результатов. По умолчанию 'Drop'

> cd CFHall_MEPhI
> mkdir Drop

Затем можно попробовать запустить расчет с теми параметрами, которые заданы по умолчанию

> srun ./runHall.py

Убедившись, что работает, можно остановить расчет, нажав 2 раза "crtl+c", и перейти к настройке параметров.

Если что-либо не получается, лучше отправить вопросы на указанный выше email.

Так как программа продолжает развиваться, чтобы убедиться в том, что у Вас последняя версия, можно набрать

> bzr update

!UPDATE! Для запуска из аудитории в ИПМ

Для начала работы:

  • Зайти как пользователь MEPhI_plasma
  • Открыть терминал (Konsole)
  • Создать свою директорию
> mkdir myname
  • перейти в эту директорию
> cd myname
  • получить программу
> git clone /mnt/Photon/Save/GIT-for-all/nastya/CFHall.git
  • перейти в директорию программы
> cd CFHall
  • Перейти к стабильной версии программы
> git checkout c9e79b2f2
  • создать директорию для вывода результатов
> ln -s /Run/MEPhI/Drop
  • попробовать запустить
> ./runHall.py

Если запуск успешен, начнется расчет. Его можно остановить, нажав Ctrl+c. Признаки того, что идет расчет:

  • в терминале появляются строчки, начинающиеся с #Step ..
  • в директории Drop появляются результаты.

Можно посмотреть результаты, и при помощи интерфейса задать свои параметры для расчета.

Интерфейс

При установке по вышеописанной процедуре Вы получаете его исходный код в виде набора файлов, содержащих исходный код программы.

Для того, чтобы использовать предоставленный код в своих целях, можно адаптировать код под свои задачи. Но, так как это, вероятно, довольно сложно, некоторое количество параметров, отвечающих постановке задачи, вынесены в файл 'runHall.py'. Этот файл и является "интерфейсом" программного комплекса CFhall.

Чтобы начать его править, откройте его в любом графическом редакторе, например, так

> kde-open runHall.py
  • Если вы открываете папку с кодом через проводник, не пытайтесь открыть 'runHall.py' двойным щелчком мыши. Программа запустится без вывода в консоль.

Синтаксис файла отвечает языку Python 2.7. Для понимания интерфейса знание этого языка не обязательно, но полезно. Основные необходимые сведения:

  • Регистр учитывается (tPulse ≠ TPulse, причем в коде C тоже)
  • Количество отступов от начала строки имеет значение. При неправильном количестве отступов может возникать ошибка
  • Изучение Python несложно. Достаточно прочитать 3-ю и 4-ю главы на этой странице [[1]]

Начало файла содержит импорт библиотек и настройку необходимых параметров:

#!/usr/bin/python
# -*- coding: utf-8 -*
from math import *
import os,sys
from k import *
mod = {}

Эту часть пользователю изменять не нужно.

Далее заведены некоторые переменные, которые можно и не импользовать, а рассматривать их как пример гибкости языка Python.

Далее приведена настройка состава плазмы. Задаются возможные сорта частиц:

sorts = {'e':(0,'(-1.0,1.0)'),
'p':(1,'(1.0, 1836.0)'),
'e_track':(2,'(-1.0,1.0)'),
'p_stop':(10,'stop(1.0)'),
'e_stop':(9,'stop(-1.0)')}
for k,v in sorts.iteritems(): mod[k]='zzz' #эта строчка нужна только для корректной работы интерфейса

Здесь через запятую перечислены формы вида:

'<название частицы>':(<уникальный номер от 0 до 10>, '(<заряд>, <масса>)')

Можно задавать неподвижные частицы. Им необходимо задать только заряд. Например, неподвижный протон задается так:

'p_stop':(10,'stop(1.0)')

Особенные свойства некоторых сортов:

  • Импульсы и положения частиц с уникальным номером сорта 0 на каждом Drop записываются в файл el***.pts. Все остальные - в файл io***.pts.
  • Частицы с сортом 2 отслеживаются: вдоль всей траектории на каждом шаге по времени в файл pts.trj записываются их координаты. Не стоит присваивать этот сорт более чем одной частице на расчет - файл становится плохо читаемым и может вырасти до слишком большого размера.

Затем задается количество частиц в каждой ячейке области, занимаемой плазмой. Например, так можно задать по одному протону электрону и протону в каждой ячейке:

PIC = ['e']+['p']

А так 5 электронов, 4 протона и один неподвижный протон:

PIC = 5*['e']+4*['p']+['p_stop']

Затем задаются начальные значения полей, импульсов и смещений частиц

E = (0.1,0,'pars.a*0.004')
B = (0,'2*sin(x)',0)
r = (-0.2,0,0)
p = (0,0,0)

Через запятую перечислены x, y и z компоненты. В приведенном примере показано, как можно задавать постоянное значение поля по области равное числу (здесь Ex = 0.1), равное некоторой переменной (Ez = 0.004*pars.a, что значит " pars.* " - см. список доступных переменных ниже). Можно задавать значения, зависящие от координат ячейки (x,y, а для z надо писать izs*dz) и ее номера вдоль соответствующей оси (ix,iy, izs)

Такой интерфейс подставляет эти числа в присваивания вида 'Ex = ... ' в коде программы. В более сложных случаях требуется вместо подстановки какого-либо значения в правую часть присваивания подставить более сложную функцию для начальных распределений. Тогда в начале подставляемой строки должен стоять символ обратной кавычки ' ` '. Вот пример для иллюстрации, где поле имеет вид, похожий на ступеньку:

s1 = '` if (ix == 0) {F.Ex.sRef(izs) = - pars.delta*(ka*pars.delta+1)*0.5*pars.dx; }\n\
else if (ix==64) {F.Ex.sRef(izs) = - pars.delta*0.5*(1.0-pars.delta)*pars.dx;}\n\
else if (ix>64) {F.Ex.sRef(izs) = - pars.delta * pars.dx;}\n\
else {F.Ex.sRef(izs) = 0.0;}'
E = (s1,0,0)

Здесь для удобной записи переменной Python 's1' присваивается некоторая строка, которая подставляется для вычисления Ex.

В строке после символа ' ` ' стоит код на языке С для вычисления Ex. При помощи этой возможности Вы фактически получаете возможность вставлять свои выражения (на С) в исполняемый код (объявлять свои переменные, производить вычисления, выводить на экран значения и т.п.).

В выражениях можно использовать глобальные переменные " pars.* ", а также координаты (x,y, izs*dz) и номера ячеек (ix,iy, izs)

Значения полей в коде записываются по подобию 'F.Ex.sRef(izs)' слева от знака присваивания ( = ) и 'F.Ex.sVal(izs)' справа от него. (вместо Ex может стоять любая компонента поля)

Символы '\n' нужны лишь для более аккуратной записи строки в код средствами Python, а символ '\' в конце строчки необходим, чтобы Python правильно обработал перенос строки.

Для задания распределения импульсов частиц по Максвеллу подготовлена переменная 'maxwell'. Такое распределение задается строкой

p = maxwell

Если в составе указаны частицы разных сортов, можно придавать им различные импульсы и смещения. По умолчанию всем будут присвоены значения из строк p=(..,..,..) r=(..,..,..). Чтобы, например, придать только электронам некоторое смещение, можно написать

mod['e'] = '{rx} = 0.1;'

Слева в квадратных скобках стоит название частицы. Справа может стоять произвольный код на C. В фигурных скобках могут стоять также величины rx,ry,rz,px,py,pz. В выражениях можно использовать глобальные переменные " pars.* ", а также координаты (x,y, izs*dz) и номера ячеек (ix,iy, izs). Вот еще пример:

mod['p'] = 'double kappa=0.2; {px} = kappa; {ry}=0.1;'

Далее задается форма электромагнитного импульса, заходящего в область с правой границы по оси x. По умолчанию импульс содержится в компонентах Bz и Ey. Его форма на этой границе представима в виде

(APulse/oPulse)*ShapeR (y,z)*ShapeT(t)*sin(oPulse * tPulse(y,z))

Это выражение передаюется в компоненту B_z на границе. Здесь APulse и oPulse - численные параметры - они задаются ниже. Импульс будет дейcтвовать в течение времени, определяемого значением параметра TPulse (Также задается среди численных параметров).

tPulse - это 'время' импульса. По умолчанию оно задается tPulse(y,z)=0 и течет, на каждой итерации расчета становясь больше на dt (шаг по времени). Если в начале поставить значение, не равное нулю, а некоторой зависимости от координат, можно задать пространственную зависимость фазы импульса. ( sin(oPulse*(t-phi(y,z)) ). Это можно использовать, например, для описания фокусирующихся импульсов.

Форма задается выражением вида:

  Pulse = {'R':' ShapeR = 1;',
           'T':'ShapeT=1;'
           'ON':0}

Здесь 'R' и 'T' - ключи, после которых можно написать строку, представляющую собой код на C, в результате которого из параметров y и z (для 'R') и из t (для 'T') присваивается значение переменным ShapeR и ShapeT.

В поле 'ON' надо написать 0, если импульс не нужен, 'x', чтобы импульc пошел с правой границы по оси x, 'y', чтобы импульc пошел с правой границы по оси y.

Более сложный пример:

Pulse = {'R':
'''    double R = 0.25*pars.b;
      double r2=y*y+z*z;
      double arg=M_PI*r2/(4.*R*R);
      ShapeR = (arg>1)?0:exp(-arg);,
     'T':'ShapeT =sin(tPulse*(M_PI/pars.TPulse));',
     'ON':'x'}}

Тройные кавычки использованы для того, чтобы Python видел выражение на нескольких строчках как одну строку. В этой строке объявляются некоторые переменные (R,r2,arg) чтобы более наглядно записать через них форму импульса. Внутри выражения под ключом 'R' при необходимости можно настроить и 'tPulse'.

В выражениях для формы можно использовать глобальные переменные " pars.* ", а также координаты (y,z) и время (tPulse).

mark - позволяет задать точку, в которой значения полей будут записываться в каждый момент времени. Если такая точка не нужна, нужно задать этот параметр пустым

 mark = ('','','',[''],['']) 

А так задается точка с координатами (17,16,Nz/2+1), в которой будут записываться значения компонент 'Ex','Ey', 'Ez', 'Bx','By', 'Bz', 'divE','divB','rho'.

mark = ('17','16','Nz/2+1',['Ex','Ey', 'Ez', 'Bx','By', 'Bz', 'divE','divB','rho'],
['rx','ry','rz', 'px', 'py', 'pz'])

Плазма по умолчанию занимает плоский слой перпендикулярно оси x. Границы слоя задаются ниже (см. xI, xX). Возможна более точная настройка. Если такая настройка не нужна, нужно оставить

 plasmaLayer = '' 

Если нужно задать блок, ограниченный по осям x и y, необходимо в эту строку вписать код на С, задающий границы блока, то есть параметры aL, aR, bL, bR. Тогда плазма будет занимать область aL<x<aR, bL<y<bR.

plasmaLayer = ' double aL=pars.xI*pars.a, aR=(pars.NxBrick-pars.xX)*pars.a, bL = pars.b*1./16., bR = pars.b*15./16.; '

Настроенные вышеописанным образом начальные и граничные распределения записываются в код программы следующей строкой:

gen_r0inithpp(sorts,PIC,mod,E,B,r,p,Pulse,mark,plasmaLayer)

Затем задаются ключевые параметры

Nz=120
MaxRank = 7
Xbrick = 2

Nz — количество ячеек по оси z (Nz). MaxRank определяет размеры по осям x и y. Размеры области составляет (2MaxRank * Xbrick) ячеек по оси x, 2MaxRank ячеек по оси y, Nz ячеек по оси z. Так, Xbrick - это количество 'brick' — кусков размером (2MaxRank x 2MaxRank x Nz) вдоль оси x.

Boundary определяет граничные условия при помощи четырех символов. m означает идеальный проводник, на границе с которым тангенциальные компоненты электрического поля равны нулю. Символом 'o' обазначение условия симметрии — на этой границе нормальные компоненты электрического поля равны нулю. Такое условие также аналогично существованию такой же области по ту сторону границы (поэтому - симметрия). В строке из четырех символов каждый отвечает определенной границе в следующем порядке (граница по x ближняя к началу координат (н.к.) - граница по y дальняя от н.к. - граница по х дальняя от н.к. - граница по y ближняя к н.к.)

Следующие 2 строки менять не нужно

genPICcode(Boundary, Pulse['ON'], mark = mark[2:], DivCalc=False) 
os.system("make NZ="+`Nz`+" MR="+`MaxRank`) 
  • Когда Вы изменяете MaxRank или Nz лучше пересобрать код заново. Для этого можно набрать 'make clean' в командной строке.

После этого компилируется код, и оставшаяся часть файла отвечает за настройку числовых параметров и запуск расчета. Эта команда настраивает шаги:

pars.setSteps(dt=0.004,dx=0.05,dy=0.05,dz=0.05)

Здесь шаги можно записывать либо напрямую числами, либо подговить заранее переменные dt, dx,dy,dz в Python. Возможны варианты записи

dx=dy=dz=0.05
dt=0.004
pars.setSteps(dt,dx,dy,dz)

а также

pars.setSteps(0.004,0.05,0.05,0.05)

Выставлены значения по умолчанию (равные dx=dy=dz=0.05, dt=0.004), поэтому возможен вариант

pars.setSteps(dt=0.03,dx=1)

в котором задаются только dt и dx, а остальные остаются равными значениям по умолчанию.

При выборе шагов необходимо соблюдать ограничения:

  • Условие Куранта
  • В длину волны проходящего импульса должно укладываться не менее 20 шагов сетки
  • В период колебаний любого изучаемого процесса должно укладываться не менее 20 шагов по времени
  • За время 4*dt частицы не должны преодолевать расстояние, большее dx. (Ограничение программной реализации — см. ниже описание численной схемы)

Далее задается директория, в которую сбрасываются результаты расчетов.

pars.DropDir = "Drop"

Директория должна быть создана заранее. Другие примеры: "Drop/var1", "Drop/var2"... Так можно запустить одновременно пару расчетов с разными параметрами с выводом результатов в различные директории.

Далее настраиваются параметры области

 pars.setPlasma(Nbr=XBrick, xX=0.9, xI=0.)

Nbr - количество 'brick' размером (2MaxRank x 2MaxRank x Nz) ячеек вдоль оси x. В итоге размер расчетной области составляет (2MaxRank·Nbr x 2MaxRank x Nz) ячеек.

Следующие параметры необходимо задавать в случае, если плазма занимает не всю расчетную область. xI, xX - размер области, занимаемой вакуумом слева и справа от плазмы вдоль оси x соответственно. Этот размер задается в долях от 'brick'. Чтобы задать область из четырех 'brick' так, чтобы только 2-й содержал плазму, можно написать так:

 pars.setPlasma(Nbr=4, xX=2, xI=1)

Значения по умолчанию Nbr=1, xX=0, xI=0.

Можно задать внешние поля. Эти поля влияют только на ускорение частиц и не участвуют в эволюции поля согласно уравнениям Максвелла. В отличие от самосогласованной части (задается выше начальными условиями) они не изменяются с течением времени.

pars.setExtFields(E0x = 0,E0y = 0,E0z = 0,B0x = 0,B0y = 0,B0z = 0)

В строке

pars.setv(vT=dx,sigma=dx,delta=0.125)

задаются различные (v - от слова various) параметры. Среди них sigma определяет температуру плазмы (импульсы задаются с распределением ~exp(p2/sigma2) ).

Остальные два параметра можно использовать на свое усмотрение, например, обращаясь к ним через 'pars.vT' и 'pars.delta' в написании граничных условий. (" pars.* " - см. список доступных переменных ниже)

Функция в следующей строке настраивает численные параметры импульсов.

pars.setPulses(o=2*pi, A=40.0/Focus**0.5,  T=200)

o - частота (oPulse), А - амплитуда (APulse), T - длительность (TPulse).

Далее следует обратить внимание на конструкцию

for step in xrange(20):
 if step:
   hall.runStep(step)
   hall.dropDia(step)
 if (step>5 and step%1 ==0): hall.dropFlds(step, "Jz,Ez,Ey,By,Bx")
 if (step%10 ==0): hall.dropPts(step)
 #if (step==16): hall.write_pulse();

Совершив 12*2MaxRank шагов по времени, программа делает 'Drop' - дает возможность закончить расчет и/или вывести результаты в указанную выше директорию.

В приведенном примере расчет будет идти, пока не будет совершено 20 'Drop'-ов. Каждый первый (step%1) шаг после пятого будут выводиться значения частиц полей указанные в строке

"Jz,Ez,Ey,By,Bx"

и каждый 10-й шаг будут выводиться в отдельные файлы данные о частицах.

Так как эти данные занимают много места на диске, иногда лучше уменьшить количество выводимых данных. Для этого можно выводить их реже (step%2, step%5 - каждый второй или пятый шаг), и можно выводить не все компоненты, удалив их из указанной выше строки.

Глобальные переменные

В выражениях для начальных и граничных условий может быть удобно использовать глобальные переменные кода

  • pars.dt - шаг по времени, задается в pars.setSteps(..) в интерфейсе.
  • pars.dx, pars.dy, pars.dz - шаги по координатам. Задаются в pars.setSteps(..)
  • pars.a, pars.b - размеры области по x и y. Равны 2MaxRank dx и 2MaxRank dy
  • pars.E0x, ... - внешние поля. Задаются в pars.setExtFields(..) в интерфейсе.
  • pars.oPulse, pars.APulse, pars.Tpulse - частота, амплитуда и длительность импульса. Задаются в pars.setPulses(..) в интерфейсе.
  • MaxRank
  • Nz

Константа π=3.1415926... записывается "M_PI" в коде C, и "pi" в коде Python.

Просмотр результатов

В указанной в интерфейсе директории (по умолчанию - 'Drop' ) по мере счета появляются файлы вида 'dropEx00000.lrc'. Их можно просматривать программой arr3D, набрав в командной строке, например, если директория для сброса результатов называется 'Drop':

> arr3D Drop/dropBy*

Об использовании arr3D можно узнать, набрав

> arr3D ?

В аудитории ИПМ им. М.В.Келдыша установлен еще одна программа для просмотра таких файлов. Она находится в директории /home/MEPhI/plasma/im3D Для запуска из директории CFHall можно запустить ее так:

> im3D Drop/dropEy*

После запуска можно нажать 'h', тогда в терминал выведется справка по этой программе. Например,

«ESC»  Выход из программы 
<Enter¦BackSpace>       Переход к следующему¦предыдущему массиву
a    Установка пределов палитры из пределов массива fMin..fMax

Чтобы картинка была больше, нужно поставить опцию --zoom:

> im3D --zoom "3 3 3" Drop/dropEy*

Диагностика средних значений

В папке Drop самостоятельно начинает формироваться файл avedia.dat.

Этот файл можно открыть почитать текстовым редактором, а можно строить графики при помощи, например gnuplot. На компьютерах в ИПМ им. М.В. Келдыша можно использовать gplt.

Столбцы файла

  • t --- время в безразмерных единицах
  • Np --- количество частиц
  • Nex --- сколько частиц передвинулось через границы ячеек

Эти величины усредняются по ячейкам:

  • E2x, E2y, E2z, B2x, B2y, B2z --- квадраты полей.
  • divE, divB --- дивергенции полей. Вероятно, не считаются, если отдельно не настроить.
  • Emx, Emy, Emz, Bmx, Bmy, Bmz --- максимальные по области значения модулей поля.

Следующие величины усредняются по частицам:

  • P2x P2y P2z --- импульсы.
  • Е --- энергия. Корень из квадрата импульса плюс один.
  • Wkx, Wky, Wkz --- приближение для кинетической энергии по компонентам. Квадрат соответствующего импульса, деленный на релятивистскую массу.

Содержание файлов avediaS#.dat, аналогично. Суммирование происходит только по одному сорту частиц. # означает номер сорта частиц.

Сведения о CFhall

Размер расчетной области составляет (2MaxRank·Nbr x 2MaxRank x Nz) ячеек.

Система единиц

Использована безразмерная система единиц, удобная для численного моделирования.

  • Масса измеряется в массах электрона me;
  • заряд — в зарядах электрона e,
  • скорости в скоростях света c,
  • время в обратных плазменных частотах электронов ωp-1,
  • единица для полей meωpc/e,
  • единица для размеров c/ωp.

Смещения частиц r (rx,ry,rz) измеряются относительно центра ячейки в долях от размера ячейки.

Численная схема

Для эволюции полей используетс схема FDTD. Поэтому при задании начальных условий лучше помнить, что значений полей на сетке задаются в соответствии с ячейкой Yee. Электрическое поле задается на гранях ячеек, магнитное — на ребрах. На рисунке ниже — одна ячейка размером dx*dy*dz.

Yee.png

Кроме этого, значения магнитного и элекрического полей разнесены на полшага по времени.

Для продвижения частиц и подсчета токов использован метод "частица-в-ячейке". Для учета разномасштабности задачи обновление токов происходит каждые 4 шага dt.