[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Tetsuo Koyama |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sun, 3 Nov 2019 03:01:21 -0500 (EST) |
branch: devel-tetsuo-hht_method
commit 84d2b7b37b1d529433b7226bd6fe709ad995c5f5
Author: Tetsuo Koyama <address@hidden>
Date: Sun Nov 3 14:21:13 2019 +0900
:package: ADD HHT scheme
---
.../source/userdoc/model_time_integration.rst | 13 +-
interface/src/gf_model_set.cc | 17 ++-
interface/tests/python/demo_wave_equation.py | 21 ++-
src/getfem/getfem_models.h | 3 +
src/getfem_models.cc | 151 ++++++++++++++++++++-
5 files changed, 195 insertions(+), 10 deletions(-)
diff --git a/doc/sphinx/source/userdoc/model_time_integration.rst
b/doc/sphinx/source/userdoc/model_time_integration.rst
index 5180564..d0213ca 100644
--- a/doc/sphinx/source/userdoc/model_time_integration.rst
+++ b/doc/sphinx/source/userdoc/model_time_integration.rst
@@ -265,7 +265,7 @@ For a problem which reads
.. math::
- \left(\left(1+\alpha\right)K+\left(1+\alpha\right)\dfrac{\gamma}{\beta\Delta
t}C+\dfrac{1}{\beta\Delta
t^{2}}M\right)U^{n}=F^{n}+M\left(\left(\dfrac{1}{2\beta}-1\right)A^{n-1}+\dfrac{1}{\beta\Delta
t}V^{n-1}+\dfrac{1}{\beta\Delta
t^{2}}U^{n-1}\right)+\left(1+\alpha\right)C\left[\left(\dfrac{\gamma}{2\beta}-1\right)\Delta
tA^{n-1}+\left(\dfrac{\gamma}{\beta}-1\right)V^{n-1}+\dfrac{\gamma}{\beta\Delta
t}U^{n-1}\right]+\alpha KU^{n-1}+\alpha CV^{n-1}
+ \left(\left(1+\alpha\right)K+\left(1+\alpha\right)\dfrac{\gamma}{\beta\Delta
t}C+\dfrac{1}{\beta\Delta
t^{2}}M\right)U^{t}=F^{t+\alpha*dt}+M\left(\left(\dfrac{1}{2\beta}-1\right)A^{t-dt}+\dfrac{1}{\beta\Delta
t}V^{t-dt}+\dfrac{1}{\beta\Delta
t^{2}}U^{t-dt}\right)+\left(1+\alpha\right)C\left[\left(\dfrac{\gamma}{2\beta}-1\right)\Delta
tA^{t-dt}+\left(\dfrac{\gamma}{\beta}-1\right)V^{t-dt}+\dfrac{\gamma}{\beta\Delta
t}U^{t-dt}\right]+\alpha KU^{t-dt}+\alpha CV^{t-dt}
where :math:`dt` means a time step, :math:`M` the matrix in term of "Dot2_u",
:math:`C` the matrix in term of "Dot_u" and :math:`K` the matrix in term of "u".
@@ -286,7 +286,7 @@ The following data are also added::
which correspond to the values of "u", "Dot_u" and "Dot2_u" at the previous
time step.
-Before the first solve, the data "Previous_u" and "Previous_Dot_u"
(corresponding to :math:`U^0` in the example) have to be initialized. The data
"Previous_Dot2_u" is to be given or precomputed (see
:ref:`precomp_time_der_section` and except for :math:`\alpha = 0, \beta = 1/2,
\gamma = 1`).
+Before the first solve, the data "Previous_u" and "Previous_Dot_u"
(corresponding to :math:`U^0` in the example) have to be initialized. The data
"Previous_Dot2_u" is to be given or precomputed (see
:ref:`precomp_time_der_section` and except for :math:`\beta = 1/2, \gamma = 1`).
The addition of this scheme to a variable is to be done thanks to::
@@ -294,6 +294,15 @@ The addition of this scheme to a variable is to be done
thanks to::
add_Hilber_Hughes_Taylor_scheme(model &md, const std::string &varname,
scalar_type alpha, scalar_type beta,
scalar_type gamma);
+.. warning::
+ You have to set exact time of external force :math:`F^{t+\alpha*dt}` when
you use this scheme.
+
+such as::
+
+ // Volumic source term.
+ getfem::add_source_term_generic_assembly_brick(model, mim,
"sin(t+alpha*dt)*Test_u");
+
+
Transient terms
***************
diff --git a/interface/src/gf_model_set.cc b/interface/src/gf_model_set.cc
index 8f1ad05..711a65b 100644
--- a/interface/src/gf_model_set.cc
+++ b/interface/src/gf_model_set.cc
@@ -2860,7 +2860,7 @@ void gf_model_set(getfemint::mexargs_in& m_in,
);
/*@SET ('add Newmark scheme', @str varname, @scalar beta, @scalar gamma)
- Attach a theta method for the time discretization of the variable
+ Attach a Newmark method for the time discretization of the variable
`varname`. Valid only if there is at most second order time derivative
of the variable. @*/
sub_command
@@ -2871,7 +2871,7 @@ void gf_model_set(getfemint::mexargs_in& m_in,
getfem::add_Newmark_scheme(*md, varname, beta, gamma);
);
- /*@SET ('add_Houbolt_scheme', @str varname)
+ /*@SET ('add Houbolt scheme', @str varname)
Attach a Houbolt method for the time discretization of the variable
`varname`. Valid only if there is at most second order time derivative
of the variable @*/
@@ -2881,6 +2881,19 @@ void gf_model_set(getfemint::mexargs_in& m_in,
getfem::add_Houbolt_scheme(*md, varname);
);
+ /*@SET ('add Hilber Hughes Taylor scheme', @str varname, @scalar alpha,
@scalar beta, @scalar gamma)
+ Attach a Hilber Hughes Taylor method for the time discretization of the
variable
+ `varname`. Valid only if there is at most second order time derivative
+ of the variable. @*/
+ sub_command
+ ("add Hilber Hughes Taylor scheme", 4, 4, 0, 0,
+ std::string varname = in.pop().to_string();
+ double alpha = in.pop().to_scalar();
+ double beta = in.pop().to_scalar();
+ double gamma = in.pop().to_scalar();
+ getfem::add_Hilber_Hughes_Taylor_scheme(*md, varname, alpha, beta,
gamma);
+ );
+
/*@SET ('disable bricks', @ivec bricks_indices)
Disable a brick (the brick will no longer participate to the
building of the tangent linear system).@*/
diff --git a/interface/tests/python/demo_wave_equation.py
b/interface/tests/python/demo_wave_equation.py
index a56505e..eaf2c28 100644
--- a/interface/tests/python/demo_wave_equation.py
+++ b/interface/tests/python/demo_wave_equation.py
@@ -27,9 +27,9 @@
also a good example of use of GetFEM++.
"""
+import os
import numpy as np
import getfem as gf
-import os
NX = 10
m = gf.Mesh('cartesian', np.arange(0., 1.+1./NX,1./NX),
@@ -70,6 +70,7 @@ T = 5.0;
# Set dt smaller to fix Newmark and Houbolt method.
# dt = 0.001;
dt = 0.025;
+alpha = 0.00;
beta = 0.25;
gamma = 0.5;
@@ -78,6 +79,7 @@ md.add_Houbolt_scheme('u1')
md.add_Hilber_Hughes_Taylor_scheme('u2', alpha, beta, gamma)
md.add_mass_brick(mim, 'Dot2_u')
md.add_mass_brick(mim, 'Dot2_u1')
+md.add_mass_brick(mim, 'Dot2_u2')
md.set_time_step(dt)
## Initial data.
@@ -86,6 +88,8 @@ md.set_variable('Previous_Dot_u', V0)
md.set_variable('Previous_u1', U0)
md.set_variable('Previous2_u1', U0)
md.set_variable('Previous3_u1', U0)
+md.set_variable('Previous_u2', U0)
+md.set_variable('Previous_Dot_u', V0)
## Initialisation of the acceleration 'Previous_Dot2_u'
md.perform_init_time_derivative(dt/2.)
@@ -93,7 +97,6 @@ md.solve()
A0 = md.variable('Previous_Dot2_u')
-
os.system('mkdir results');
mf.export_to_vtk('results/displacement_0.vtk', U0)
mf.export_to_vtk('results/velocity_0.vtk', V0)
@@ -106,6 +109,13 @@ mf.export_to_vtk('results1/displacement_0.vtk', U0)
mf.export_to_vtk('results1/velocity_0.vtk', V0)
mf.export_to_vtk('results1/acceleration_0.vtk', A0)
+A0 = md.variable('Previous_Dot2_u2')
+
+os.system('mkdir results2');
+mf.export_to_vtk('results2/displacement_0.vtk', U0)
+mf.export_to_vtk('results2/velocity_0.vtk', V0)
+mf.export_to_vtk('results2/acceleration_0.vtk', A0)
+
## Iterations
n = 1;
for t in np.arange(0.,T,dt):
@@ -127,6 +137,13 @@ for t in np.arange(0.,T,dt):
mf.export_to_vtk('results1/velocity_%d.vtk' % n, V)
mf.export_to_vtk('results1/acceleration_%d.vtk' % n, A)
+ U = md.variable('u2')
+ V = md.variable('Dot_u2')
+ A = md.variable('Dot2_u2')
+ mf.export_to_vtk('results2/displacement_%d.vtk' % n, U)
+ mf.export_to_vtk('results2/velocity_%d.vtk' % n, V)
+ mf.export_to_vtk('results2/acceleration_%d.vtk' % n, A)
+
n += 1
md.shift_variables_for_time_integration()
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 0d72462..370fe84 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -1182,6 +1182,9 @@ namespace getfem {
void add_Houbolt_scheme(model &md, const std::string &varname);
+ void add_Hilber_Hughes_Taylor_scheme(model &md, const std::string &varname,
+ scalar_type alpha, scalar_type beta,
+ scalar_type gamma);
//=========================================================================
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index 8f6a42a..f17b904 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -1395,8 +1395,8 @@ namespace getfem {
scalar_type beta, gamma;
public:
- // V = (U-U0)/(theta*dt) - dt*(1-theta)*V0
- // A = (U-U0)/(theta^2*dt^2) - V0/(theta^2*dt) - dt*(1-theta)*A0
+ // V = (1-gamma/beta)*V0+(1-gamma/(2*beta))*dt*A0+gamma/(beta*dt)*(U-U0)
+ // A = (1-1/(2.0*beta))*A0-1/(beta*dt)*V0+1/(beta*dt*dt)*(U-U0)
virtual void init_affine_dependent_variables(model &md) const {
scalar_type dt = md.get_time_step();
scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
@@ -1429,7 +1429,6 @@ namespace getfem {
md.set_real_constant_part(A));
gmm::add(gmm::scaled(md.real_variable(A0), -a2),
md.set_real_constant_part(A));
-
}
}
@@ -1513,7 +1512,6 @@ namespace getfem {
}
}
-
};
void add_Newmark_scheme(model &md, const std::string &varname,
@@ -1636,7 +1634,152 @@ namespace getfem {
md.add_time_integration_scheme(varname, ptsc);
}
+ // ----------------------------------------------------------------------
+ //
+ // Hilber Hughes Taylor (HHT) method
+ //
+ // ----------------------------------------------------------------------
+
+ class APIDECL Hilber_Hughes_Taylor_scheme
+ : public virtual_time_scheme {
+
+ std::string U, U0, V, V0, A, A0;
+ scalar_type alpha, beta, gamma;
+ public:
+ // V = (1-gamma/beta)*V0+(1-gamma/(2*beta))*dt*A0+gamma/(beta*dt)*(U-U0)
+ // A = (1-1/(2.0*beta))*A0-1/(beta*dt)*V0+1/(beta*dt*dt)*(U-U0)
+ virtual void init_affine_dependent_variables(model &md) const {
+ scalar_type dt = md.get_time_step();
+ scalar_type a0 = scalar_type(1)/(beta*dt*dt), a1 = dt*a0;
+ scalar_type a2 = (scalar_type(1) - scalar_type(2)*beta)
+ / (scalar_type(2)*beta);
+ scalar_type b0 = (scalar_type(1)+alpha)*gamma/(beta*dt);
+ scalar_type b1 = (scalar_type(1)+alpha)*(beta-gamma)/beta - alpha;
+ scalar_type b2 =
(scalar_type(1)+alpha)*dt*(1-gamma/(scalar_type(2)*beta));
+ scalar_type c0 = (scalar_type(1)+alpha);
+ scalar_type c1 = alpha;
+
+ md.set_factor_of_variable(U, c0);
+ md.set_factor_of_variable(V, b0);
+ md.set_factor_of_variable(A, a0);
+ if (md.is_complex()) {
+ gmm::copy(gmm::scaled(md.complex_variable(U0), complex_type(c1)),
+ md.set_complex_constant_part(U));
+ gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(b0)),
+ gmm::scaled(md.complex_variable(V0), complex_type(b1)),
+ md.set_complex_constant_part(V));
+ gmm::add(gmm::scaled(md.complex_variable(A0), complex_type(b2)),
+ md.set_complex_constant_part(V));
+ gmm::add(gmm::scaled(md.complex_variable(U0), -complex_type(a0)),
+ gmm::scaled(md.complex_variable(V0), -complex_type(a1)),
+ md.set_complex_constant_part(A));
+ gmm::add(gmm::scaled(md.complex_variable(A0), -complex_type(a2)),
+ md.set_complex_constant_part(A));
+ } else {
+ gmm::add(gmm::scaled(md.real_variable(U0), -b0),
+ gmm::scaled(md.real_variable(V0), b1),
+ md.set_real_constant_part(V));
+ gmm::add(gmm::scaled(md.real_variable(A0), b2),
+ md.set_real_constant_part(V));
+ gmm::add(gmm::scaled(md.real_variable(U0), -a0),
+ gmm::scaled(md.real_variable(V0), -a1),
+ md.set_real_constant_part(A));
+ gmm::add(gmm::scaled(md.real_variable(A0), -a2),
+ md.set_real_constant_part(A));
+ }
+ }
+
+ // V = (U-U0)/dt (backward Euler for precomputation)
+ // A = (U-U0)/(dt^2) - V0/dt
+ virtual void init_affine_dependent_variables_precomputation(model &md)
+ const {
+ scalar_type dt = md.get_time_step();
+ md.set_factor_of_variable(V, scalar_type(1)/dt);
+ md.set_factor_of_variable(A, scalar_type(1)/(dt*dt));
+ if (md.is_complex()) {
+ gmm::copy(gmm::scaled(md.complex_variable(U0),
+ -complex_type(1)/dt),
+ md.set_complex_constant_part(V));
+ gmm::add(gmm::scaled(md.complex_variable(U0),
+ -complex_type(1)/(dt*dt)),
+ gmm::scaled(md.complex_variable(V0),
+ -complex_type(1)/dt),
+ md.set_complex_constant_part(A));
+ } else {
+ gmm::copy(gmm::scaled(md.real_variable(U0),
+ -scalar_type(1)/dt),
+ md.set_real_constant_part(V));
+ gmm::add(gmm::scaled(md.real_variable(U0),
+ -scalar_type(1)/(dt*dt)),
+ gmm::scaled(md.real_variable(V0),
+ -scalar_type(1)/dt),
+ md.set_real_constant_part(A));
+ }
+ }
+
+ virtual void time_derivative_to_be_initialized
+ (std::string &name_v, std::string &name_previous_v) const {
+ if (beta != scalar_type(0.5) || gamma != scalar_type(1))
+ { name_v = A; name_previous_v = A0; }
+ }
+
+ virtual void shift_variables(model &md) const {
+ if (md.is_complex()) {
+ gmm::copy(md.complex_variable(U), md.set_complex_variable(U0));
+ gmm::copy(md.complex_variable(V), md.set_complex_variable(V0));
+ gmm::copy(md.complex_variable(A), md.set_complex_variable(A0));
+ } else {
+ gmm::copy(md.real_variable(U), md.set_real_variable(U0));
+ gmm::copy(md.real_variable(V), md.set_real_variable(V0));
+ gmm::copy(md.real_variable(A), md.set_real_variable(A0));
+ }
+ }
+
+ Hilber_Hughes_Taylor_scheme(model &md, std::string varname,
+ scalar_type al, scalar_type be,
+ scalar_type ga) {
+ U = varname;
+ U0 = "Previous_" + U;
+ V = "Dot_" + U;
+ V0 = "Previous_Dot_" + U;
+ A = "Dot2_" + U;
+ A0 = "Previous_Dot2_" + U;
+ alpha = al; beta = be; gamma = ga;
+ GMM_ASSERT1(alpha >= -scalar_type(1) && alpha <= scalar_type(0)
+ && beta > scalar_type(0) && beta <= scalar_type(1)
+ && gamma >= scalar_type(0.5) && gamma <= scalar_type(1),
+ "Invalid parameter values for the Hilber Hughes Taylor
scheme");
+
+ if (!(md.variable_exists(V)))
+ md.add_affine_dependent_variable(V, U);
+ if (!(md.variable_exists(A)))
+ md.add_affine_dependent_variable(A, U);
+
+ const mesh_fem *mf = md.pmesh_fem_of_variable(U);
+ size_type s = md.is_complex() ? gmm::vect_size(md.complex_variable(U))
+ : gmm::vect_size(md.real_variable(U));
+
+ if (mf) {
+ if (!(md.variable_exists(U0))) md.add_fem_data(U0, *mf);
+ if (!(md.variable_exists(V0))) md.add_fem_data(V0, *mf);
+ if (!(md.variable_exists(A0))) md.add_fem_data(A0, *mf);
+ } else {
+ if (!(md.variable_exists(U0))) md.add_fixed_size_data(U0, s);
+ if (!(md.variable_exists(V0))) md.add_fixed_size_data(V0, s);
+ if (!(md.variable_exists(A0))) md.add_fixed_size_data(A0, s);
+ }
+ }
+
+ };
+
+ void add_Hilber_Hughes_Taylor_scheme(model &md, const std::string &varname,
+ scalar_type alpha, scalar_type beta,
+ scalar_type gamma) {
+ ptime_scheme ptsc = std::make_shared<Hilber_Hughes_Taylor_scheme>
+ (md, varname, alpha, beta, gamma);
+ md.add_time_integration_scheme(varname, ptsc);
+ }