Is this a good way to implement the method of characteristics for advection equation?

Question

Dear everyone,

I tried to solve the advection equation
$$ u_t - c u_x = 0, \quad x \in[x_l, x_r] $$
with initial condition
$$ u0 = 0.5(u_l - u_r) \tanh(100(0.5 - x)) + u_r $$
with the method of characteristic.

I don't know how will you implement this method using FENICS. My approach is to define a function below

 class CharacteristicMethod(Expression):
    def eval(self, values, x):
        if ((x+ c*dt) > x_r):
            tmp = ur
        else:
            tmp = u(x+c*dt)
        values[0] = tmp
    def value_shape(self):
        return (1,)

then in every time step I use the interpolate method of the FEMfunction

u_init = InitialConditions(degree=1)
u.interpolate(u_init)

dt = 0.1
t = 0
Tend = 1
while (t < Tend):
    t += dt
    u_characteristic = CharacteristicMethod(degree=1)
    u.interpolate(u_characteristic)
    # I need to do some other things with u in every time step ...

However this approach is not very efficient. How will you implement the method of characteristic for the linear advection equation?

Thank you in advance!
Jesse

Comments

This has an explicit solution and so Fenics wouldn't be ncessary. Why do you want to use it?

---

Hi nabarnaf,

Thank you very much.
I have to solve another elliptic equation with the solution obtained after using method of characteristic, thus I want to use implement the method of characteristic in the Fenics environment.

Answer 1

Hi, I guess your problem is about speed, and that happens because dolfin generates c++ code, and when it needs to evaluate your Expression, it goes back to python (in every quadrature point!). To solve this, you can create the expression in c++ code. I have done it, and the increase is unbelievable. First define the expression:

    cpp_code = '''
    class test : public Expression
        {
        public:

        AdvSol() : Expression() { }

        void eval(Array<double>& values, const Array<double>& x) const {                
            double temp = 0;
            if (x[0] + c*dt >  x_r){
             values[0] = u_r;
             }
            else{
            Array<double> xcdt (1);
            xcdt[0] = x[0] + c*dt;
            u->eval(values, xcdt);
            }
        }

        std::shared_ptr<const Function> u; 
        double c;
        double dt;
        double u_r;
        double x_r;
    };
    '''

And then create the expression and initiate the requiered fields:

f = Expression(cpp_code, FiniteElement('CG', mesh.ufl_cell(), 2))
f.u = u0 #u0 must be a Function
f.c = 2.
f.dt = dt
f.u_r = u_r
f.x_r =c_r

Comments

Dear nabarnaf,

Thank you very much for your help. I didn't know that cpp code can be used in FENICS. Your code really helps me a lot.

Best wishes.

---

I just did something similar and thought it could help. If it works, you should close the question. Best regards!