過渡問題の積分のためのモデルツール

時間積分法は各反復で求解する問題を記述することによってmodelオブジェクトを使用して直接書くことができますが，modelオブジェクトはそのような方法の記述を容易にするためのいくつかの基本的なツールを備えています．これらのツールは，次の基本原則に基づいています．

• モデルの元の変数は，現在の時間ステップ（例えば，ステップn）で解かれるべきシステムの状態を表します．これは，中点法の場合でも，主にシステムの異なる変数に異なる方法を適用する必要がある場合，すべての変数は固有の時間ステップでシステムを記述する必要があるためです．

• 前のタイムステップにおけるシステムの状態を表すために，一部のデータがモデルに追加されます．従来のワンステップ方式(現時点では，1ステップのスキームのみが提供されています)では，直前のタイムステップのみが保存されます．例えば， u が変数である場合(このようにステップnで表される)， Previous_u ， Previous2_u ， Previous3_u は，前の時間ステップ(ステップn-1，n-2およびn-3)における変数の状態を表すデータとなります．

• いくつかの中間変数は，時間導関数（および2次問題の2次時間導関数）を表すためにモデルに追加されます．例えば， u が変数であれば， Dot_u は u の一次時間導関数を表し， Dot2_u は2次のものを表します．モデル内のこれらの変数を参照して，要素を追加するか，GWFL，汎用弱形式言語で使用することができます．しかし，これらは独立変数とはみなされず，時間積分法によって対応する元の変数に（Affine方式で）リンクされます．アルゴリズムの大部分は前の時間ステップで時間微分を必要とし， Previous_Dot_u と Previous_Dot2_u データをモデルに加えることが可能です．

• 異なる時間積分法は，モデルの各変数に適用することができます．ほとんどの場合，乗数変数，より一般的には時間微分が使用されない変数は時間積分法を必要としないことに注意してください．

• データ t は時間パラメータを表し，（GWFLで，またはいくつかの要素のパラメータとして）使用することができます．システムの構築前に，データ t は自動的に時間ステップ n に更新されます．

• 各イテレーションで解くべき問題は，速度と加速度が導入された中間変数によって表されるその自然（弱）定式化における過渡問題の定式化に対応します．例えば，問題をGWFLに訳すと

  \[\dot{u}(t,x) - \Delta u(t,x) = \sin(t)\]

  の弱形式言語による表現は以下のように記述することができます．

  Dot_u*Test_u + Grad_u.Grad_Test_u - sin(t)*Test_u
  

  （もちろん，このような状況では，効率の理由から線形要素の使用が望ましいです．）

• 実装されたすべての1ステップ法では，時間の問題は1次と2次の両方の問題（または準静的問題）でも反復から変更することができます．

• 2次時間問題（1次時間）に対する計算法は，2次または1次の時間に，あるいは準静的問題（1次または準静的問題へ）に適用することもできます．速度/変位に対応する初期データを計算法の次数に関して初期化しなければならないという点を除いて，1次問題の理論は，2次問題の時間問題に適用することはできません．

1次の陰的Theta法問題

次の問題を見てください．

\[M\dot{U} = F(U)\]

ここで， F(U) は非線形である可能性があります（結合問題の他の変数に依存するかもしれません） \(dt\) は時間ステップ， \(V = \dot{U}\) と \(U^n, V^n\) は時刻 \(ndt\) における \(U, V\) の近似です．Theta法は次の式になります．

\[\begin{split}\left\{ \begin{array}{l} U^n = U^{n-1} + dt(\theta V^n + (1-\theta) V^{n-1}), \\ MV^n = F(U^n), \end{array}\right.\end{split}\]

Theta法のパラメータである \(\theta \in (0, 1]\) （ \(\theta = 0\) 法は前方Euler法に対応し，陰解法ではありません）そして， \(U^{n-1}, V^{n-1}\) が与えられます．

最初の時間ステップの前に \(U^0\) を初期化する必要があります，しかし， \(V^0\) も必要です（ \(\theta = 1\) の場合を除く）．この例では，それは \(M^{-1}F(U^0)\) に対応するはずです．汎用的な結合問題 \(M\) は特異であり，汎用的な事前計算 \(V^0\) を得ることは困難です．したがって \(V^0\) もまた準備しなければいけません．代替的に（以下を参照），modelオブジェクト（および標準解）は（非常に）小さな時間ステップで Backward Euler 法を適用することによりそれらを評価するための平均値を計算します．

次の式は時間微分により導き出すことができます：

\[V^n = \frac{U^n - U^{n-1}}{\theta dt} - \frac{1-\theta}{\theta}V^{n-1}\]

この方法をモデルの変数 "u" に適用すると，次の Affine 依存変数がモデルに追加されます:

"Dot_u"

これは変数の時間微分を表し，いくつかのブリック要素定義で使用することができます．

以下のデータも追加されています．

"Previous_u", "Previous_Dot_u"

これは前回の時間ステップにおける "u" と "Dot_u" の値に対応します．

解く前にデータ "Previous_u" （例 \(U^0\) を除く）を初期化する必要があります．ここでも， "Previous_Dot_u" は，次のセクションで説明するように，初期化されるか事前計算される必要があります．したがって， "Dot_u" のAffine依存性は，次のように与えられます．

Dot_u = (u - Previous_u)/(theta*dt) - Previous_Dot_u*(1-theta)/theta

つまり， "Dot_u" は構築時に "u" （ \(1/(\theta*dt)\) を掛けたもの）と前の時間ステップに依存する残りの定数項による表現で置き換えられます．変数へのこの方法の追加は次のように行います．:

add_theta_method_for_first_order(model &md, const std::string &varname, scalar_type theta);

速度/加速度の事前計算

大部分の時間積分法（例えば，後方Euler法を除く）は，最初の時間ステップの前に第1または第2の時間微分の事前計算を必要とします（例えば，1次元問題のためのTheta法，第2次問題のための \(A^0\) ...）．

ユーザーはこれらの導関数を初期化するか，モデルに自動的に近似するよう依頼するか選択できます．

補完的な導関数を近似するために（現時点で）使用されている方法は次の通りです．

\[M\dot{U} = F(U)\]

Theta法（前のセクションを参照）を使用します． \(V_0\) を近似するために，Theta法は，非常に小さな時間ステップで， \(\theta = 1\) （すなわち，後方Euler法）に適用されます．これは，後退Eulerが初期時間微分を必要としないので可能です．その後，非常に小さな時間ステップの終了時の後方Eulerのおかげで計算された時間微分は，単純に初期時間微分の近似値として使用されます．

モデル md について，以下の命令は

model.perform_init_time_derivative(ddt);
standard_solve(model, iter);

初期時間導関数の近似を自動的に実行することを可能にします．パラメータ ddt は，近似を実行するために使用される小さな時間ステップに対応します．典型的には，dtt = dt / 20 を使用します，ここで， dt は過渡問題を近似するために使用される時間ステップです（以下の例を参照）．

2次問題の陰的Theta法

次の問題を見てください．

\[M\ddot{U} = F(U)\]

\(dt\) は時間ステップ， \(V = \dot{U}\), \(A = \ddot{U}\) と \(U^n, V^n, A^n\) は時刻 \(ndt\) における \(U, V, A\) の近似です，Theta法は次の式になります．

\[\begin{split}\left\{ \begin{array}{l} U^n = U^{n-1} + dt(\theta V^n + (1-\theta) V^{n-1}), \\ V^n = V^{n-1} + dt(\theta A^n + (1-\theta) A^{n-1}), \\ MA^n = F(U^n), \end{array}\right.\end{split}\]

Theath法のパラメータは \(\theta \in (0, 1]\) であり( \(\theta = 0\) ，計算法は前方Euler法に対応し，陰解法ではありません） \(U^{n-1}, V^{n-1}, A^{n-1}\) が与えられています．

最初のステップでは， \(U^0, V^0\) が与えられ， \(A^0\) が与えられるか，あらかじめ計算されている必要があります（ \(\theta = 1\) を除く）．

次の式は時間微分により導き出すことができます：

\[ \begin{align}\begin{aligned}V^n = \frac{U^n - U^{n-1}}{\theta dt} - \frac{1-\theta}{\theta}V^{n-1}\\A^n = \frac{U^n - U^{n-1}}{\theta^2 dt^2} - \frac{1}{\theta^2dt}V^{n-1} - \frac{1-\theta}{\theta}A^{n-1}\end{aligned}\end{align} \]

この法をモデルの変数 "u"に適用すると，以下のAffine依存変数がモデルに追加されます:

"Dot_u", "Dot2_u"

変数の1次および2次の時間微分を表し，いくつかのブリック要素定義で使用できます．

以下のデータも追加されています．

"Previous_u", "Previous_Dot_u", "Previous_Dot2_u"

前回の時間ステップにおける "u"， "Dot_u"， "Dot2_u"の値に対応します．

解く前に，データ "Previous_u"と "Previous_Dot_u"（例では \(U^0\) ）を初期化し，" Previous_Dot2_u "を初期化または事前計算する必要があります \(\theta = 1\) ）．したがって，Affine依存性は次のように与えられます．

Dot_u = (u - Previous_u)/(theta*dt) - Previous_Dot_u*(1-theta)/theta
Dot2_u = (u - Previous_u)/(theta*theta*dt*dt) - Previous_Dot_u/(theta*theta*dt) - Previous_Dot2_u*(1-theta)/theta

変数へのこの計算法の追加は次のように行います．

add_theta_method_for_second_order(model &md, const std::string &varname,
                                  scalar_type theta);

2次問題の陰的Newmark法

次の問題を見てください．

\[M\ddot{U} = F(U)\]

\(dt\) は時間ステップ， \(V = \dot{U}\), \(A = \ddot{U}\) と \(U^n, V^n, A^n\) は時刻 \(ndt\) における \(U, V, A\) の近似です，Theta法は次の式になります．

\[\begin{split}\left\{ \begin{array}{l} U^n = U^{n-1} + dtV^n + \frac{dt^2}{2}(2\beta V^n + (1-2\beta) V^{n-1}), \\ V^n = V^{n-1} + dt(\gamma A^n + (1-\gamma) A^{n-1}), \\ MA^n = F(U^n), \end{array}\right.\end{split}\]

\(\beta \in (0, 1]\) と \(\gamma \in [1/2, 1]\) はNewmark法のパラメータであり， \(U^{n-1}, V^{n-1}, A^{n-1}\) が与えられています．

最初のステップでは， \(U^0, V^0\) が与えられているとします．そして( \(\beta = 1/2, \gamma = 1\) の場合を除いて) \(A^0\) が与えられているか事前計算されているとします．

次の式は時間微分により導き出すことができます：

\[ \begin{align}\begin{aligned}V^n = \frac{\gamma}{\beta dt}(U^n - U^{n-1}) + \frac{\beta-\gamma}{\beta}V^{n-1} + dt(1-\frac{\gamma}{2\beta})A^{n-1}\\A^n = \frac{U^n - U^{n-1}}{\beta dt^2} - \frac{1}{\beta dt}V^{n-1} - (1/2-\beta)A^{n-1}\end{aligned}\end{align} \]

この法をモデルの変数 "u"に適用すると，以下のAffine依存変数がモデルに追加されます:

"Dot_u", "Dot2_u"

変数の1次および2次の時間微分を表し，いくつかのブリック要素定義で使用できます．

以下のデータも追加されています．

"Previous_u", "Previous_Dot_u", "Previous_Dot2_u"

前回の時間ステップにおける "u"， "Dot_u"， "Dot2_u"の値に対応します．

最初の球解の前に，データ "Previous_u"と "Previous_Dot_u" （この例では \(U^0\) に対応）を初期化する必要があります．データ "Previous_Dot2_u" は与えられているか，事前計算されています（ 速度/加速度の事前計算 を参照してください \(\beta = 1/2, \gamma = 1\) を除きます）．

変数へのこの計算法の追加は次のように行います．

add_Newmark_scheme(model &md, const std::string &varname,
                   scalar_type beta, scalar_type gamma);

Houbolt陰解法

次の問題を見てください．

\[(K+\frac{11}{6 dt}C+\frac{2}{dt^2}M) u_{n} = F_{n} + (\frac{5}{dt^2} M + \frac{3}{ dt} C) u_{n-1} - (\frac{4}{dt^2} M + \frac{3}{2 dt} C) u_{n-2} + (\frac{1}{dt^2} M + \frac{1}{3 dt} C) u_{n-3}\]

ここで， \(dt\) は時間ステップを意味し， \(M\) は "Dot2_u" の項の行列， \(C\) は "Dot_u" の項の行列， \(K\) は "u" の項の行列を意味します．したがって，affine 依存関係は次のようになります．

Dot_u  = 1/(6*dt)*(11*u-18*Previous_u+9*Previous2_u-2*Previous3_u)
Dot2_u = 1/(dt**2)*(2*u-5*Previous_u+4*Previous2_u-Previous3_u)

この法をモデルの変数 "u"に適用すると，以下のAffine依存変数がモデルに追加されます:

"Dot_u", "Dot2_u"

変数の1次および2次の時間微分を表し，いくつかのブリック要素定義で使用できます．

以下のデータも追加されています．

"Previous_u", "Previous2_u", "Previous3_u"

これは，時間ステップ n-1, n-2 n-3 における "u" の値に対応します．

ソルバを実行する前に，(例では \(U^0\) に対応する)データ "Previous_u" ， "Previous2_u" ， "Previous3_u" を初期化する必要があります．

変数へのこの計算法の追加は次のように行います．

add_Houbolt_scheme(model &md, const std::string &varname);

過渡項

これまでのセクションで説明したように，手法が適用される変数の時間微分を表すために，いくつかの中間変数がモデルに追加されています．再度， "u" がそのような変数である場合， "Dot_u" は，使用された法によって近似された "u" の時間微分を表します．

これはまた，過渡項を表現するために "Dot_u" ( および2次の時間問題において "Dot2_u" ）を使用できることを意味します．GWFLでは，

\[\int_{\Omega} \dot{u} v dx\]

次のように単純に表現することができます．

Dot_u*Test_u

同様に， GetFEM は "Dot_u" に適用できます．たとえば，次のような場合です．

getfem::add_mass_brick(model, mim, "Dot_u");

同じ過渡項を追加します．

非常に重要： "Dot_u" のようなAffine依存変数に適用された既存のモデルブリック要素を追加する場合，対応する試験関数が，対応する元の変数（ここでは "Test_u" ）のものであると常に仮定されます．言い換えれば，速度に対応する試行変数である "Test_Dot_u" は使用されません．これは，試験関数が元の変数に対応するように，元の変数の項で問題を解決するために行われた選択に対応します．

Kelvin-Voigt線形化粘度項を説明するために使用できるモデルブリック要素の別の例は，線形弾性ブリック要素となります．

getfem::add_isotropic_linearized_elasticity_brick(model, mim, "Dot_u", "lambda_viscosity", "mu_viscosity");

2次の過渡の弾性の問題に適用します．

時間ステップのシーケンスに対する計算

一般に，異なる時間ステップでの解法は次の形式をとります．

for (scalar_type t = 0.; t < T; t += dt) { // time loop

  // Eventually compute here some time dependent terms

  iter.init();
  getfem::standard_solve(model, iter);

  // + Do something with the solution (plot or store it)

  model.shift_variables_for_time_integration();
}

手法の呼び出しは次の通りです．

model.shift_variables_for_time_integration();

変数の現在の値（例えば u と Dot_u ）を前ステップの値（ Previous_u と Previous_Dot_u ）にコピーするので，2つの時間ステップの間に呼び出す必要があります．

境界条件

もちろん，標準的な境界条件は，未知の様々な変数に通常適用することができます．デフォルトでは，Dirichlet，Neumannまたは接触境界条件を未知数に適用するということは，現在の時間ステップnで変数に条件が規定されていることを意味します．

小さな例：熱方程式

GetFEM ディストリビューションの完全なコンパイル可能なプログラム tests/heat_equation.cc がテストプログラムに対応しています．matlabインタフェースの2次時間ステップ問題の例については， /interface/tests/matlab/demo_wave_equation.m を参照してください．

mf_u と mim が有効なメッシュに定義された有効な有限要素と積分法であると仮定すると，次のコードは単位拡散係数を仮定したメッシュ上の温度の進展を近似します．

getfem::model model;
model.add_fem_variable("u", mf_u, 2); // Main unknown of the problem

getfem::add_generic_elliptic_brick(model, mim, "u"); // Laplace term

// Volumic source term.
getfem::add_source_term_generic_assembly_brick(model, mim, "sin(t)*Test_u");


// Dirichlet condition.
getfem::add_Dirichlet_condition_with_multipliers
    (model, mim, "u", mf_u, DIRICHLET_BOUNDARY_NUM);

// transient part.
getfem::add_theta_method_for_first_order(model, "u", theta);
getfem::add_mass_brick(model, mim, "Dot_u");

gmm::iteration iter(residual, 0, 40000);

model.set_time(0.);        // Init time is 0 (not mandatory)
model.set_time_step(dt);   // Init of the time step.

// Null initial value for the temperature.
gmm::clear(model.set_real_variable("Previous_u"));

// Automatic computation of Previous_Dot_u
model.perform_init_time_derivative(dt/20.);
iter.init();
standard_solve(model, iter);


// Iterations in time
for (scalar_type t = 0.; t < T; t += dt) {

  iter.init();
  getfem::standard_solve(model, iter);

  // + Do something with the solution (plot or store it)

  // Copy the current variables "u" and "Dot_u" into "Previous_u"
  // and "Previous_Dot_u".
  model.shift_variables_for_time_integration();
}

陰的/陽的な項

...

陽解法

...

時間ステップ適応

...

準静的問題

...