Nieliniowe zagadnienie różniczkowe w PETSc
W tym rozdziale zajmiemy się (wciąż bardzo elementarnym)
przykładem nieliniowego równania pochodzącego z dyskretyzacji równania
konwekcji-dyfuzji
Dla prostoty, to równanie będziemy aproksymować
metodą różnic skończonych, dostając problem dyskretny postaci:
Problem.
Dla danego wektora znaleźć spełniający
dla ,
.
Jak widzimy, mamy tym razem do czynienia z układem
N równań nieliniowych, gdzie N oznacza liczbę wewnętrznych węzłów
siatki (warunki brzegowe Dirichleta wyznaczają od razu u0
i uN+1, więc możemy ograniczyć się jedynie
do obliczenia u w pozostałych węzłach). Jeśli więc siatka jest gęsta,
to rozmiar zadania nieliniowego jest bardzo duży, a macierz pochodnej rozrzedzona.
Właśnie to powoduje, że stosowanie PETSc do problemów pochodzących z dyskretyzacji
nieliniowych równań różniczkowych cząstkowych zaczyna być sensowne. Narzędzia
dostępne w PETSc (a więc także w SNESie) były tworzone właśnie z myślą
o efektywnych metodach rozwiązywania problemów powstałych z dyskretyzacji
równań różniczkowych cząstkowych.
Wstępna analiza problemu: wybór strategii
Zobaczymy (nareszcie!) jak, używając PETSc, napisać
zwarty, prosty, a do tego elegancki kod, rozwiązujący nieliniowe równanie
różniczkowe. Na naszym przykładzie nauczymy się pewnych technik, skutecznych
również w bardziej skomplikowanych przypadkach.
Zacznijmy od przeanalizowania problemu. To jest
bardzo ważne, albowiem od tego, w jaki sposób podejdziemy do zagadnienia,
będzie zależało wiele naszych następnych posunięć. Mamy do czynienia z
nieliniowym zagadnieniem dyskretnym, które zamierzamy rozwiązywać metodą
Newtona. Musimy więc najpierw przeformułować to zagadnienie jako zero pewnego
odwzorowania
. Oznaczając przez
uh wektor rozwiązania:
,
możemy zapisać zagadnienie dyskretne w postaci wektorowej:
,
gdzie Ah oznacza macierz (minus)
laplasjanu, a Bh macierz dyskretyzacji pierwszej pochodnej
różnicą centralną, tzn.:
Zwróćmy uwagę, że tym razem operację mnożenia wektorów
rozumiemy nie w sensie iloczynu skalarnego, tylko jako operację mnożenia
odpowiadających współrzędnych.
Przy takim zapisie łatwo już określić funkcję, której
zerem będzie rozwiązanie naszego zagadnienia:
Jej pochodną będzie naturalnie taki operator liniowy ,
który w działaniu na wektor w da wektor F'(u)w, określony
przez . Określiliśmy
więc oraz wyprowadziliśmy
wzór na F'. Jak pamiętamy, SNES będzie potrzebował procedur obliczających i.
Musimy teraz zastanowić się, jak zaimplementować te funkcje w PETSc. Tutaj
pojawia się następna zaleta zapisu wektorowego: daje się on bardzo łatwo
"przetłumaczyć" na kod źródłowy PETSc!
Definicja funkcji F
Popatrzmy, z jaką łatwością i elegancją możemy
zaimplementować procedurę obliczającą wartości funkcji
F. Elegancja
jest tu równie ważna, albowiem, wbrew twierdzeniom rozmaitych domorosłych
hackerów i innym pozorom twierdzimy, że czytelny i zrozumiały, choć być
może nieoptymalny kod jest bardziej wartościowy od pełnego sztuczek superniezrozumiałego
superkodu. Z naszego punktu widzenia, przy pisaniu dużej aplikacji gigantyczną
porcję czasu zajmuje odpluskwianie programu. Szczególnie niewdzięcznie
wygląda to zadanie w przypadku aplikacji numerycznych, takich, jakimi się
tu zajmujemy (cóż, przykre doświadczenia życiowe...) Elegancki, zwarty,
modułowy kod pozwala na lepsze śledzenie działania programu i szybszą lokalizację
błędów. Poza tym, należy zdawać sobie sprawę z tego, że dobry zapis wektorowy
algorytmów numerycznych ma istotne znaczenie, nie tylko dla przejrzystości
kodu, ale także - jego efektywności na nowoczesnych architekturach, zob.
[17], [5].
Efektywność programu to nie tylko szybkość jego działania, ale także
szybkość doprowadzenia go do działania. Z drugiej strony, nie należy
wpadać w skrajności i np. poświęcać całą szybkość programu dla piękna kodu
źródłowego...
Otóż przypuśćmy, że już zdefiniowaliśmy gdzieś w programie macierze
Ah i Bh dyskretyzacji laplasjanu i
pierwszej pochodnej (skądinąd doskonale wiemy, jak to zrobić!). Wtedy reszta
to fraszka: wystarczy wygenerować wektor prawej strony i w procedurze zapisać,
że wartość F(X) oblicza się w/g następującego algorytmu:
|
|
macierz laplasjanu pomnożona przez X funkcją MatMult() |
|
|
wynik mnożenia macierzy pierwszej pochodnej przez X,
razy (po współrzędnych) X - funkcje MatMult() oraz VecPointwiseMult(), |
|
|
prawa strona |
Co prawda widać, że w tym celu będzie nam trzeba
kilku wektorów pomocniczych, tak, abyśmy mieli gdzie przechowywać wyniki
pośrednie, ale prostota definicji F jest uderzająca. Części liniowe
F określają nam odpowiednie macierze, nieliniowości implementujemy
wykonując odpowiednie operacje na wektorach.
Czy z pochodną
F'(u) pójdzie nam równie
łatwo?...
Pochodna F
Jak wiadomo, pochodną odwzorowania
możemy reprezentować jako macierz
.
Zwróćmy jednak uwagę na pewien detal: było nam wygodniej zapisać pochodną
F nie jako macierz pochodnych cząstkowych, lecz jako operator, którego
wartości na zadanym wektorze
w są zdefiniowane pewnym wzorem. W
taki właśnie sposób definiuje się w PETSc macierze powłokowe (
shell
matrices), o których była mowa na stronie 18. Idea macierzy powłokowych
wyjątkowo dobrze sprawdzi się w naszym zagadnieniu:
-
Pochodną było nam po prostu łatwiej określić operatorowo, niż macierzowo.
Niby mała rzecz, bo i tak umiemy różniczkować, ale każdy, kto doświadczył
drobnej pomyłki w trakcie obliczeń na papierze, doceni fakt, że nie musimy
tego robić, o ile używamy macierzy powłokowej.
-
Przypomnijmy sobie, że w każdym kroku metody Newtona, macierz jakobianu
F'(X) zmienia się, ponieważ zmienia się punkt, w którym jest wyznaczana.
Gdybyśmy używali postaci macierzowej, za każdym razem musielibyśmy uaktualniać
elementy tej macierzy. Mając macierz powłokową, nie musimy robić nic, ponieważ
zmiany są uwzględniane na bieżąco.
-
Do rozwiązywania układów z macierzą pochodnej F'(X) w PETSc będziemy
zapewne stosować jakąś metodę iteracyjną, a te wymagają jedynie operacji
mnożenia przez macierz. Ponieważ i tak na początku procesu Newtona nie
mamy dobrego przybliżenia, będziemy mogli pozwolić sobie na mniej dokładne
rozwiązywanie układów z macierzą jakobianu w pierwszych kilku krokach,
co dramatycznie obniży koszt całego procesu. Mówimy o tym więcej w rozdziale
"Inne opcje".
-
Wreszcie, co też jest ważne, wykorzystamy już istniejące obiekty.
W definicji funkcji F używaliśmy macierzy laplasjanu i pochodnej,
i teraz te same macierze posłużą nam do zdefiniowania F'. Nie będziemy
w bólach definiować element po elemencie macierzy jakobianu!
Wobec powyższego, wystarczy odpowiednio zdefiniować
macierz powłokową.
Koncept kontekstu
Pomyślmy chwilkę o tym, jak w praktyce implementowalibyśmy
funkcję F. Jak wiadomo, SNES wymaga od nas funkcji formatu int
Funkcja(SNES MojSnes,Vec X,Vec F,void * Ctx). My zaś dla wyliczenia
F potrzebujemy odwołać się do mnożenia przez macierze Ah
i Bh. Jak to zrealizować w praktyce? Ktoś bardzo początkujący
mógłby pomyśleć o takiej realizacji (tak jak podczas rzeczywistego procesu
zastanawiania się nad implementacją, nie dbamy o szczegóły syntaktyczne
wywołań poszczególnych funkcji, skupiając się na ich istotnych argumentach):
int Funkcja( SNES MojSnes, Vec X, Vec F, void
* Ctx )
{
Mat A, B;
/* generuj A i B, a nastepnie przypisz
im wartosci */
MatCreate(A); MatCreate(B);
MatSetValues(A); ...itd....
/* uzyj A i B do obliczenia F */
MatMult(A); VecAXPY(); MatMult(B);
itd...
/* zwolnij macierze A i B */
MatDestroy(A); MatDestroy(B);
}
Zaproponowana funkcja w zasadzie nie wygląda źle:
macierze A i B są używane lokalnie przez F, więc
są też zadeklarowane lokalnie. Po wykorzystaniu są niszczone, bo przy następnym
wywołaniu funkcja zainicjuje inne zmienne wewnętrzne. Niestety, ta implementacja
jest nie do przyjęcia, bo jest kompletnie nieefektywna. Zdajmy sobie sprawę
z tego, że przy każdym obliczeniu wartości F(X), będą generowane
dwie (całkiem duże) macierze, wpisywane do nich wartości - i za każdym
razem dokładnie te same! Jest to szalona rozrzutność, dająca w konsekwencji
owszem, działający kod, ale potwornie wolny...
Podstawową wadą zaproponowanej procedury jest
wielokrotne robienie tego samego: generowanie wciąż tych samych macierzy.
Najprostszym sposobem usunięcia tego problemu byłoby jednokrotne utworzenie
i wygenerowanie elementów macierzy A i B (raz, a dobrze!)
gdzieś na zewnątrz procedury Funkcja(). Świetnie, ale jak
skorzystać z tych macierzy wewnątrz procedury? Znów, najszybciej
przychodzi nam do głowy rozwiązanie najprostsze, ale jednocześnie najbardziej
bałaganiarskie: po prostu, trzeba zadeklarować te macierze jako zmienne
globalne!
Mat A, B; /* nasze macierze jako zmienne globalne:
deklarowane na zewnatrz wszystkich funkcji */
/******************************************************
Program glowny
******************************************************/
int main(int argc, char **argv)
{
/* generuj A i B, a nastepnie przypisz
im wartosci */
MatCreate(A); MatCreate(B);
MatSetValues(A); ...itd...
/* ......... */
SNESSetFunction( Funkcja );
/* ......... */
/* zwolnij macierze A i B */
MatDestroy(A); MatDestroy(B);
return(0);
}
/******************************************************
Procedura obliczajaca F(X)
******************************************************/
int Funkcja( SNES MojSnes, Vec X, Vec F, void
* Ctx )
{
/* uzyj A i B do obliczenia F */
MatMult(A); VecAXPY(); MatMult(B);
itd...
return(0);
}
To już jest całkiem dobre, a zwłaszcza - znacznie bardziej efektywne,
niż poprzednia propozycja. Niemniej, kiedy wiemy, jak to zrobić efektywnie,
powinniśmy spróbować zrobić to elegancko i temu właśnie będzie służył tytułowy
koncept kontekstu użytkownika (pomysł na grę słów zaczerpnęliśmy
od B.Smitha). Cóż, eleganckim rozwiązaniem byłoby przekazanie do Funkcji()
macierzy Ah i Bh (uprzednio wygenerowanych
na zewnątrz procedury mnożenia!) w postaci parametrów. Z tym jednak
pozornie jest kłopot: format Funkcji() jest sztywny i nie możemy
dowolnie modyfikować listy jej argumentów. Na szczęście, ostatnim dopuszczalnym
argumentem funkcji jest wskaźnik do kontekstu użytkownika, zawierającego
pomocnicze dla procedury obiekty i właśnie w ten sposób przekażemy nasze
macierze!
Jak taki kontekst powinien wyglądać? Zasadniczo,
będzie to struktura w C, zawierająca parametry, które zechcemy przekazać
do Funkcji(). Tymi parametrami mogą być macierze, wektory robocze
(o tak, przydadzą się!), a nawet wskaźniki do funkcji udostępniających
zmienne lokalne funkcjom używającym naszego kontekstu. Taka struktura,
ze zmiennymi lokalnymi, funkcjami dostępu itp. mocno pachnie porządną klasą
w C++ i tak jest w rzeczywistości: programowanie z użyciem kontekstów ma
wiele cech programowania obiektowego. Osoby dociekliwe z pewnością zauważyły,
że typy (konteksty!) Vec, Mat,SLES itp. maja
bardzo wiele cech regularnych klas w C++.
Wracając do naszej konkretnej sytuacji, kontekst,
jaki chcielibyśmy przekazać do Funkcji(), powinien zawierać dwie
macierze: dyskretyzacji laplasjanu i pierwszej pochodnej. Kontekst ma być
strukturą, więc wygodnie nam będzie najpierw zdefiniować odpowiedni typ
strukturalny:
typedef struct {
Mat A;
Mat B;
} USERCTX;
a następnie zadeklarować zmienną tego typu, zaincjalizować
kontekst, wykorzystać i porzucić:
int main(int argc, char **argv)
{
USERCTX UserCtx;
/* inicjalizacja kontekstu */
USERCTXCreate( &UserCtx );
/* ....wykorzystanie kontekstu....
*/
SNESSetFunction( MojSnes, X, R, Funkcja,
(void *) &UserCtx );
/* zwolnienie kontekstu */
USERCTXDestroy( &UserCtx );
}
Funkcje USERCTXCreate() oraz USERCTXDestroy(),
także zdefiniowane przez nas, ładnie oddzielą od reszty kodu etap "tworzenia''
kontekstu (u nas: wygenerowanie macierzy i ich elementów) i jego usunięcia:
int USERCTXCreate( USERCTX *UserCtx )
{
MatCreate( &UserCtx->A ); MatCreate(
&UserCtx->B );
MatSetValues( UserCtx->A );
...itd....
}
int USERCTXDestroy( USERCTX *UserCtx )
{
MatDestroy( UserCtx->A ); ...itd...
}
Teraz użycie naszego kontekstu w Funkcji()
wyglądałoby mniej więcej tak:
int Funkcja( SNES MojSnes, Vec X, Vec F, void
* Ctx )
{
USERCTX *UserCtx;
/* najpierw rzutujemy wskaznik na
odpowiedni typ */
UserCtx = (USERCTX *) Ctx;
/* uzyj A i B do obliczenia F */
MatMult( UserCtx->A, X, F ); VecAXPY();
MatMult(); itd...
return(0);
}
Implementacja
Podobnie jak w przypadku układu 2x2, funkcję F
i pochodną F' zdefiniujemy w osobnych plikach. Kod źródłowy funkcji
main(), zapisany w kolejnym pliku, będzie zawierał zwarty kod
realizujący ogólny schemat użycia SNESu.
Zaczynamy od kontekstu użytkownika
W kilku miejscach (w pliku definiującym funkcję
F, w pliku definiującym pochodną F', w funkcji main())
będziemy się odwoływać do naszego kontekstu. Dobrze więc byłoby, gdybyśmy
dołączali definicję tego typu z jednego, wspólnego pliku nagłówkowego,
powiedzmy, user.h. Oto jego zawartość:
user.h
static Scalar zero = 0.0;
static Scalar one = 1.0;
static Scalar mone = -1.0;
typedef struct
{
Mat Lap; /* macierz laplasjanu
*/
Mat Poch; /* macierz pierwszej pochodnej */
Vec X; /* wektor rozwiazania
*/
Vec F; /*
prawej strony */
Vec R; /*
residuum
*/
Vec y; /*
pomocniczy
*/
Vec z; /*
pomocniczy
*/
} USERCTX;
Oprócz macierzy Ah i Bh,
które tutaj nazwaliśmy bardziej opisowo Lap i Poch, w
kontekście umieściliśmy jeszcze kilka wektorów roboczych, które wykorzystamy
potem w procedurach Funkcja() i Jakobian().
Dodatkowo, zdefiniowaliśmy parę stałych typu Scalar,
bardzo przydatnych do wykonania operacji obliczenia sumy lub różnicy dwóch
wektorów przy pomocy funkcji VecAXPY(). Będą one dostępne w każdym
pliku zawierającym dyrektywę #include "user.h".
Implementacja funkcjiF(X)
Procedura Funkcja() definiuje sposób
wyliczania wartości F(X). Oto plik zawierający jej kod źródłowy:
cdfun.c
#include "snes.h" /* plik z definicjami funkcji
i obiektow SNES; zawiera takze
dyrektywe #include "sles.h" */
#include "user.h" /* kontekst uzytkownika
i przydatne stale */
/******************************************************
Procedura obliczajaca wartosci funkcji F
******************************************************/
int Funkcja(SNES MojSnes,Vec X,Vec F,void
*Ctx)
{
int ierr;
USERCTX *UserCtx = (USERCTX *)Ctx;
/* Pomnoz przez pochodna centralna...
*/
ierr = MatMult( UserCtx->Poch, X, UserCtx->y
); CHKERRQ(ierr);
/* ...i domnoz przez u */
ierr = VecPointwiseMult( X, UserCtx->y,
F ); CHKERRQ(ierr);
/* Pomnoz przez laplasjan */
ierr = MatMult( UserCtx->Lap, X, UserCtx->y
); CHKERRQ(ierr);
/* Zsumuj! */
ierr = VecAXPY( &one, UserCtx->y,
F ); CHKERRQ(ierr);
ierr = VecAXPY( &mone, UserCtx->F,
F ); CHKERRQ(ierr);
return(ierr);
}
Wykorzystujemy wektory robocze kontekstu
UserCtx
do przechowywania wyników pośrednich, obliczanych według algorytmu zaproponowanego
w sekcji
"Definicja funkcji F" w
tym rozdziale. Zwróćmy ponadto uwagę na użycie funkcji
VecAXPY()
ze stałymi
one i
mone (definiowanymi w pliku
"user.h")
dla obliczenia sumy i różnicy dwóch wektorów.
Implementacja macierzy pochodnej F'(X)
Jak pamiętamy, zdecydowaliśmy się implementować
F'(X) w postaci macierzy powłokowej. Definiowanie macierzy powłokowej
odbywa się w dwóch etapach: najpierw tworzymy powłokę komendą
MatShellCreate(),
a następnie przekazujemy do niej funkcje, realizujące operacje dopuszczalne
na naszej macierzy, np. operację mnożenia:
MatShellSetOperation()
z odpowiednimi parametrami (patrz rozdział
"Macierze
powłokowe"). Zakładając, że powłokowa macierz jakobianu jest już dobrze
zdefiniowana - robimy to w programie głównym -możemy przystąpić do implementacji
funkcji
Jakobian(), przekazującą SNESowi macierz pochodnej: i
tu czeka nas bardzo miła niespodzianka!
cdjac.c
/******************************************************
Definicja pochodnej F
******************************************************/
#include "snes.h"
#include "user.h"
/* Implementacja funkcji Jakobian(), przekazywanej
do SNESu */
int Jakobian(SNES MojSnes, Vec X, Mat *Jac,
Mat *Precond, MatStructure *flag, void *Ctx)
{
*flag = SAME_NONZERO_PATTERN;
return(0);
}
/* Implementacja mnozenia przez macierz jakobianu
*/
int MultJakobian( Mat Jac, Vec VecIn, Vec
VecOut )
{
int ierr;
USERCTX *UserCtx;
/* Najpierw musimy wydobyc kontekst
*/
ierr = MatShellGetContext( Jac, (void
**)&UserCtx );CHKERRQ(ierr);
/* Pomnoz przez pochodna centralna...
*/
ierr = MatMult( UserCtx->Poch, VecIn,
UserCtx->y ); CHKERRQ(ierr);
/* ...i domnoz przez u */
ierr = VecPointwiseMult( UserCtx->X,
UserCtx->y, VecOut ); CHKERRQ(ierr);
/* Pomnoz przez laplasjan */
ierr = MatMult( UserCtx->Lap, VecIn,
UserCtx->y ); CHKERRQ(ierr);
/* Zsumuj! */
ierr = VecAXPY( &one, UserCtx->y,
VecOut ); CHKERRQ(ierr);
ierr = MatMult( UserCtx->Poch, UserCtx->X,
UserCtx->y ); CHKERRQ(ierr);
ierr = VecPointwiseMult( VecIn, UserCtx->y,
UserCtx->z ); CHKERRQ(ierr);
ierr = VecAXPY( &one, UserCtx->z,
VecOut ); CHKERRQ(ierr);
return(ierr);
}
Tak jest, funkcja Jakobian() właściwie
nie robi nic. Rzeczywiście, skoro macierz jakobianu zadana jest
przez operację mnożenia przez wektor - a to jest określone wcześniej przy
definiowaniu macierzy powłokowej - wyznaczanie macierzy jakobianu w specjalnej
procedurze jest zbędne. Niemniej, musimy przekazać jakąś sensowną
funkcję Jakobian() do SNESu, więc stąd nasza definicja.
Ponieważ wszystko załatwia operacja mnożenia,
jest to dobre miejsce na zdefiniowanie procedury mnożenia wektora przez
macierz jakobianu. Ta procedura zostanie potem przekazana jako argument
MatShellSetOperation(), dopełniając definiowanie powłokowej macierzy
F'(X). Procedurę MultJakobian(), definiującą de facto
macierz F'(X), wykorzystamy przy okazji tworzenia powłokowej macierzy
jakobianu w programie głównym.
Program główny: wykorzystanie SNESu
Główny plik, poza funkcją main(), zawiera także
funkcje tworzenia i niszczenia kontekstu użytkownika, USERCTXCreate()
i USERCTXDestroy().
cd.c
#include <stdio.h>
#include <math.h> /* zawiera
funkcje sin(x) */
#include "snes.h"
#include "user.h"
static char help[] = "Rozwiazuje nieliniowe
jednowymiarowe rownanie\n\
konwekcji-dyfuzji na odcinku [0,L]\n\
Opcje: -N [int] liczba wewnetrznych wezlow
siatki\n\
-L [Scalar]
dlugosc odcinka\n\n";
/* Deklaracje procedur obliczajacych funkcje
i jakobian */
extern int MultJakobian( Mat, Vec, Vec);
int Jakobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
int Funkcja(SNES,Vec,Vec,void*);
int USERCTXCreate(USERCTX*,int,Scalar);
int USERCTXDestroy(USERCTX*);
/**************************************************************
Program glowny
**************************************************************/
int main( int argc, char **argv )
{
SNES MojSnes;
/* obiekt SNES, ktorym rozwiazemy zadanie */
Mat
MatJ; /* macierz
jakobianu */
USERCTX UserCtx;
/* kontekst uzytkownika */
int
n;
/* liczba wezlow siatki */
static Scalar L = M_PI; /* dlugosc
przedzialu calkowania,
zatem krok siatki h = L/n */
int
ierr, flg, its;
PetscInitialize( &argc, &argv,PETSC_NULL,help
);
/* Interpretacja opcji uzytkownika
*/
ierr = OptionsGetInt( PETSC_NULL,
"-N", &n, &flg ); CHKERRA(ierr);
if( !flg ) n = 3;
ierr = OptionsGetScalar( PETSC_NULL,
"-L", &L, &flg );
CHKERRA(ierr);
/* Tworzenie potrzebnych obiektow:
SNESu, wektorow, macierzy pomocniczych */
ierr = SNESCreate(MPI_COMM_WORLD,SNES_NONLINEAR_EQUATIONS,&MojSnes);
CHKERRA(ierr);
ierr = USERCTXCreate( &UserCtx,
n, L ); CHKERRA(ierr);
VecView( UserCtx.F, VIEWER_DRAWX_WORLD
); PetscSleep(5);
/* Przekaz SNESowi funkcje i wektor
residualny */
ierr = SNESSetFunction(MojSnes, UserCtx.R,
Funkcja, &UserCtx); CHKERRA(ierr);
/* Przekaz SNESowi macierz jakobianu
i funkcje obliczajaca jakobian */
ierr = MatCreateShell(MPI_COMM_WORLD,n,n,n,n,&UserCtx,&MatJ);
CHKERRA(ierr);
ierr = MatShellSetOperation(MatJ,
MATOP_MULT, MultJakobian ); CHKERRA(ierr);
ierr = SNESSetJacobian(MojSnes,MatJ,MatJ,Jakobian,&UserCtx);
CHKERRA(ierr);
/* Uwzglednij opcje uruchomienia SNESu
typu -snes_view -snes_monitor */
ierr = SNESSetFromOptions(MojSnes);
CHKERRA(ierr);
/* Rozwiaz! */
ierr = SNESSolve(MojSnes,UserCtx.X,&its);
CHKERRA(ierr);
PetscPrintf(MPI_COMM_WORLD, "Liczba
iteracji = %d\n", its);
PetscPrintf(MPI_COMM_WORLD,"Wektor
rozwiazania:\n");
ierr = VecView(UserCtx.X, VIEWER_DRAWX_WORLD);
CHKERRA(ierr);
PetscSleep( 10 );
/* Zwolnij pamiec */
ierr = USERCTXDestroy(&UserCtx);
CHKERRA(ierr);
ierr = MatDestroy(MatJ); CHKERRA(ierr);
ierr = SNESDestroy(MojSnes); CHKERRA(ierr);
PetscFinalize();
return 0;
}
/**************************************************************
Funkcje zwiazane z kontekstem uzytkownika
**************************************************************/
int USERCTXCreate( USERCTX * UserCtx, int
n, Scalar L )
{
extern int GenerujLaplasjan1D( Mat, Scalar
);
extern int GenerujPochodnaCentralna1D( Mat,
Scalar );
int ierr=0;
int i;
double val;
double h = L/(n+1); /* krok siatki */
/* tworzymy wektory uzywane w kontekscie
naszego zagadnienia */
ierr = VecCreate(MPI_COMM_WORLD,PETSC_DECIDE,n,&UserCtx->X);
CHKERRQ(ierr);
ierr = VecDuplicate(UserCtx->X,&UserCtx->R);
CHKERRQ(ierr);
ierr = VecDuplicate(UserCtx->X,&UserCtx->F);
CHKERRQ(ierr);
ierr = VecDuplicate(UserCtx->X,&UserCtx->y);
CHKERRQ(ierr);
ierr = VecDuplicate(UserCtx->X,&UserCtx->z);
CHKERRQ(ierr);
/* tworzymy macierze uzywane w kontekscie
naszego zagadnienia */
ierr = MatCreate(MPI_COMM_WORLD,n,n,&UserCtx->Lap);
CHKERRQ(ierr);
ierr = MatCreate(MPI_COMM_WORLD,n,n,&UserCtx->Poch);
CHKERRQ(ierr);
/* generujemy elementy macierzy pomocniczych
*/
ierr = GenerujLaplasjan1D( UserCtx->Lap,
(Scalar) h ); CHKERRQ(ierr);
ierr = MatView( UserCtx->Lap, VIEWER_DRAWX_WORLD
); CHKERRQ(ierr);
PetscSleep( 3 );
ierr = GenerujPochodnaCentralna1D(
UserCtx->Poch, (Scalar) h );CHKERRQ(ierr);
ierr = MatView( UserCtx->Poch, VIEWER_DRAWX_WORLD
); CHKERRQ(ierr);
PetscSleep( 3 );
/* generujemy prawa strone */
for( i = 1 ; i <= n ; i++ )
{
val = sin(i*h) + sin(i*h)*cos(i*h); /* rozwiazanie u(x) = sin(x) */
VecSetValue( UserCtx->F, i-1, (Scalar)val, INSERT_VALUES );
}
/* Przyblizenie poczatkowe */
/* Zakladamy, ze zerowe... Od razu
wpisujemy do wektora rozwiazania X */
ierr = VecSet(&zero, UserCtx->X);
CHKERRQ(ierr);
return(ierr);
}
int USERCTXDestroy( USERCTX * UserCtx )
{
int ierr;
ierr = VecDestroy(UserCtx->X); CHKERRQ(ierr);
ierr = VecDestroy(UserCtx->R); CHKERRQ(ierr);
ierr = VecDestroy(UserCtx->F); CHKERRQ(ierr);
ierr = VecDestroy(UserCtx->z); CHKERRQ(ierr);
ierr = VecDestroy(UserCtx->y); CHKERRQ(ierr);
ierr = MatDestroy(UserCtx->Lap); CHKERRQ(ierr);
ierr = MatDestroy(UserCtx->Poch);
CHKERRQ(ierr);
return(ierr);
}
Program główny - funkcja main() - jest
znów prościutki, co więcej, jest on prawie identyczny z programem rozwiązującym
problem 2x2 z poprzedniego przykładu. To bardzo ważne spostrzeżenie, albowiem
oznacza, że raz nauczywszy się posługiwać SNESem, będziemy mogli używać
go bez kłopotu w innych sytuacjach.
Po wstępnych #include następuje lista
prototypów procedur - zdefiniowanych gdzie indziej, np. w innych plikach
- a wykorzystywanych poniżej. Przydaje się to kompilatorowi dla pełniejszego
sprawdzania poprawności wywołań tych funkcji. W samej funkcji main()
deklarujemy na początku kilka zmiennych lokalnych, w tym obiekt SNESu MojSnes,
macierz (za chwilę okaże się, że powłokową) jakobianu MatJ oraz
kontekst UserCtx, zawierający m.in. macierze laplasjanu i pierwszej
pochodnej, a także wektory rozwiązania i residuum. Zaczynamy klasycznie,
od inicjalizacji PETSc oraz ściągnięcia opcji użytkownika:
PetscInitialize( &argc, &argv,PETSC_NULL,help );
/* Interpretacja opcji uzytkownika */
OptionsGetInt( PETSC_NULL, "-N", &n, &flg );
if( !flg ) n = 3;
OptionsGetScalar( PETSC_NULL, "-L", &L, &flg );
Dalej, tworzymy SNES i kontekst użytkownika
/* Tworzenie potrzebnych obiektow: SNESu, wektorow, macierzy pomocniczych */
SNESCreate(MPI_COMM_WORLD,SNES_NONLINEAR_EQUATIONS,&MojSnes);
USERCTXCreate( &UserCtx, n, L );
Następnie, jak zwykle, "nakarmimy" SNES danymi
zadania. Najpierw przekażemy mu funkcję:
/* Przekaz SNESowi funkcje i wektor residualny */
SNESSetFunction(MojSnes, UserCtx.R, Funkcja, &UserCtx);
Zanim przekażemy SNESowi jakobian, musimy najpierw
utworzyć macierz (powłokową) jakobianu:
MatCreateShell(MPI_COMM_SELF,n,n,n,n,&UserCtx,&MatJ);
MatShellSetOperation(MatJ, MATOP_MULT, MultJakobian );
Kontekstem wykorzystywanym przez macierz jakobianu
będzie oczywiście nasz
UserCtx. Jako procedurę mnożenia przekazujęmy
procedurę
MultJakobian(). Działa ona zgodnie z tym, co wcześniej
wymyśliliśmy w sekcji
"Implementacja
macierzy pochodnej F'(X)" w tym rozdziale.W końcu, możemy przekazać
SNESowi (teraz już dobrze określoną) macierz jakobianu
MatJ i
(prawie pustą) funkcję
Jakobian():
SNESSetJacobian(MojSnes,MatJ,MatJ,Jakobian,&UserCtx);
Dalszy ciąg programu to pełny standard: uwzględniamy
w SNESie opcje linii komend, rozwiązujemy, oglądamy i na koniec sprzątamy
po sobie. Dlatego na tym zakończymy analizę funkcji main(). Podkreślmy
raz jeszcze, że praktycznie nie różni się ona od funkcji main() z
poprzedniego przykładu!
Funkcje związane z kontekstem USERCTX
Są to USERCTXCreate() i USERCTXDestroy()
- mało skomplikowane funkcje, więc komentarze w pliku źródłowym w
wystarczającym stopniu opisują ich działanie. Zwróćmy jedynie uwagę na
to, że generowaniem elementów macierzy laplasjanu i pierwszej pochodnej
zajmują się zewnętrzne funkcje GenerujLaplasjan1D() i GenerujPochodnaCentralna1D(),
zdefiniowane gdzie indziej zapewne w jeszcze jednym pliku, którego tu wcale
nie zamieszczamy, zostawiając opracowanie tych procedur Czytelnikowi jako
ćwiczenie.
Ćwiczenia
-
Uzupełnij powyższy program o funkcję GenerujLaplasjan1D() i przetestuj
dla rozmaitych prawych stron zagadnienia konwekcji-dyfuzji.
-
Rozwiąż zagadnienie konwekcji-dyfuzji w obszarze dwuwymiarowym.
-
Rozwiąż zagadnienie , z jednorodnymi
warunkami brzegowymi Dirichleta.
-
Zastanów się, jak wykorzystać moc SNESu do implementacji metody typu Predyktor-Korektor
rozwiązywania układu równań różniczkowych zwyczajnych. (Na każdym kroku
schematu PC będziemy musieli rozwiązać równanie nieliniowe!)
-
Trudne, ale warto... Napisz program, rozwiązujący równanie Naviera-Stokesa
w obszarze L-kształtnym:
Wybór schematu numerycznego dla aproksymacji tego
równania nie jest trywialny.
Copyright (C) Marcin Bojko i Piotr Krzyzanowski,
1998.
Ostatnia modyfikacja: 12.X.98.