Surface and Volume Exports in HydrogenTransportProblem DiscontinuousChangeVar

User Support
1

Hi all,

In HydrogenTransportProblemDiscontinuousChangeVar, the concentration fields themselves appear to be exported as I would expect (atoms/m \!^{3-d} , where d is the dimension of the simulation). However, the TotalVolume and SurfaceFlux exports appear to use the chemical potential instead of the concentration, and therefore are off by a factor of the solubility. Was this intended?

Here is an MWE that compares the behavior to the base HydrogenTransportProblem:

L = 1.0
T0 = 500.0
D0, ED = 1.0, 0.0
KS0 = 5.0
c_left, c_right = 1.0, 0.0
vertices = np.concatenate([np.linspace(0.0, L, 200)])


def build_common(model):
    model.mesh = F.Mesh1D(vertices)

    mat = F.Material(D_0=D0, E_D=ED, K_S_0=KS0, E_K_S=0.0)

    vol = F.VolumeSubdomain1D(id=3, borders=[0.0, L], material=mat)
    left_surf = F.SurfaceSubdomain1D(id=1, x=vertices[0])
    right_surf = F.SurfaceSubdomain1D(id=2, x=vertices[-1])

    model.subdomains = [vol, left_surf, right_surf]

    H = F.Species("H")
    H.subdomains = [vol]
    model.species = [H]

    model.boundary_conditions = [
        F.FixedConcentrationBC(left_surf, value=c_left, species=H),
        F.FixedConcentrationBC(right_surf, value=c_right, species=H),
    ]
    model.temperature = T0

    flux = F.SurfaceFlux(field=H, surface=right_surf)
    inv = F.TotalVolume(field=H, volume=vol)
    model.exports = [flux, inv]
    return model, flux, inv


def run_continuous():
    model = F.HydrogenTransportProblem()
    model, flux, inv = build_common(model)
    model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
    model.initialise()
    model.run()
    return flux.value, inv.value


def run_changevar():
    model = F.HydrogenTransportProblemDiscontinuousChangeVar()
    model, flux, inv = build_common(model)
    model.settings = F.Settings(atol=1e-10, rtol=1e-10, transient=False)
    model.initialise()
    model.run()
    return flux.value, inv.value


flux_disc, inv_disc = run_continuous()
flux_cv, inv_cv= run_changevar()


def report(label, disc, cv):
    print(f"{label}")
    print(f"  Continuous : {disc: .6e}")
    print(f"  ChangeVar  : {cv: .6e}")
    print(f"  ratio      : {disc / cv: .6e}")
    print()

report("Permeation flux downstream (x = L):", flux_disc, flux_cv)
report("Inventory:", inv_disc, inv_cv)

This prints

Permeation flux downstream (x = L):
  Continuous :  1.000000e+00
  ChangeVar  :  2.000000e-01
  ratio      :  5.000000e+00

Inventory:
  Continuous :  5.000000e-01
  ChangeVar  :  1.000000e-01
  ratio      :  5.000000e+00

I think I could work around this by using the CustomQuantity class, but was curious to get some clarification on whether or not this is was intentional. If it’s a bug, I can open an issue on the repo. Thanks!