festim-dev / FESTIM

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

No need to project solute*S when temperature is not type "expression" #296

Closed RemDelaporteMathurin closed 3 years ago

RemDelaporteMathurin commented 3 years ago

When running sims with conservation of chemical potential, the computational cost increases due to these projections in the post processing when computing surface fluxes and exporting the retention to XDMF:

https://github.com/RemDelaporteMathurin/FESTIM/blob/41f11c588366d0555e0f3c3e712e67c10d88bab2/FESTIM/post_processing.py#L80-L82

and

https://github.com/RemDelaporteMathurin/FESTIM/blob/41f11c588366d0555e0f3c3e712e67c10d88bab2/FESTIM/post_processing.py#L411-L413

The latter could be avoided when the type of the temperature is not "expression". See the MWE below:

import FESTIM
from fenics import *

parameters = {
    "mesh_parameters": {
        "vertices": [0, 0.5, 1],
        },
    "materials": [
        {
            "borders": [0, 1],
            "D_0": 1,
            "E_D": 1,
            "S_0": 1,
            "E_S": 1,
            "thermal_cond": 1,
            "id": 1,
        }
        ],
    "traps": [],

    "boundary_conditions": [],

    "solving_parameters": {
        "type": "solve_stationary",
        "newton_solver": {
            "absolute_tolerance": 1e10,
            "relative_tolerance": 1e-8,
            "maximum_iterations": 50,
        }
    },

    "exports": {
        "derived_quantities": {
            "surface_flux": [
                {
                    "surfaces": [2],
                    "field": "solute"
                }
                ],
            "file": "derived_quantities.csv",
            "folder": 'Results/1D_Results',
            "nb_iterations_between_exports": 1,
        }
    }
}

# make the temperature type expression

parameters["temperature"] = {
    "type": "expression",
    "value": 1,
}

my_sim = FESTIM.Simulation(parameters)
my_sim.initialise()

solute = my_sim.u*my_sim.S
n = FacetNormal(my_sim.mesh)
# print(assemble(dot(grad(solute), n)*ds))  # doesn't work

# make the temperature type solve_stationary
parameters["temperature"] = {
    "type": "solve_stationary",
    "boundary_conditions": [
        {
            "type": "dc",
            "value": 1,
            "surfaces": [1, 2]
        }
    ]
}

my_sim = FESTIM.Simulation(parameters)
my_sim.initialise()

solute = my_sim.u*my_sim.S
n = FacetNormal(my_sim.mesh)
print(assemble(dot(grad(solute), n)*ds))  # works!

A simple test checking the type of the temperature should be sufficient to greatly optimise the computational time in some cases. Even better, it'd be super to sort out how to make it work with the expression type....

RemDelaporteMathurin commented 3 years ago

This is related to #234 @ehodille