|
|
|
## Motivation
|
|
|
|
Beim Umsetzen einer Differentialgleichung ist der Standardweg in AMDiS diese zu zerlegen in einzelne Terme, zu klassifizieren als 0.-Orderung, 1.-Ordnung oder 2.-Ordnungs Term und dann alle in Operatoren zu packen, die wiederum dem zu lösenden Problem zugewiesen werden. Beispielsweise findet man folgenden Ausdruck im Navier-Stokes Beispiel:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
Operator *opUGradU2 = new Operator(getFeSpace(i), getFeSpace(i));
|
|
|
|
opUGradU2->addTerm(new Vec2AndPartialDerivative_FOT(
|
|
|
|
densityPhase, prob->getSolution()->getDOFVector(j), j), GRD_PHI);
|
|
|
|
prob->addMatrixOperator(*opUGradU2, i, i);
|
|
|
|
```
|
|
|
|
|
|
|
|
Was genau dieser Operator nun implementiert lässt sich am Namen der Terme nur noch erahnen. Als Ausweg daraus und aus der Notwendigkeit für immer neue Termkonstallationen habe ich folgende Variante implementiert, die die Sytax zumindest einer Zeile etwas aufräumt:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
Operator *opUGradU2 = new Operator(getFeSpace(i), getFeSpace(i));
|
|
|
|
addFOT(opUGradU2, valueOf(densityPhase) *
|
|
|
|
valueOf(prob->getSolution()->getDOFVector(j)), j, GRD_PHI);
|
|
|
|
prob->addMatrixOperator(*opUGradU2, i, i);
|
|
|
|
```
|
|
|
|
|
|
|
|
Dies ist zwar nur unwesentlich kürzer, gibt dafür direkte Auskunft darüber, wie sich der Koeffizient des Terms zusammensetzt, nämlich als Produkt von 2 DOFVectoren (`densityPhase` und `solution[j]`). Die Syntax lässt sich aber nicht nur anwenden beim Füllen von Operatoren, sondern auch für die Integration über Ausdrücke von DOFVectoren und bei der Transformation von DOFVectoren. Dazu aber später mehr. Nachfolgend ist die neue Syntax aufgelistet.
|
|
|
|
|
|
|
|
## Syntaxbeschreibung
|
|
|
|
### Header
|
|
|
|
|
|
|
|
Um die neue Sytax nutzen zu können muss folgende Header-datei eingebunden werden:
|
|
|
|
```c++
|
|
|
|
#include "Expressions.h"
|
|
|
|
```
|
|
|
|
|
|
|
|
### Elementare Ausdrücke
|
|
|
|
In der Folgenden Tabelle ist DV ein Pointer/Reference auf DOFVector<T> mit beliebigem T:
|
|
|
|
|
|
|
|
| Expression | Semantics |
|
|
|
|
| --------------- | ----------------------------- |
|
|
|
|
| `valueOf(DV)` | Auswertung des DOFVectors an Quadraturpunkten |
|
|
|
|
| `gradientOf(DV)`| Auswertung des Gradienten eines DOFVectors an Quadraturpunkten |
|
|
|
|
| `derivativeOf(DV, i)` | Auswertung der partiellen Ableitung des DOFVectors nach der i-ten Koordinate an Quadraturpunkten |
|
|
|
|
| `hessianOf(DV)` | Auswertung der Hesse-Matrix des DOFVectors: `H_ij = d_i(d_j(DV))` (kann nur in `transformDOF`, bzw. `<< ` verwendet werden) |
|
|
|
|
| `laplacianOf(DV)` | Auswertung des Laplace des DOFVectors: `L = sum_i(d_i(d_i(DV)))` (kann nur in `transformDOF`, bzw. `<< ` verwendet werden) |
|
|
|
|
| `componentOf(DV, i)` | Auswertung der i-ten Komponente eines DOFVectors an Quadraturpunkten. <br />Dies macht nur Sinn für DOFVector<Vector<T> > Objekten. |
|
|
|
|
| `X(), X<i>()` | Der Koordinatenvektor, bzw. die i-te Komponente des Koordinatenvektors an den QP. |
|
|
|
|
| `N(b), N(b,i)` | Der äußere Normalenvektor zurgörig zur Fläche mit dem Index b, bzw. deren i-te Komponente. |
|
|
|
|
| `M(), M(i)` | Der Elementnormalenvektor (bei Oberflächengittern), bzw. deren i-te Komponente |
|
|
|
|
| `constant(c), c` | Konstanten (numeric values) sind z.B. double, float, int. Entweder eingeschlossen in constant(), oder direkt |
|
|
|
|
| `constant<C>()` | (Compiletime) Konstanten sind Integer, die beim Differenzierten von Ausdrücken verwendet werden können (s.u.). |
|
|
|
|
| `ref_(c)` | Eine Referenz oder ein Pointer auf einen Wert, der sich dynamisch verändern kann. |
|
|
|
|
|
|
|
|
Der Ausdruck `valueOf(DV)` kann außerdem auf Vektoren und Matrizen von DOFVectoren angewendet werden:
|
|
|
|
|
|
|
|
| Expression | Semantics |
|
|
|
|
| ---------- | --------- |
|
|
|
|
| `valueOf(DV)` | Auswertung des DOFVectors an Quadraturpunkten |
|
|
|
|
| `valueOf(vector<DV<T>>)` | Ergibt an den Quadraturpunkten einen vector<T> mit den Werten der einzelnen DOFVektoren and dieser Stelle. |
|
|
|
|
| `valueOf(matrix<DV<T>>)` | Ergibt an den Quadraturpunkten eine matrix<T> mit den Werten der einzelnen DOFVektoren and dieser Stelle. |
|
|
|
|
| `valueOf<Name>(X)` | Eine Expression, der ein Name zugeordnet ist. `Name` ist dabei ein C++-Typ. Ist kein Name explizit angegeben, wird `_unknown` verwendet. (siehe auch Differenzieren von Ausdrücken) |
|
|
|
|
|
|
|
|
### Operationen
|
|
|
|
In der folgenden Tabelle sind Operationen aufgelistet, über die die elementaren Term (s.o.) kombiniert werden können. Dabei steht `T`, bzw. `Tn` für einen Term und `F` für ein Funktor-Objekt (siehe unten).
|
|
|
|
|
|
|
|
| Expression | Semantics |
|
|
|
|
| ---------- | --------- |
|
|
|
|
| `+, -, *, /` | Elementare arithmetische Ausdrücke |
|
|
|
|
| `pow<p>(T)` | Die p-te Potenz eines Term, wenn Term skalar ist. |
|
|
|
|
| `sqrt(T)` | Die Wurzel des Ausdrucks T |
|
|
|
|
| `exp(T)` | Die Exponentialfunktion angewendet auf das Ergebnis der Auswertung des Terms T. |
|
|
|
|
| `log(T)` | Der natürlich Logarithmus der Auswertung des Terms T. |
|
|
|
|
| `cos(T), sin(T), tan(T)` | Der Cosinus, Sinus, bzw. Tangens der Auswertung des Terms T. |
|
|
|
|
| `acos(T), asin(T), atan(T)` | Der Arkus-Cosinus, Arkus-Sinus, bzw. Arkus Tangens der Auswertung des Terms T. |
|
|
|
|
| `atan2(T1, T2)` | Der Arkus-Tangens 2 = atan(T1 / T2). |
|
|
|
|
| `cosh(T), sinh(T), tanh(T)` | Der Cosinus-Hyperbolicus, Sinus-Hyperbolicus, bzw. Tangens-Hyperbolicus der Auswertung des Terms T. |
|
|
|
|
| `acosh(T), asinh(T), atanh(T)` | Der Arkus Cosinus-Hyperbolicus, Arkus Sinus-Hyperbolicus, Arkus Tangens-Hyperbolicus der Auswertung des Terms T. |
|
|
|
|
| `max(T1, T2), min(T1, T2)` | Das Maximum/Minimum der Auswertungen von T1, bzw. T2. |
|
|
|
|
| `abs_(T), signum(T)` | Der Absolutbetrag von T, bzw. das Vorzeichen (`T < 0 => signum(T) := -1, T>0 => signum(T) := 1, T==0 => signum(T) := 0`) |
|
|
|
|
| `clamp(T, lo, hi)` |
|
|
|
|
| `ceil(T), floor(T)` | die kleines ganze Zahl größer, bzw. die größte ganze Zahl kleiner als Wert von T |
|
|
|
|
| `function_(F, T1, T2, T3, ...)` | Die Anwendung eines Funktors F auf die Auswertung der Terme T1,T2,... | If `T` compares less than `lo`, returns `lo`; otherwise if `hi` compares less than `T`, returns `hi`; otherwise returns `T`
|
|
|
|
|
|
|
|
Ein Funktor ist dabei eine Klasse mit folgender Struktur:
|
|
|
|
```c++
|
|
|
|
struct Functor : public FunctorBase
|
|
|
|
{
|
|
|
|
typedef (...) value_type;
|
|
|
|
int getDegree(int d1, int d2, ...) const { return (...); }
|
|
|
|
|
|
|
|
value_type operator()(const T1::value_type&, const T2::value_type&, ...) const { return (...); }
|
|
|
|
};
|
|
|
|
```
|
|
|
|
|
|
|
|
Die Argumente, die an die Funktion `getDegree` übergeben werden sind die Polynomgrade der Terme:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
getDegree(T1.getDegree(), T2.getDegree(), ...)
|
|
|
|
```
|
|
|
|
|
|
|
|
|
|
|
|
### Vektor-/Matrix-Ausdrücke
|
|
|
|
Für vektor- bzw. matrixwertige Term, wie z.B. `gradientOf(DV)` gibt es eine Reihe von Ausdrücken / Operationen, die im Folgenden aufgelistet sind. Dabei bezeichnet `V` eine vektorwerte Expression und `M` eine matrixwertige Expressions.
|
|
|
|
|
|
|
|
| Expression | Semantics |
|
|
|
|
| ---------- | --------- |
|
|
|
|
| `+, -` | Elementare arithmetische Ausdrücke (Elementweise) |
|
|
|
|
| `unary_dot(V)` | Skalarprodukt eines vektorwertigen Terms mit sich selbst: `result = V^H * V`. |
|
|
|
|
| `dot(V1, V2)` | Skalarprodukt zweier vektorwertigen Terme: `result = V1^H * V2`. |
|
|
|
|
| `one_norm(V)` | Die 1-Norm eines Vektors `result = sum_i(abs(V_i))`. |
|
|
|
|
| `one_norm(M)` | Die 1-Norm einer Matrix `result = max_j(sum_i(abs(M_ij)))`. |
|
|
|
|
| `two_norm(V)` | Die 2-Norm eines Vektors `result = sqrt(V^H * V)`. |
|
|
|
|
| `p_norm<p>(V)` | Die p-Norm eines Vektors `result = [sum_i(abs(v_i)^p)]^(1/p)`. |
|
|
|
|
| `cross(V1, V2)` | Kreuzprodukt zweier vektorwertigen Terme: `result = V1 x V2`. |
|
|
|
|
| `diagonal(V)` | Diagonalmatrix aus einträgen eines vektorwertigen Terms: `matrix = diagonal(V)`. |
|
|
|
|
| `outer(V1, V2)` | Äußeres Produkt (dyadisches Produkt / Tensorprodukt) zweier vektorwertigen Terme: `matrix = V1 * V2^T`. |
|
|
|
|
| `trans(M)` | Das Transponierte eines matrixwertigen Terms: `result = M^T`. |
|
|
|
|
| `at(V, i)` | Zugriff auf eine Komponente eines Vektors: `result = V_i`. |
|
|
|
|
| `at(M, i, j)` | Zugriff auf eine Komponente einer Matrix: `result = M_ij`. |
|
|
|
|
|
|
|
|
### Differenzieren von Ausdrücken
|
|
|
|
Terme / Expressions sind Ausdrücke mit Variablen und Operationssymbolen, die Formal nach einzelnen Variablen differenziert werden können. Variablen sind hierbei Ausdrücke `valueOf<Name>(DV)` denen ein Variablenname zugeordnet ist. Jedem `valueOf` Ausdruck ist standardmäßig der Name `_unknown` zugeordnet. Für skalare Terme (also Terme ohne Vektorausdrücke und vektorwertige Variablen) lässt sich die Differentiation automatisch durchführen.
|
|
|
|
|
|
|
|
Um dies nutzen zu können, muss die Datei `diff_expr.hpp` eingebunden werden:
|
|
|
|
```c++
|
|
|
|
#include "expressions/diff_expr.hpp"
|
|
|
|
```
|
|
|
|
|
|
|
|
| Expression | Semantics |
|
|
|
|
| ---------- | --------- |
|
|
|
|
| `diff<Name>(T)` | Differentiation nach Variable "Name": `d_Name (Term)` |
|
|
|
|
| `diff<Name>(T, U)` | Richtungsableitung nach Variable "Name", in Richtung `U`: `d_Name (Term)[U]` |
|
|
|
|
| `simplify(T)` | Vereinfachen eines Terms. Insbesondere werden Multiplikationen mit 0 und mit 1 entfernt und verschachtelte Operationen zwischen Konstanten zusammengefasst. |
|
|
|
|
|
|
|
|
### Anwendung
|
|
|
|
Diese verallgemeinerten Terme können nun beim Assemblieren, integrieren und transformieren von DOFVectoren verwendet werden: In der Folgenden Tabelle wird `Operator` für einen Pointer/Reference auf ein Operator Objekt verwendet und Term steht für einen Ausdruck, der sich aus den obigen Zutaten zusammensetzt. `Flag` ist entweder `GRD_PHI`, oder `GRD_PSI`, je nachdem, ob die Ableitung bei einem 1.-Ordnungs Term auf der Ansatzfunktion `u`, oder Testfunktion `v` steht.
|
|
|
|
|
|
|
|
| Expression | Semantics |
|
|
|
|
| ---------- | --------- |
|
|
|
|
| `addZOT(Operator, Term)` | Ein 0.-Ordnungs Term `(Term * u, v)` |
|
|
|
|
| `addFOT(Operator, Term, Flag)` | Ein 1.-Ordnungs Term `(Term * grad(u), v)`, bzw. `Term * u, grad(v))`, wenn der `value_type` von `Term` einem `WorldVector<double>` entspricht, bzw. `(_1_*Term * grad(u), v)`, bzw. `(_1_ * Term * u, grad(v))`, mit `_1_=(1,1,1)^T`, wenn der `value_type` von `Term` einem `double` entspricht. |
|
|
|
|
| `addFOT(Operator, Term, i, Flag)` | Ein 1.-Ordnungs Term `(Term * d_i(u), v)`, bzw. `(Term * u, d_i(v))` |
|
|
|
|
| `addSOT(Operator, Term)` | Ein 2.-Ordnungs Term `(Term * grad(u), grad(v))`, wobei `value_type` von `Term` einem `WorldMatrix<double>`, bzw. `double` entsprechen muss. |
|
|
|
|
| `addSOT(Operator, Term, i, j)` | Ein 2.-Ordnungs Term `(Term * d_j(u), d_i(v))`>, wobei `value_type` von `Term` einem `double` entsprechen muss. |
|
|
|
|
|
|
|
|
Für die Berechnung eines Integrals über einen Ausdruck von DOFVectoren kann folgender Befehl verwendet werden:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
integrate(Term);
|
|
|
|
```
|
|
|
|
|
|
|
|
und für die Transformation eines DOFVectors folgendes:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
DV << Term;
|
|
|
|
```
|
|
|
|
bzw.
|
|
|
|
```c++
|
|
|
|
transformDOF(Term, DV);
|
|
|
|
```
|
|
|
|
Wobei jeweils DV kein DOFVector-Pointer sein kann.
|
|
|
|
|
|
|
|
## Beispiele
|
|
|
|
### Bilinearform
|
|
|
|
|
|
|
|
Folgende Bilinearform soll assembliert werden:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Mit Hilfe der oben definierten Ausdrücke lässt sich dies nun wie folgt schreiben:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
Operator *op = new Operator(feSpace, feSpace);
|
|
|
|
addZOT(op, (1.0/epsilon) * (pow<2>(valueOf(phi)) - 1.0) );
|
|
|
|
addSOT(op, max(1.e-5, valueOf(phi) + 1.0) )
|
|
|
|
prob->addMatrixOperator(*op, 0, 0);
|
|
|
|
```
|
|
|
|
|
|
|
|
### Integration und Interpolation
|
|
|
|
Danach sollen folgende Integrale und ein Double-well berechnet werden:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Dazu können wieder die Expressions verwendet werden:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
double mean = integrate( 0.5 * (valueOf(c) + 1.0) );
|
|
|
|
double length = integrate( two_norm(gradientOf(c)) );
|
|
|
|
|
|
|
|
DOFVector<double> DW(c.getFeSpace(), "Double-Well");
|
|
|
|
DW << 0.5*( pow<2>(valueOf(c)) * pow<2>(1.0 - valueOf(c)) + unary_dot(gradientOf(c)) );
|
|
|
|
```
|
|
|
|
|
|
|
|
### Funktoren
|
|
|
|
Komplexe Ausdrücken können auch zu einem Funktor zusammen gefasst werden. Dazu wird erst das Funktor-Objekt erstellt und dieses dann über den Term `function_` in einem Expression eingebunden:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Mit `a_i` sind Koeffizienten und `U={U_i}` ein `DOFVector<WorldVector<double> >`. Dieser Term kann umgesetzt werden, mittels
|
|
|
|
|
|
|
|
```c++
|
|
|
|
struct F : FunctorBase {
|
|
|
|
typedef double value_type;
|
|
|
|
F(WorldVector<double>& alpha_) : alpha(alpha_) {}
|
|
|
|
|
|
|
|
int getDegree(int d0) const {
|
|
|
|
return d0;
|
|
|
|
}
|
|
|
|
value_type operator()(const WorldVector<double>& U) const {
|
|
|
|
double result = 0.0;
|
|
|
|
for (int i = 0; i < U.getSize(); i++)
|
|
|
|
result += alpha[i]*U[i];
|
|
|
|
return result;
|
|
|
|
}
|
|
|
|
private:
|
|
|
|
WorldVector<double> alpha;
|
|
|
|
};
|
|
|
|
|
|
|
|
int main() {
|
|
|
|
// ... //
|
|
|
|
DOFVector<WorldVector<double> > U(feSpace, "U");
|
|
|
|
DOFVector<double> V(feSpace, "V");
|
|
|
|
|
|
|
|
WorldVector<double> alpha = (...);
|
|
|
|
U << (...)
|
|
|
|
|
|
|
|
V << function_(F(alpha), valueOf(U));
|
|
|
|
}
|
|
|
|
```
|
|
|
|
|
|
|
|
Bisher wird in AMDiS anstelle von Objekten zur Basis `FunctorBase` die AbstractFunction verwendet. Um diese nutzen zu können, gibt es Wrapper, die diese implizit zu einem Funktor Konvertieren. Für den Spezialfall der Auswertung an Koordinaten gibt es zusätzlich einen Shortcut:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
struct F : AbstractFunction<double, WorldVector<double> > {
|
|
|
|
F() : AbstractFunction<double, WorldVector<double> >(1) {}
|
|
|
|
|
|
|
|
double operator()(const WorldVector<double>& x) const {
|
|
|
|
return cos(x[0]);
|
|
|
|
}
|
|
|
|
};
|
|
|
|
|
|
|
|
int main() {
|
|
|
|
// ... //
|
|
|
|
DOFVector<double> V(feSpace, "V");
|
|
|
|
|
|
|
|
V << function_(wrap(new F), X());
|
|
|
|
// bzw.
|
|
|
|
V << eval(new F)
|
|
|
|
|
|
|
|
// und auch verschobene Auswertungen sind möglich:
|
|
|
|
WorldVector<double> shift = (...)
|
|
|
|
V << function_(wrap(new F), X() + shift);
|
|
|
|
}
|
|
|
|
```
|
|
|
|
|
|
|
|
### Differentiation
|
|
|
|
Insbesondere in Newtonverfahren für nichtlineare PDEs und Rosenbrock-Verfahren für die Zeitintegration (nichtlinearer Gleichungen), wird die Jacobimatrix eines Differentialoperators benötigt. Den kann man häufig auf die Ableitung der Koeffizientefunktionen zurückführen. Siehe aus Motivation auch das Nonlinear Demo [Nonlinear example](nonlinear_example). Grundlage der Implementierung könnte sein, dass die Terme automatisch differenziert werden Als Beispiel soll hier ein einfaches Polynom als Term dienen:
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Um nach `u` zu differenzieren, müssen wir dem zugrundeliegenden DOFVector einen Namen geben:
|
|
|
|
|
|
|
|
```c++
|
|
|
|
struct u_ {};
|
|
|
|
struct v_ {};
|
|
|
|
...
|
|
|
|
DOFVector<double> U(...);
|
|
|
|
DOFVector<double> V(...);
|
|
|
|
auto expr = pow<4>(valueOf<u_>(U));
|
|
|
|
|
|
|
|
auto diff_expr = diff<u_>(expr); // = 4 * u^3
|
|
|
|
auto diff_expr2 = diff<u_>(expr, valueOf<v_>(V)); // 4 * u^3*v
|
|
|
|
``` |
|
|
|
\ No newline at end of file |
