1D Time of flight

[RaphaelPile]

Hello there !

I am working on 1D version of time of flight example. I have not modified the beginning and the end of the code, just the boundaries and the mesh at first. I got a first problem with the velocity vector, which wa solved by using w_vec = Expression(('1.7e5', 0.), element=W.ufl_element()) instead of w = interpolate(Constant(('0','1.7e5')), W) in the flux Gamma calculation.

But now I am a bit stuck with the error I get below. Maybe my previous modification was a mistake, I don't really know....

Here is my code (the modified part of time-of-flight.py) and the error log.

Any help is welcome.

# ============================================================================
# Defining the geometry of the problem and corresponding boundaries.
# ============================================================================
if coordinates == 'cylindrical':
    r = Expression('x[0]', degree = 1)
    z = Expression('x[1]', degree = 1)

# Gap length and width
box_height = 1e-3 # [m]

boundaries = [['point', 0.0, 0.0],\
                ['point', 0.0, box_height]] # Defining list of boundaries lying between z0 and z1

# ============================================================================
# Mesh setup. Structured mesh is generated using built-in mesh generator.
# ============================================================================
mesh = IntervalMesh(80, 0. , box_height) # Generating structured mesh.
mesh_statistics(mesh) # Prints number of elements, minimal and maximal cell diameter.
h = MPI.max(MPI.comm_world, mesh.hmax()) # Maximuml cell size in mesh.

log('conditions', files.model_log, dt.time_step, 'None', p0, box_height, N0, Tgas)
log('initial time', files.model_log, t)

# ============================================================================
# Defining type of elements and function space, test functions, trial functions and functions for storing variables, and weak form
# of equation.
# ============================================================================
V = FunctionSpace(mesh, 'P', 2) # Defining function space
W = VectorFunctionSpace(mesh, 'P', 1) # Defining vector function space

u = TrialFunction(V) # Defining trial function
v = TestFunction(V) # Defining test function
u_old = Function(V) # Defining function for storing the data at k-1 time step
u_old1 = Function(V) # Defining function for storing the data at k-2 time step
u_new = Function(V) # Defining function for storing the data at k time step

u_analytical  = Expression('std::log(exp(-(pow(x[0]-w*t, 2))/(4.0*D*t)+alpha*w*t)/pow(4*D*t*pi,1.5))', D = 0.12, w = 1.7e5, alpha = 5009.51, t = t, pi=pi,  degree = 3) # Analytical solution of the particle balance equation.
u_old.assign(interpolate(u_analytical , V)) # Setting up value at k-1 time step
u_old1.assign(interpolate(u_analytical , V)) # Setting up value at k-2 time step

w_vec = Expression(('1.7e5', 0.), element=W.ufl_element())
D = interpolate(Constant(0.12), V) # Diffusion coefficient [m^2/s]
alpha_eff = interpolate(Constant(5009.51), V) #Effective ionization coefficient [1/m]

Gamma = -grad(D*exp(u)) + w_vec*exp(u) # Defining electron flux [m^{-2} s^{-1}]

Traceback (most recent call last):
File "fedm-tof_1d.py", line 143, in
nonlinear_solver.solve(problem, u_new.vector()) # Solving the system of equations
File "/home/fenics/.local/lib/python3.6/site-packages/fedm/functions.py", line 200, in J
df.assemble(self.bilinear_form, tensor=A)
File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/assembling.py", line 213, in assemble
assembler.assemble(tensor, dolfin_form)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

---

*** [email protected]

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** Error: Unable to assemble form.
*** Reason: Invalid value dimension 0 for coefficient 6 (got 2 but expecting 1). You might have forgotten to specify the value dimension correctly in an Expression subclass.
*** Where: This error was encountered inside AssemblerBase.cpp.
*** Process: 0

---

*** DOLFIN version: 2019.1.0
*** Git changeset: 74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Segmentation fault

[RaphaelPile]

I got some help from @AleksandarJ1984 : solved removing "r" (radius) and using w_vec = interpolate(Constant((1.7e5, )), W). However, I don't get the expected analytical results.

By the way, when I have finished this part, I will try to do a PR asap since it is linked to #3

[RaphaelPile]

If someone has an explanation for the following source term:

f = Expression('exp(-(pow(x[0]-w*t, 2))/(4.0*D*t)+alpha*w*t)*(w*alpha)/(8*pow(pi,1.5)*pow(D*t, 1.5))', D = 0.12, w = 1.7e5, alpha = 5009.51, t = t, pi=pi, degree = 2) # Defining source term

The first part corresponds to the analytical expression of ne, multiplied by the velocity and alpha. So why it is not in logarithmic scale for the "ne" expression ?

[markus-m-becker]

In fact, the source term definition does not change whether the log formulation is used or not. With log formulation the time derivative term changes to n_e*dN_e/dt, where n_e = exp(N_e). The equation is solved for N_e = log(n_e) but the rest of the equation, including the source term, remains unchanged (using n_e and not N_e).

[markus-m-becker]

I got some help from @AleksandarJ1984 : solved removing "r" (radius) and using w_vec = interpolate(Constant((1.7e5, )), W). However, I don't get the expected analytical results.

Make sure that you consider the proper analytical solution to compare with when the problem is reduced to 1D.

[RaphaelPile]

The equation is solved for N_e = log(n_e) but the rest of the equation, including the source term, remains unchanged (using n_e and not N_e).
Thank you very much for this information. I understand better the Logarithmic stabilization method now.

Make sure that you consider the proper analytical solution to compare with when the problem is reduced to 1D.
Yes you are right, but I will be looking in this direction.

[RaphaelPile]

This is working now, I will do a PR as soon as possible.

[AleksandarJ1984]

@RaphaelPile Thank you very much for making the 1D example. I will check it out and merge it with other examples.