From 697797fe6ea447e94164ebc589ab00f59aa5473f Mon Sep 17 00:00:00 2001 From: Pile Raphael Date: Mon, 20 Feb 2023 17:48:05 +0100 Subject: [PATCH 1/3] Keep parameters data configurable at the top of the file --- examples/time_of_flight/fedm-tof.py | 17 +++++++++++------ .../integrated_tests/time_of_flight/fedm_tof.py | 16 ++++++++++------ 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/examples/time_of_flight/fedm-tof.py b/examples/time_of_flight/fedm-tof.py index 9df3f7b..8fe614e 100644 --- a/examples/time_of_flight/fedm-tof.py +++ b/examples/time_of_flight/fedm-tof.py @@ -42,6 +42,11 @@ M = me charge = -elementary_charge equation_type = ['drift-diffusion-reaction'] # Defining the type of the equation (reaction | diffusion-reaction | drift-diffusion-reaction) +wez = 1.7e5 # electron velocity +De = 0.12 # electron diffusion coefficient +alpha_e = 5009.51 # reaction coefficient + + log('properties', files.model_log, gas, model, particle_species_type, M, charge) # Writting particle properties into a log file vtkfile_u = output_files('pvd', 'number density', particle_species_type) # Setting-up output files @@ -99,20 +104,20 @@ 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[1]-w*t, 2)+pow(x[0], 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_analytical = Expression('std::log(exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)/pow(4*D*t*pi,1.5))', D = De, w = wez, alpha = alpha_e, 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 = interpolate(Constant(('0','1.7e5')), W) # Electron drift velocity [m/s] -D = interpolate(Constant(0.12), V) # Diffusion coefficient [m^2/s] -alpha_eff = interpolate(Constant(5009.51), V) #Effective ionization coefficient [1/m] +w = interpolate(Constant(('0', wez)), W) # Electron drift velocity [m/s] +D = interpolate(Constant(De), V) # Diffusion coefficient [m^2/s] +alpha_eff = interpolate(Constant(alpha_e), V) #Effective ionization coefficient [1/m] Gamma = -grad(D*exp(u)) + w*exp(u) # Defining electron flux [m^{-2} s^{-1}] -f = Expression('exp(-(pow(x[1]-w*t, 2)+pow(x[0], 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 +f = Expression('exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)*(w*alpha)/(8*pow(pi,1.5)*pow(D*t, 1.5))', D = De, w = wez, alpha = alpha_e, t = t, pi=pi, degree = 2) # Defining source term F = weak_form_balance_equation_log_representation(equation_type[0], dt, dt_old, dx, u, u_old, u_old1, v, f, Gamma, r) # Definition of variational formulation of the balance equation for the electrons -u_new.assign(interpolate(Expression('std::log(exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)/pow(4.0*D*t*pi,1.5) + DOLFIN_EPS)', D = 0.12, w = 1.7e5, alpha = 5009.51, t = t, pi=pi, degree = 2), V)) # Setting up initial guess for nonlinear solver +u_new.assign(interpolate(Expression('std::log(exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)/pow(4.0*D*t*pi,1.5) + DOLFIN_EPS)', D = De, w = wez, alpha = alpha_e, t = t, pi=pi, degree = 2), V)) # Setting up initial guess for nonlinear solver # ============================================================================ # Setting-up nonlinear solver diff --git a/tests/integrated_tests/time_of_flight/fedm_tof.py b/tests/integrated_tests/time_of_flight/fedm_tof.py index 2b64e55..35d52eb 100644 --- a/tests/integrated_tests/time_of_flight/fedm_tof.py +++ b/tests/integrated_tests/time_of_flight/fedm_tof.py @@ -47,6 +47,10 @@ def main(output_dir = None): M = me charge = -elementary_charge equation_type = ['drift-diffusion-reaction'] # Defining the type of the equation (reaction | diffusion-reaction | drift-diffusion-reaction) + wez = 1.7e5 # Electron drift velocity z component [m/s] + De = 0.12 # Electron Diffusion coefficient [m^2/s] + alpha_e = 5009.51 # Effective ionization coefficient [1/m] + log('properties', files.model_log, gas, model, particle_species_type, M, charge) # Writting particle properties into a log file vtkfile_u = output_files('pvd', 'number density', particle_species_type) # Setting-up output files @@ -109,20 +113,20 @@ def main(output_dir = None): 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[1]-w*t, 2)+pow(x[0], 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_analytical = Expression('std::log(exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)/pow(4*D*t*pi,1.5))', D = De, w = wez, alpha = alpha_e, 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 = interpolate(Constant(('0','1.7e5')), W) # Electron drift velocity [m/s] - D = interpolate(Constant(0.12), V) # Diffusion coefficient [m^2/s] - alpha_eff = interpolate(Constant(5009.51), V) #Effective ionization coefficient [1/m] + w = interpolate(Constant(('0', wez)), W) # Electron drift velocity [m/s] + D = interpolate(Constant(De), V) # Diffusion coefficient [m^2/s] + alpha_eff = interpolate(Constant(alpha_e), V) # Effective ionization coefficient [1/m] Gamma = -grad(D*exp(u)) + w*exp(u) # Defining electron flux [m^{-2} s^{-1}] - f = Expression('exp(-(pow(x[1]-w*t, 2)+pow(x[0], 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 + f = Expression('exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)*(w*alpha)/(8*pow(pi,1.5)*pow(D*t, 1.5))', D = De, w = wez, alpha = alpha_e, t = t, pi=pi, degree = 2) # Defining source term F = weak_form_balance_equation_log_representation(equation_type[0], dt, dt_old, dx, u, u_old, u_old1, v, f, Gamma, r) # Definition of variational formulation of the balance equation for the electrons - u_new.assign(interpolate(Expression('std::log(exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)/pow(4.0*D*t*pi,1.5) + DOLFIN_EPS)', D = 0.12, w = 1.7e5, alpha = 5009.51, t = t, pi=pi, degree = 2), V)) # Setting up initial guess for nonlinear solver + u_new.assign(interpolate(Expression('std::log(exp(-(pow(x[1]-w*t, 2)+pow(x[0], 2))/(4.0*D*t)+alpha*w*t)/pow(4.0*D*t*pi,1.5) + DOLFIN_EPS)', D = De, w = wez, alpha = alpha_e, t = t, pi=pi, degree = 2), V)) # Setting up initial guess for nonlinear solver # ============================================================================ # Setting-up nonlinear solver From 1783a58259c1c6b2b0a5e3e29a626e01f3bec543 Mon Sep 17 00:00:00 2001 From: Pile Raphael Date: Mon, 20 Feb 2023 18:09:55 +0100 Subject: [PATCH 2/3] fix: Enable to use animation in paraview. Use assign when writing data (instead of recreating the function each time) so it keeps the same field_name in every vtu files. --- examples/time_of_flight/fedm-tof.py | 7 +++++-- tests/integrated_tests/time_of_flight/fedm_tof.py | 7 +++++-- .../integrated_tests/time_of_flight/test_time_of_flight.py | 2 +- 3 files changed, 11 insertions(+), 5 deletions(-) diff --git a/examples/time_of_flight/fedm-tof.py b/examples/time_of_flight/fedm-tof.py index 8fe614e..b9e1a74 100644 --- a/examples/time_of_flight/fedm-tof.py +++ b/examples/time_of_flight/fedm-tof.py @@ -134,6 +134,9 @@ nonlinear_solver.parameters['maximum_iterations'] = maximum_iterations # Setting up maximum number of iterations # nonlinear_solver.parameters["preconditioner"]="hypre_amg" # Setting the preconditioner, uncomment if iterative solver is used +n_exact = Function(V) # Defining function for storing the data (analytical solution) +n_num = Function(V) # Defining function for storing the data (numerical solution) + while abs(t-T_final)/T_final > 1e-6: t_old = t # Updating old time steps u_old1.assign(u_old) # Updating variable value in k-2 time step @@ -149,8 +152,8 @@ nonlinear_solver.solve(problem, u_new.vector()) # Solving the system of equations if abs(t-t_output)/t_output <= 1e-6: - n_exact = project(exp(u_analytical), V, solver_type='mumps') - n_num = project(exp(u_new), V, solver_type='mumps') + n_exact.assign(project(exp(u_analytical), V, solver_type='mumps')) + n_num.assign(project(exp(u_new), V, solver_type='mumps')) relative_error = errornorm(n_num, n_exact, 'l2')/norm(n_exact, 'l2') # Calculating relative difference between exact analytical and numerical solution with open(files.error_file, "a") as f_err: f_err.write('h_max = ' + str(h) + '\t dt = ' + str(dt.time_step) + '\t relative_error = ' + str(relative_error) + '\n') diff --git a/tests/integrated_tests/time_of_flight/fedm_tof.py b/tests/integrated_tests/time_of_flight/fedm_tof.py index 35d52eb..8439708 100644 --- a/tests/integrated_tests/time_of_flight/fedm_tof.py +++ b/tests/integrated_tests/time_of_flight/fedm_tof.py @@ -143,6 +143,9 @@ def main(output_dir = None): nonlinear_solver.parameters['maximum_iterations'] = maximum_iterations # Setting up maximum number of iterations # nonlinear_solver.parameters["preconditioner"]="hypre_amg" # Setting the preconditioner, uncomment if iterative solver is used + n_exact = Function(V) # Defining function for storing the data (analytical solution) + n_num = Function(V) # Defining function for storing the data (numerical solution) + while abs(t-T_final)/T_final > 1e-6: t_old = t # Updating old time steps u_old1.assign(u_old) # Updating variable value in k-2 time step @@ -158,8 +161,8 @@ def main(output_dir = None): nonlinear_solver.solve(problem, u_new.vector()) # Solving the system of equations if abs(t-t_output)/t_output <= 1e-6: - n_exact = project(exp(u_analytical), V, solver_type='mumps') - n_num = project(exp(u_new), V, solver_type='mumps') + n_exact.assign(project(exp(u_analytical), V, solver_type='mumps')) + n_num.assign(project(exp(u_new), V, solver_type='mumps')) relative_error = errornorm(n_num, n_exact, 'l2')/norm(n_exact, 'l2') # Calculating relative difference between exact analytical and numerical solution with open(files.error_file, "a") as f_err: f_err.write('h_max = ' + str(h) + '\t dt = ' + str(dt.time_step) + '\t relative_error = ' + str(relative_error) + '\n') diff --git a/tests/integrated_tests/time_of_flight/test_time_of_flight.py b/tests/integrated_tests/time_of_flight/test_time_of_flight.py index 28ce35b..78656d5 100644 --- a/tests/integrated_tests/time_of_flight/test_time_of_flight.py +++ b/tests/integrated_tests/time_of_flight/test_time_of_flight.py @@ -34,7 +34,7 @@ def data_dir(tmpdir_factory): @pytest.fixture def tof_data(data_dir): path = data_dir / "number density" / "electrons" / "electrons000000.vtu" - return read_vtu(path, field_name="f_3199") + return read_vtu(path, field_name="f_52") @pytest.fixture From aa0a11044f8db9ae3b7d962ce629808134f1f6ce Mon Sep 17 00:00:00 2001 From: Pile Raphael <23015350+RaphaelPile@users.noreply.github.com> Date: Wed, 22 Feb 2023 15:36:57 +0100 Subject: [PATCH 3/3] feat: add a 1D version of tof example --- examples/time_of_flight_1D/fedm-tof_1d.py | 176 ++++++++++++++++++++++ 1 file changed, 176 insertions(+) create mode 100644 examples/time_of_flight_1D/fedm-tof_1d.py diff --git a/examples/time_of_flight_1D/fedm-tof_1d.py b/examples/time_of_flight_1D/fedm-tof_1d.py new file mode 100644 index 0000000..6f83e14 --- /dev/null +++ b/examples/time_of_flight_1D/fedm-tof_1d.py @@ -0,0 +1,176 @@ +""" +This small example presents modelling of time-of-flight experiment and is used for verification purpose. +Namely, for this particular test case the exact analytical solution for the electron number density exists, +so it is used to verify the accuracy of the code using the method of exact solutions. L2 norm is used to +quantify the difference between the solutions and to determine the mesh and time order-of-accuracy of the +code for solving the balance equation. +""" +from dolfin import * +import numpy as np +from timeit import default_timer as timer +import time +import sys +from fedm.physical_constants import * +from fedm.file_io import * +from fedm.functions import * + +# Optimization parameters. +parameters["form_compiler"]["optimize"] = True +parameters["form_compiler"]["cpp_optimize"] = True +parameters['krylov_solver']['nonzero_initial_guess'] = True + +# Defining tye of used solver and its parameters. +linear_solver = "mumps" # Setting linear solver: mumps | gmres | lu +maximum_iterations = 50 # Setting up maximum number of nonlinear solver iterations +relative_tolerance = 1e-10 # Setting up relative tolerance + +# ============================================================================ +# Definition of the simulation conditions, model and coordinates +# ============================================================================ +model = 'Time_of_flight' # Model name +coordinates = 'cylindrical' # Coordinates choice +gas = 'Air' +Tgas = 300.0 # Gas temperature in [K] +p0 = 760.0 # Pressure in [Torr] +N0 = p0*3.21877e22 # Number density in [m^-3] + + +# ============================================================================ +# Defining number of species for which the problem is solved, their properties and creating output files. +# ============================================================================ +number_of_species = 1 +particle_species_type = ['electrons', 'analytical solution'] # Defining particle species types, where the values are used as names for output files +M = me +charge = -elementary_charge +equation_type = ['drift-diffusion-reaction'] # Defining the type of the equation (reaction | diffusion-reaction | drift-diffusion-reaction) +wez = 1.7e5 # Electron drift velocity z component [m/s] +De = 0.12 # Electron Diffusion coefficient [m^2/s] +alpha_e = 5009.51 # Effective ionization coefficient [1/m] +x0 =3e-4 +l = 0.00004 # 0.04 # Gaussian characteristic width + +log('properties', files.model_log, gas, model, particle_species_type, M, charge) # Writting particle properties into a log file +vtkfile_u = output_files('pvd', 'number density', particle_species_type) # Setting-up output files + +# ============================================================================ +# Definition of the time variables. +# ============================================================================ +t_old = None # Previous time step +t0 = 0 # Initial time step +t = t0 # Current time step +T_final = 3e-9 # Simulation time end [s] + +dt_init = 1e-11 # Initial time step size +dt = Expression("time_step", time_step = dt_init, degree = 0) # Time step size [s] +dt_old = Expression("time_step", time_step = 1e30, degree = 0) # Time step size [s] set up as a large value to reduce the order of BDF. + +t_output_step =10* dt_init #0.5e-9 # Time step intervals at which the results are printed to file +t_output = t0 + 10*dt_init # Initial output time step + + +# ============================================================================ +# Defining the geometry of the problem and corresponding boundaries. +# ============================================================================ +if coordinates == 'cylindrical': + #r = Expression('x[0]', degree = 1) + z = Expression('x[0]', 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(4000, 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] - x0 -w*t)/l, 2))/(1 + 4.0*D*t/pow(l,2))+alpha*w*t)/pow(1 + 4*D*t/pow(l,2),0.5))', x0=x0, D = De, w = wez, alpha = alpha_e, t = t, pi=pi, l=l, 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 = interpolate(Constant((wez, )), W) +D = interpolate(Constant(De), V) # Diffusion coefficient [m^2/s] +alpha_eff = interpolate(Constant(alpha_e), V) #Effective ionization coefficient [1/m] + +Gamma = -grad(D*exp(u)) + w*exp(u) # Defining electron flux [m^{-2} s^{-1}] +f = Expression('exp(-(pow((x[0] - x0 -w*t)/l, 2))/(1 + 4.0*D*t/pow(l,2))+alpha*w*t)*(w*alpha)/pow(1 + 4*D*t/pow(l,2),0.5)', x0=x0, D = De, w = wez, alpha = alpha_e, t = t, pi=pi, l=l, degree = 2) # Defining source term + +F = weak_form_balance_equation_log_representation(equation_type[0], dt, dt_old, dx, u, u_old, u_old1, v, f, Gamma) # Definition of variational formulation of the balance equation for the electrons + +#u_new.assign(interpolate(Expression('std::log(exp(-(pow((x[0] - x0)/l, 2))) + DOLFIN_EPS)', x0=x0, l=l, degree = 2), V)) # Setting up initial guess for nonlinear solver +u_new.assign(interpolate(Expression('std::log(exp(-(pow((x[0] - x0 -w*t)/l, 2))/(1 + 4.0*D*t/pow(l,2))+alpha*w*t)/pow(1 + 4*D*t/pow(l,2),0.5) + DOLFIN_EPS)', x0=x0, D = De, w = wez, alpha = alpha_e, t = t, pi=pi, l=l, degree = 3), V)) # Setting up initial guess for nonlinear solver + +# ============================================================================ +# Setting-up nonlinear solver +# ============================================================================ +# Defining the problem +F = action(F, u_new) +J = derivative(F, u_new, u) +problem = Problem(J, F, []) + +# Initializing nonlinear solver and setting up the parameters +nonlinear_solver = PETScSNESSolver() # Nonlinear solver initialization +nonlinear_solver.parameters['relative_tolerance'] = relative_tolerance # Setting up relative tolerance of the nonlinear solver +nonlinear_solver.parameters["linear_solver"]= linear_solver # Setting up linear solver +nonlinear_solver.parameters['maximum_iterations'] = maximum_iterations # Setting up maximum number of iterations +# nonlinear_solver.parameters["preconditioner"]="hypre_amg" # Setting the preconditioner, uncomment if iterative solver is used + +n_exact = Function(V) # Defining function for storing the data at k-2 time step +n_num = Function(V) # Defining function for storing the data at k time step + +while abs(t-T_final)/T_final > 1e-6: + t_old = t # Updating old time steps + u_old1.assign(u_old) # Updating variable value in k-2 time step + u_old.assign(u_new) # Updating variable value in k-1 time step + t += dt.time_step # Updating the new time steps + + log('time', files.model_log, t) # Time logging + print_time(t) # Printing out current time step + + f.t=t # Updating the source term for the current time step + u_analytical.t = t # Updating the analytical solution for the current time step + + nonlinear_solver.solve(problem, u_new.vector()) # Solving the system of equations + + if abs(t-t_output)/t_output <= 1e-6: + + n_exact.assign(project(exp(u_analytical), V, solver_type='mumps')) + n_num.assign(project(exp(u_new), V, solver_type='mumps')) + + relative_error = errornorm(n_num, n_exact, 'l2')/norm(n_exact, 'l2') # Calculating relative difference between exact analytical and numerical solution + with open(files.error_file, "a") as f_err: + f_err.write('h_max = ' + str(h) + '\t dt = ' + str(dt.time_step) + '\t relative_error = ' + str(relative_error) + '\n') + if(MPI.rank(MPI.comm_world)==0): + print(relative_error) + vtkfile_u[0] << (n_num, t) + vtkfile_u[1] << (n_exact, t) + t_output += t_output_step + print("Saved") + + if t > (t0 + dt_init): + dt_old.time_step = dt.time_step # For initialization BDF1 is used and after initial step BDF2 is activated + if(MPI.rank(MPI.comm_world)==0): + print(str(dt_old.time_step) + '\t' + str(dt.time_step) +'\n') + + +print("Finished") \ No newline at end of file