festim-dev / FESTIM

Coupled hydrogen/tritium transport and heat transfer modelling using FEniCS
https://festim.readthedocs.io/en/stable/
Apache License 2.0
88 stars 22 forks source link

[BUG] Cannot export to VTX with multispecies #622

Closed RemDelaporteMathurin closed 10 months ago

RemDelaporteMathurin commented 10 months ago

Describe the bug When exporting to VTX with multispecies, it just exports zero everywhere. The correct output is given for XDMF

To Reproduce Run the following on the fenicsx branch

import festim as F
import numpy as np

L = 3e-04
vertices = np.linspace(0, L, num=1000)

my_model = F.HydrogenTransportProblem()
my_model.mesh = F.Mesh1D(vertices)

my_mat = F.Material(D_0=1.9e-7, E_D=0.2, name="my_mat")
my_subdomain = F.VolumeSubdomain1D(id=1, borders=[0, L], material=my_mat)
left_surface = F.SurfaceSubdomain1D(id=1, x=0)
right_surface = F.SurfaceSubdomain1D(id=2, x=L)
my_model.subdomains = [
    my_subdomain,
    left_surface,
    right_surface,
]

mobile_H = F.Species("H")
trapped_H = F.Species("trapped_H")
my_model.species = [mobile_H, trapped_H]

temperature = 500.0
my_model.temperature = temperature

my_model.boundary_conditions = [
    F.DirichletBC(subdomain=right_surface, value=0, species=mobile_H),
    F.DirichletBC(subdomain=left_surface, value=1e12, species=mobile_H),
]
my_model.exports = [
    F.VTXExport("mobile_concentration_h.bp", field=mobile_H),
    F.VTXExport("trapped_concentration_h.bp", field=trapped_H),
    F.XDMFExport("mobile_concentration_h.xdmf", field=mobile_H),
    F.XDMFExport("trapped_concentration_h.xdmf", field=trapped_H),
]

my_model.settings = F.Settings(
    atol=1e10,
    rtol=1e-10,
    max_iterations=30,
    final_time=50,
)

my_model.settings.stepsize = F.Stepsize(initial_value=1 / 20)

my_model.initialise()

from petsc4py import PETSc

my_model.solver.convergence_criterion = "incremental"
ksp = my_model.solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

my_model.run()
RemDelaporteMathurin commented 10 months ago

I tracked this down to the use of .collapse() for the post_processing_solution attribute of species.

https://github.com/RemDelaporteMathurin/FESTIM/blob/1dbe47ddd5dab75c54963f1ab570c3121e542131/festim/hydrogen_transport_problem.py#L288