Fix timestep bug

[kaelyndunnell]

Proposed changes

Fixes bug in timestep implementation in accordance with issue #920.

Types of changes

What types of changes does your code introduce to FESTIM?

• Bugfix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Code refactoring
• Documentation Update (if none of the other choices apply)
• New tests

Checklist

• Black formatted
• Unit tests pass locally with my changes
• I have added tests that prove my fix is effective or that my feature works
• I have added necessary documentation (if appropriate)

Further comments

If this is a relatively large or complex change, kick off the discussion by explaining why you chose the solution you did and what alternatives you considered, etc...

[RemDelaporteMathurin]

Your test is wrong as you have the wrong expected value

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 97.49%. Comparing base (0769a6f) to head (933650d).
Report is 8 commits behind head on fenicsx.

Additional details and impacted files

@@           Coverage Diff            @@
##           fenicsx     #922   +/-   ##
========================================
  Coverage    97.49%   97.49%           
========================================
  Files           44       44           
  Lines         2193     2193           
========================================
  Hits          2138     2138           
  Misses          55       55           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[RemDelaporteMathurin]

Thanks! CI is passing.

It's failing on the nightly build of FEniCS but we'll wait a bit to see if that gets resolved.

Could you format your code please?

[jorgensd]

It’s failing because of the check added by @garth-wells in https://github.com/FEniCS/dolfinx/pull/3502/files
I am not sure why it shouldn’t be allowed to make a mixed element with a single element.

[garth-wells]

It’s failing because of the check added by @garth-wells in https://github.com/FEniCS/dolfinx/pull/3502/files I am not sure why it shouldn’t be allowed to make a mixed element with a single element.

Mixed elements do not have a well-defined value shape, whereas 'base' sub-elements of a mixed element do. An element should be what it is. If it's mixed it should be mixed, if it's not mixed it should not be mixed.

[jorgensd]

It’s failing because of the check added by @garth-wells in https://github.com/FEniCS/dolfinx/pull/3502/files I am not sure why it shouldn’t be allowed to make a mixed element with a single element.

Mixed elements do not have a well-defined value shape, whereas 'base' sub-elements of a mixed element do. An element should be what it is. If it's mixed it should be mixed, if it's not mixed it should not be mixed.

@garth-wells
I still don't see why you shouldn't be allowed to do mixed_element([single_element]).
This is a perfectly valid use-case for more generalized problems where one can have M number of species in a problem, where 1 species is a special case. It would be silly to have to have dual code paths for handling this in any application code.

[RemDelaporteMathurin]

@garth-wells I agree with @jorgensd here. For FESTIM, not being able to define this special case would mean we have to rewrite many parts of the code since mixed elements and elements behave differently on many things (enforcing dirichlet BCs, post-processing, etc.). I think this is something we had to do in FESTIM1 and we were super happy to realise we could do mixed_element([element]) in dolfinx because it made the code so much cleaner!

[garth-wells]

@RemDelaporteMathurin can you share an example of where you rely on a single (base) element to behave like a mixed element?

[jorgensd]

@RemDelaporteMathurin can you share an example of where you rely on a single (base) element to behave like a mixed element?

let’s say you want to support having multiple diffusion equations co-located in space (but with different function spaces. Then you would have a mixed element ([el_0,el_1]). In other parts of the domain (with submesh involved, you might only have a single species. One doesn’t want to treat those differently.

[RemDelaporteMathurin]

@RemDelaporteMathurin can you share an example of where you rely on a single (base) element to behave like a mixed element?

@garth-wells we do not need a single (base) element to behave like a mixed element per se. It would just simplify the logic for us. For example, in the HydrogenTransportProblem class, this is how we assign functions (solutions, test functions, previous solutions, etc) to "species":

FESTIM/src/festim/hydrogen_transport_problem.py

Lines 485 to 511 in effcbac

	def assign_functions_to_species(self):
	"""Creates the solution, prev solution, test function and
	post-processing solution for each species, if model is multispecies,
	created a collapsed function space for each species"""
	if not self.multispecies:
	sub_solutions = [self.u]
	sub_prev_solution = [self.u_n]
	sub_test_functions = [ufl.TestFunction(self.function_space)]
	self.species[0].sub_function_space = self.function_space
	self.species[0].post_processing_solution = self.u
	else:
	sub_solutions = list(ufl.split(self.u))
	sub_prev_solution = list(ufl.split(self.u_n))
	sub_test_functions = list(ufl.TestFunctions(self.function_space))
	for idx, spe in enumerate(self.species):
	spe.sub_function_space = self.function_space.sub(idx)
	spe.post_processing_solution = self.u.sub(idx)
	spe.collapsed_function_space, _ = self.function_space.sub(
	idx
	).collapse()
	for idx, spe in enumerate(self.species):
	spe.solution = sub_solutions[idx]
	spe.prev_solution = sub_prev_solution[idx]
	spe.test_function = sub_test_functions[idx]

We have this check if self.multispecies: in multiple places to check if there are more than one element and if so, treat it as a mixed element.

This is how we do it in the HTransportProblemDiscontinuous class (using a mixed domain approach):

FESTIM/src/festim/hydrogen_transport_problem.py

Lines 990 to 1021 in effcbac

	element_CG = basix.ufl.element(
	basix.ElementFamily.P,
	subdomain.submesh.basix_cell(),
	degree,
	basix.LagrangeVariant.equispaced,
	)
	element = basix.ufl.mixed_element([element_CG] * nb_species)
	V = dolfinx.fem.functionspace(subdomain.submesh, element)
	u = dolfinx.fem.Function(V)
	u_n = dolfinx.fem.Function(V)
	# store attributes in the subdomain object
	subdomain.u = u
	subdomain.u_n = u_n
	# split the functions and assign the subfunctions to the species
	us = list(ufl.split(u))
	u_ns = list(ufl.split(u_n))
	vs = list(ufl.TestFunctions(V))
	for i, species in enumerate(unique_species):
	species.subdomain_to_solution[subdomain] = us[i]
	species.subdomain_to_prev_solution[subdomain] = u_ns[i]
	species.subdomain_to_test_function[subdomain] = vs[i]
	species.subdomain_to_function_space[subdomain] = V.sub(i)
	species.subdomain_to_post_processing_solution[subdomain] = u.sub(
	i
	).collapse()
	species.subdomain_to_collapsed_function_space[subdomain] = V.sub(
	i
	).collapse()
	name = f"{species.name}_{subdomain.id}"
	species.subdomain_to_post_processing_solution[subdomain].name = name

Here we can have mixed element with only one element and we don't have to have these checks everywhere.

This is an example for assigning functions to species but we have other places in the code where this happens, for instance assigning BCs:

FESTIM/src/festim/hydrogen_transport_problem.py

Lines 537 to 586 in effcbac

	def create_dirichletbc_form(self, bc):
	"""Creates a dirichlet boundary condition form
	Args:
	bc (festim.DirichletBC): the boundary condition
	Returns:
	dolfinx.fem.bcs.DirichletBC: A representation of
	the boundary condition for modifying linear systems.
	"""
	# create value_fenics
	if not self.multispecies:
	function_space_value = bc.species.sub_function_space
	else:
	function_space_value = bc.species.collapsed_function_space
	bc.create_value(
	temperature=self.temperature_fenics,
	function_space=function_space_value,
	t=self.t,
	)
	# get dofs
	if self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
	function_space_dofs = (
	bc.species.sub_function_space,
	bc.species.collapsed_function_space,
	)
	else:
	function_space_dofs = bc.species.sub_function_space
	bc_dofs = bc.define_surface_subdomain_dofs(
	facet_meshtags=self.facet_meshtags,
	function_space=function_space_dofs,
	)
	# create form
	if not self.multispecies and isinstance(bc.value_fenics, (fem.Function)):
	# no need to pass the functionspace since value_fenics is already a Function
	function_space_form = None
	else:
	function_space_form = bc.species.sub_function_space
	form = fem.dirichletbc(
	value=bc.value_fenics,
	dofs=bc_dofs,
	V=function_space_form,
	)
	return form

Assigning ICs:

FESTIM/src/festim/hydrogen_transport_problem.py

Lines 610 to 642 in effcbac

	def create_initial_conditions(self):
	"""For each initial condition, create the value_fenics and assign it to
	the previous solution of the condition's species"""
	if len(self.initial_conditions) > 0 and not self.settings.transient:
	raise ValueError(
	"Initial conditions can only be defined for transient simulations"
	)
	function_space_value = None
	for condition in self.initial_conditions:
	# create value_fenics for condition
	function_space_value = None
	if callable(condition.value):
	# if bc.value is a callable then need to provide a functionspace
	if not self.multispecies:
	function_space_value = condition.species.sub_function_space
	else:
	function_space_value = condition.species.collapsed_function_space
	condition.create_expr_fenics(
	mesh=self.mesh.mesh,
	temperature=self.temperature_fenics,
	function_space=function_space_value,
	)
	# assign to previous solution of species
	if not self.multispecies:
	condition.species.prev_solution.interpolate(condition.expr_fenics)
	else:
	idx = self.species.index(condition.species)
	self.u_n.sub(idx).interpolate(condition.expr_fenics)

Having everything as mixed elements would greatly simplify the logic on our side.