Time-stepping propagation in C++ wave problem

Question

Hello,

I want to propagate a wave in a square mesh going from left to right.

What I got so far is the left boundary changing according the the sinusoidal law I imposed, but the solution is not propagating. Clearly I must be doing something wrong in the time-stepping, like I'm probably not saving the current solution to be reused as initial condition.

I have no idea how to do so. Any help is appreciated.

---

#include <dolfin.h>
#include "Wave.h"
#include <mshr.h>

using namespace dolfin;

// Source term for boundary condition
class Source_BC : public Expression
{
    public:

    // Constructor
    Source_BC() : t(0) {}

    // Evaluate wave at left boundary
  void eval(Array<double>& values, const Array<double>& x) const
  {
    double eps = 1e-1;
    values[0] = -std::cos(eps*t);
  }

  // Current time
  double t;
};

// Sub domain for Dirichlet boundary condition
class DirichletBCLeft : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[0] < -1.0 + DOLFIN_EPS;
  }
};

class DirichletBCUp : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[1] > 1.0 - DOLFIN_EPS ;
  }
};

class DirichletBCRight : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[0] > 1.0 - DOLFIN_EPS;
  }
};

class DirichletBCDown : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[1] < -1.0 + DOLFIN_EPS;
  }
};


int main()
{
    // Parameters for the Newmark beta method
    double beta = 0.25;
    double gamma = 0.5;

    auto dt = std::make_shared<Constant>(1); 
  double T = 10*(*dt);

    // Mesh definition
    auto rectangle = mshr::Rectangle(dolfin::Point(-1.0, 1.0, 0.0), dolfin::Point(1.0, -1.0, 0.0));
    auto mesh = mshr::generate_mesh(rectangle, 50);

    auto V = std::make_shared<Wave::FunctionSpace>(mesh);

    auto b1 = std::make_shared<Source_BC>();
  auto boundary1 = std::make_shared<DirichletBCLeft>();
  auto boundary2 = std::make_shared<DirichletBCUp>();
  auto boundary3 = std::make_shared<DirichletBCRight>();
  auto boundary4 = std::make_shared<DirichletBCDown>();
  DirichletBC bc1(V, b1, boundary1);
  DirichletBC bc2(V, std::make_shared<Constant>(0.0), boundary2);
  DirichletBC bc3(V, std::make_shared<Constant>(0.0), boundary3);
  DirichletBC bc4(V, std::make_shared<Constant>(0.0), boundary4);

  std::vector<const DirichletBC*> bcs = {&bc1, &bc2, &bc3, &bc4};

  Function u(V);
  Function u0(V);

  auto f = std::make_shared<Constant>(0.0);

  // Create forms
  Wave::BilinearForm a(V, V);
  Wave::LinearForm L(V);

  L.f = f;

  //solve(a == L, u, bcs);

  File file("Wave.pvd");

  double t;
  while (t < T + DOLFIN_EPS)
  {
    b1->t = t;
    solve(a == L, u, bcs);
    *u0.vector() = *u.vector();
    t += *dt;
    file << u;
    }

  // Plot solution
  plot(u0);
  interactive();

  return 0;
}

---

# Compile this form with FFC: ffc -l dolfin Wave.ufl.

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)

a = u*v*dx
L = f*v*dx

---

Animaton of the results so far:

Comments

Can you link to a document which states the equations you're trying to solve, and possibly even their variational formulation? Examining your ufl file, you just have a mass matrix and load vector. I see no convection term.

---

Thank for your comment, that's something I haven't yet started to think about.
What I'm trying to do is to replicate the wave behaviour seen here, but without the obstacle with different density inside, so I didn't think about any equations yet.
The fenics file they used, however, seems to be very old and not up to date anymore.

I have found a working propagating wave example in Python and I think I might be able to "translate" it in C++, even if the right boundary makes the wave "bounce", and I have no idea yet how to solve that.

---

# From http://fenicsproject.org/qa/3441/time-evolution-of-the-scalar-wave-equation
from dolfin import *
#import plot

mesh = RectangleMesh(Point(0,1), Point(1, 0), 30, 30, "right/left")
Va = FunctionSpace(mesh, "Lagrange", 1)
Vb = FunctionSpace(mesh, "Lagrange", 1)
V  = Va*Vb

t  = 0.0
dt = 0.01

u0 = Expression(('x[0]', 'cos(10*t)'),t=0)
u0.t = t

def boundary(x, on_boundary):
  return x[0] > 1 - DOLFIN_EPS #or x[1] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS

def boundary2(x, on_boundary):
  return x[0] < DOLFIN_EPS and on_boundary

#dbc1 = DirichletBC(V, Constant(('0.0', '0.0')), boundary)
dbc2 = DirichletBC(V, u0, boundary2)

#dbc = [dbc1, dbc2]

class VolumeInitialCondition(Expression):
  def eval(self, value, x):
    value[0] = 0
    value[1] = 0

  def value_shape(self):
    return (2,)

(du, u)           = TrialFunctions(V)
(test_du, test_u) = TestFunctions(V)

w_prev = interpolate(VolumeInitialCondition(), V)
du_prev, u_prev = w_prev.split()

lhs1 = (test_u*u + test_du*du)*dx
lhs2 = dt*inner(grad(u), grad(test_du))*dx - dt*du*test_u*dx
rhs1 = test_u*u_prev*dx
rhs2 = test_du*du_prev*dx

A = assemble(lhs1 + lhs2)
b = assemble(rhs1 + rhs2)
dbc2.apply(A)
dbc2.apply(b)
#dbc1.apply(A)
#dbc1.apply(b)

w = Function(V)
solve(A, w.vector(), b)
output_file = File("scalar_results_euler_back.pvd", "compressed")
wave = w.split()[1]
wave.rename("WaveFunction", wave.name())
output_file << (wave, t)

for step in range(300):
  # store previous values, increment time step
    u0.t = t
    dut, ut = w.split(deepcopy = True)
    du_prev.assign(dut)
    u_prev.assign(ut)

  # assemble forms and solve
    A = assemble(lhs1 + lhs2)
    b = assemble(rhs1 + rhs2)
    dbc2.apply(A)
    dbc2.apply(b)
    #dbc1.apply(A)
    #dbc1.apply(b)
    solve(A, w.vector(), b)

  # save information
    wave = w.split()[1]
    wave.rename("WaveFunction", wave.name())
    output_file << (wave, t)
    t += dt

---

I will update the post as soon as I have some results.

---

So, I have been trying and trying to translate the above Python code in C++ but I have some problems with the code.
Namely, I am in troubles with the definition of the functions v and v_prec, representing new and last solution respectively.

Before entering the time loop, I want the first solution to be the initial solution, so

auto v = std::make_shared<Function>(V);
auto v_prev = std::make_shared<InitialConditions>();

where InitialConditions is an expression in a class (see full code below), so that I can impose

L.v_prev = v_prev;

Problems arise inside the time loop when i try to update v_prev with

v_prev = v;

or

*v_prev = *v;

or

v_prev.vector() = v.vector();

and so on, since

error: no match for ‘operator=’ (operand types are ‘InitialConditions’ and ‘dolfin::Function’)
   *v_prev = *v;

If i however had written

Function v_prev(V);
InitialConditions In;
v_prev.interpolate(In);

then I wouldn't be able to set the coefficient v_prev in the linear form L since

error: no match for ‘operator=’ (operand types are ‘dolfin::CoefficientAssigner’ and ‘dolfin::Function’)
  L.v_prev = v_prev;

So what can I do to overcome this problem?

Full code

---

#ffc -l dolfin Wave1.ufl

# Define function spaces
P1 = FiniteElement("Lagrange", triangle, 1)

V = P1*P1

# Define trial and test functions
(u, du)  = TrialFunction(V)
(test_u, test_du) = TestFunction(V)

# Define coefficients
k = Constant(triangle)
v_prev = Coefficient(V)

(u_prev, du_prev) = split(v_prev)

# Define bilinear and linear forms
a = inner(test_u,u)*dx + inner(test_du,du)*dx + k*inner(grad(u), grad(test_du))*dx - k*du*test_u*dx
L = test_u*u_prev*dx + test_du*du_prev*dx

---

#include <dolfin.h>
#include "Wave1.h"
#include <mshr.h>

using namespace dolfin;

// Source term for boundary condition
class Source_BC : public Expression
{
    public:

    // Constructor
    Source_BC() : Expression(2) {}

    // Evaluate wave at left boundary
  void eval(Array<double>& values, const Array<double>& x) const
  {
    double eps = 5e-2;
    values[0] = std::cos(eps*t);
    values[1] = 0.0;
  }

  // Current time
  double t;
};

class InitialConditions : public Expression
{
public:

    InitialConditions() : Expression(2) {}

  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0]= 0.0;
    values[1]= 0.0;
  }
};


// Sub domain for Dirichlet boundary condition
class DirichletBCLeft : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[0] < -1.0 + DOLFIN_EPS;
  }
};

class DirichletBCRight : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[0] > 1.0 - DOLFIN_EPS;
  }
};

int main()
{
    double dt = 0.05; 
  double T = 10*(dt);

  auto rectangle = mshr::Rectangle(dolfin::Point(-1.0, 1.0, 0.0), dolfin::Point(1.0, -1.0, 0.0));
    auto mesh = mshr::generate_mesh(rectangle, 20);

    auto V = std::make_shared<Wave1::FunctionSpace>(mesh);
    auto P1 = std::make_shared<Wave1::FunctionSpace>(mesh);

    auto bound1 = std::make_shared<Source_BC>();
  auto boundary1 = std::make_shared<DirichletBCLeft>();
  auto boundary3 = std::make_shared<DirichletBCRight>();

  DirichletBC bc1(V, bound1, boundary1);
  DirichletBC bc3(V, std::make_shared<const dolfin::Constant>(0.0, 0.0), boundary3);
  std::vector<const DirichletBC*> bcs = {&bc1, &bc3};
  //bcs.push_back(&bc1); bcs.push_back(&bc3);


  auto v = std::make_shared<Function>(V);
  auto v_prev = std::make_shared<InitialConditions>(); 

  //Function v_prev(V);
  //InitialConditions In;
  //v_prev.interpolate(In);


  auto k = std::make_shared<Constant>(dt);

  Wave1::BilinearForm a(V, V);
  Wave1::LinearForm L(V);

    a.k = k;
    L.v_prev = v_prev;

    File file1("results/wave.pvd");

    double t = dt;

    while (t < T + DOLFIN_EPS)
    {
        bound1->t = t;

        solve(a == L, *v, bcs);

        *v_prev = *v;

        //Function& u = v[0];
        //Function& du = v[1];

        file1 << *v;

        t += dt;
        std::cout << "t = " << t << "\n";
    }

    //plot(*v, "Wave1");
    //interactive();

    return 0;
}