Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BUG] SurfaceFlux doesn't account for Soret effect #830

Closed
RemDelaporteMathurin opened this issue Jul 30, 2024 · 0 comments · Fixed by #831
Closed

[BUG] SurfaceFlux doesn't account for Soret effect #830

RemDelaporteMathurin opened this issue Jul 30, 2024 · 0 comments · Fixed by #831
Labels
bug Something isn't working
Milestone

Comments

@RemDelaporteMathurin
Copy link
Collaborator

Describe the bug
When running a simulation with the Soret effect, I realised that setting a SurfaceFlux derived quantity would result in the wrong values. This is because the soret flag in compute seems to be ignored.

To Reproduce
Simple 1D steady state simulation, checking the flux balance with soret activated.

import numpy as np
import festim as F

my_model = F.Simulation()
my_model.mesh = F.MeshFromVertices(vertices=np.linspace(0, 1, num=100))
my_model.materials = F.Material(id=1, D_0=1, E_D=0, Q=2)
my_model.T = F.Temperature(value=700 + 1000 * F.x)

my_model.boundary_conditions = [
    F.DirichletBC(surfaces=[1], value=10, field=0),
    F.DirichletBC(surfaces=[2], value=1, field=0),
]
my_model.settings = F.Settings(
    absolute_tolerance=1e-10, relative_tolerance=1e-10, transient=False, soret=True
)

flux_left = F.SurfaceFlux(field=0, surface=1)
flux_right = F.SurfaceFlux(field=0, surface=2)
my_model.exports = [F.DerivedQuantities([flux_left, flux_right])]

my_model.initialise()
my_model.run()


# these two should be equal in steady state
print(flux_left.data)
print(flux_right.data)

Produces:

[370.7625458744958]
[1.1058666646404112]

Expected behavior
The two fluxes should have the same absolute value.

To convince us even further, try and add assert False in the compute method with soret:

    def compute(self, soret=False):
        flux = f.assemble(
            self.prop * f.dot(f.grad(self.function), self.n) * self.ds(self.surface)
        )
        if soret and self.field in [0, "0", "solute"]:
            assert False
            flux += f.assemble(
                self.prop
                * self.function
                * self.Q
                / (k_B * self.T**2)
                * f.dot(f.grad(self.T), self.n)
                * self.ds(self.surface)
            )
        return flux

And you will see everything runs fine.

Moreover, I found this TODO comment:

# TODO need to support for soret flag in surface flux

@RemDelaporteMathurin RemDelaporteMathurin added the bug Something isn't working label Jul 30, 2024
@RemDelaporteMathurin RemDelaporteMathurin added this to the 1.3.1 milestone Jul 30, 2024
@RemDelaporteMathurin RemDelaporteMathurin linked a pull request Jul 30, 2024 that will close this issue
10 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant