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:
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

  1. Uzupełnij powyższy program o funkcję GenerujLaplasjan1D() i przetestuj dla rozmaitych prawych stron zagadnienia konwekcji-dyfuzji.
  2. Rozwiąż zagadnienie konwekcji-dyfuzji w obszarze dwuwymiarowym.
  3. Rozwiąż zagadnienie , z jednorodnymi warunkami brzegowymi Dirichleta.
  4. 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!)
  5. Trudne, ale warto... Napisz program, rozwiązujący równanie Naviera-Stokesa w obszarze L-kształtnym:
  6. Wybór schematu numerycznego dla aproksymacji tego równania nie jest trywialny.

 
Pierwsze zadanie nieliniowe w PETSc Spis treści Co jeszcze warto wiedzieć o SNESie?
 

Copyright (C) Marcin Bojko i Piotr Krzyzanowski, 1998.
Ostatnia modyfikacja: 12.X.98.