6.4 SNES: нелинейные решатели

Модуль SNES из библиотеки PETSC решает уравнения вида:

F(x) = 0.
(6.7)

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

xk+1 = x0 + [F ′(xk)]- 1F (xk)?k = 0,1...
(6.8)

где x0 – начальное решение, а F(xk) – Якобиан.

Метод Ньютона можно записать в виде следующих двух шагов:

1.
(приблизительно) решить F(xk)δxk = -F(xk)
2.
xk+1 = xk + δxk.

При создании решателя SNES вызываем следующую функцию:

    SNESCreate(MPI_Comm comm,SNES *snes);

Затем решаем систему:

    SNESSolve(SNES snes,Vec b,Vec x);

и вызываем

    SNESDestroy(SNES snes);

После создания решателя пользователь должен задать и определить функцию:

    SNESSetFunction(SNES snes, Vec f, PetscErrorCode (*FormFunction)(SNES snes, Vec x, Vec f, void *ctx), void *ctx);

Здесь ctx – некая структура, определенная пользователем, в которой хранятся параметры модели.

Также для решения нужно задавать Якобиан, его можно задать как вручную, подобно заданию функции formFunction, как и параметром, при указании которого решатель сам численно вычисляет якобиан, используя конечные разности вида:

F′≈ F-(u+-h-*dxi)--F(u).
 i           h
(6.9)

DMMG – высокоуровневая поддержка мультигрида. PETSC предоставляет простой высокоуровневый интерфейс для многосеточного метода. Рассмотрим небольшой пример с использованием DMMG:

    DMMG *dmmg; 
    DA da; 
    Vec soln;

Создаем DA, которая хранит информацию о сетке

    ierr = DACreate3d(PETSC_COMM_WORLD, DA_NONPERIODIC, DA_STENCIL_STAR, 
                      3, 3, 3, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 
                      1, 1, 0, 0, 0, da); CHKERRQ(ierr);

Создаем структуры данных DMMG

    ierr = DMMGCreate(PETSC_COMM_WORLD, 3, PETSC_NULL, dmmg); CHKERRQ(ierr); 
    ierr = DMMGSetDM(dmmg, (DM)da); CHKERRQ(ierr);

Задаем SNES решателю пользовательских функции, характеризующих решаемую задачу

    ierr = DMMGSetSNES(dmmg, myFunction, myJacobian); CHKERRQ(ierr);

Решаем задачу

    ierr = DMMGSolve(dmmg); CHKERRQ(ierr);

Получаем доступ к решению

    soln = DMMGGetx(dmmg); 
    ierr = DMMGDestroy(dmmg); CHKERRQ(ierr); 
    ierr = DADestroy(da); CHKERRQ(ierr);

Рассмотрим задаваемую пользователем функцию. В случае одной неизвестной в каждом узле сетки функция задается в следующем виде:

    int localfunction(DALocalInfo *info,double ***x, double ***f,void *ctx)

Если создаем двумерную задачу с тремя неизвестными в одном узле сетки (u,v,p), то для этого мы сначала должны создать структуру для их хранения:

typedef struct{ 
    PetscScalar u,v,p; 
} Field;

а затем и функцию, использующую это структуру

int localfunction(DALocalInfo *info, Field **x, Field **f, void *ctx)