Fenics don't work in few threads.

Hi. I use FENICS 1.0.0 on Windows. it's version have uBLAS or STL linear_algebra_backend only.

Whenever I add:

parameters["num_threads"] = 4;

to my code, I get this recurring on the terminal:

*** Warning: Form::coloring does not properly consider form type.

It also doesn't seem to be using all the cores for some reasons.

poisson.ufl

cell = tetrahedron
element = FiniteElement("Lagrange" , cell , 1)

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

d = Coefficient(element)
w = Coefficient(element*element*element)
S = Coefficient(element)
u0 = Coefficient(element)

dt1 = Constant(cell)
lam = Constant(cell)

a = (((dt1+lam)*u*v) + inner(w,nabla_grad(u))*v+ inner(d*nabla_grad(u),nabla_grad(v)))*dx
L = (S+dt1*u0)*v*dx

main.cpp

#include <dolfin.h>
#include "Poisson.h"
#include <iostream>
#include <windows.h>

using namespace dolfin;

class Source1 : public Expression
{
public:

  Source1() : Expression() {}

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

};
class Source2 : public Expression
{
public:

  Source2() : Expression() {}

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

};
class Source3 : public Expression
{
public:

  Source3(double& t) : Expression(), t(t) {}

  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = (x[0]*x[0]+x[1]*x[1]+(x[2]-0.5)*(x[2]-0.5)<1)?1/(0.5*t*t+1):0.0;
  }

private:

  double& t;

};
class Source4 : public Expression
{
public:

  Source4() : Expression(3) {}

  void eval(Array<double>& values, const Array<double>& x) const
  {
    values[0] = 10.0/5.0*x[2];
    values[1] = 0.0;
    values[2] = 0.1;
  }

};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary1 : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const 
  {
    return ((on_boundary)and((x[0]<=0.) or (x[2]>4.9)));
  }
};

int main()
{ parameters["num_threads"] = 4;

  // Create mesh and function space
  double dt = 1;
  double te = 10;
  double t = dt;

  Mesh mesh ("untitled.xml");

  Poisson::FunctionSpace V(mesh);
  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);
  Function u(V);
  Function u0(V);
  L.u0 = u0;
  // Define boundary condition

  Source2 d;
  a.d = d;
  Source3 S(t);
  L.S = S;
  Source4 w;
  a.w = w;

  Source1 g;
  DirichletBoundary1 boundary;
  DirichletBC bc1(V, g, boundary);


  Constant dt1(1/dt);
  a.dt1 = dt1;
  L.dt1 = dt1;
  Constant lam(0.01);
  a.lam = lam;

  // Define variational forms

  Matrix A;
  assemble (A, a);

  Vector b;

  File file("poisson.pvd");
  int start = GetTickCount();

  while (t <= te) 
{
  // Compute solution
  assemble (b,L);
  bc1.apply(A,b);
  solve(A, *u.vector(), b); 
  file << u;
  u0 = u;
  t = t + dt;  
}
  int end = GetTickCount();
  cout<<(end - start)/t*dt<<endl;

  plot(u); 

  return 0;
}

Comments (8)