diff --git a/mwe.py b/mwe.py new file mode 100644 index 00000000..c52da7f8 --- /dev/null +++ b/mwe.py @@ -0,0 +1,32 @@ +from src.custom_pathsim_blocks import FestimWall + +from pathsim import Simulation, Connection +import pathsim.blocks +import numpy as np +import matplotlib.pyplot as plt + + +# Create blocks +blocks, events = [], [] + +# Create Festim wall +festim_wall = FestimWall(thickness=1, D_0=1, E_D=0, temperature=300, n_vertices=100) +source_c0 = pathsim.blocks.Constant(0) +source_cL = pathsim.blocks.Constant(1) + +scope = pathsim.blocks.Scope(labels=["flux_0", "flux_L"]) + +blocks = [festim_wall, source_c0, source_cL, scope] + +connections = [ + Connection(source_c0, festim_wall[0]), + Connection(source_cL, festim_wall[1]), + Connection(festim_wall[0], scope[0]), + Connection(festim_wall[1], scope[1]), +] + +simulation = Simulation(blocks=blocks, connections=connections, dt=0.005) +simulation.run(0.1) + +scope.plot(marker="o") +plt.show() diff --git a/saved_graphs/festim_example.json b/saved_graphs/festim_example.json new file mode 100644 index 00000000..b6b0abde --- /dev/null +++ b/saved_graphs/festim_example.json @@ -0,0 +1,208 @@ +{ + "nodes": [ + { + "id": "4", + "type": "constant", + "position": { + "x": 1093.1885987981248, + "y": 271.20138523630243 + }, + "data": { + "label": "cL", + "value": "0" + }, + "measured": { + "width": 205, + "height": 53 + }, + "selected": false, + "dragging": false + }, + { + "id": "5", + "type": "scope", + "position": { + "x": 1271, + "y": 457 + }, + "data": { + "label": "scope 5" + }, + "measured": { + "width": 120, + "height": 140 + }, + "selected": false, + "dragging": false + }, + { + "id": "9", + "type": "scope", + "position": { + "x": 967.9637725351604, + "y": 81.60628755860085 + }, + "data": { + "label": "scope 9" + }, + "measured": { + "width": 120, + "height": 140 + }, + "selected": false, + "dragging": false + }, + { + "id": "10", + "type": "stepsource", + "position": { + "x": 589, + "y": 224 + }, + "data": { + "label": "stepsource 10", + "amplitude": "1", + "delay": "0.2" + }, + "measured": { + "width": 120, + "height": 120 + }, + "selected": false, + "dragging": false + }, + { + "id": "12", + "type": "wall", + "position": { + "x": 952, + "y": 417 + }, + "data": { + "label": "wall", + "thickness": "1", + "surface_area": "1", + "temperature": "300", + "D_0": "0.05", + "E_D": "0", + "n_vertices": "100" + }, + "measured": { + "width": 70, + "height": 200 + }, + "selected": true, + "dragging": false + } + ], + "edges": [ + { + "id": "e4-9", + "source": "4", + "target": "9", + "sourceHandle": null, + "targetHandle": null, + "type": "smoothstep", + "data": {}, + "style": { + "strokeWidth": 2, + "stroke": "#ECDFCC" + }, + "markerEnd": { + "type": "arrowclosed", + "width": 20, + "height": 20, + "color": "#ECDFCC" + } + }, + { + "id": "e10-9", + "source": "10", + "target": "9", + "sourceHandle": null, + "targetHandle": null, + "type": "smoothstep", + "data": {}, + "style": { + "strokeWidth": 2, + "stroke": "#ECDFCC" + }, + "markerEnd": { + "type": "arrowclosed", + "width": 20, + "height": 20, + "color": "#ECDFCC" + } + }, + { + "id": "e10-12-to_c_0", + "source": "10", + "target": "12", + "sourceHandle": null, + "targetHandle": "c_0", + "type": "smoothstep", + "data": {}, + "style": { + "strokeWidth": 2, + "stroke": "#ECDFCC" + }, + "markerEnd": { + "type": "arrowclosed", + "width": 20, + "height": 20, + "color": "#ECDFCC" + } + }, + { + "id": "e4-12-to_c_L", + "source": "4", + "target": "12", + "sourceHandle": null, + "targetHandle": "c_L", + "type": "smoothstep", + "data": {}, + "style": { + "strokeWidth": 2, + "stroke": "#ECDFCC" + }, + "markerEnd": { + "type": "arrowclosed", + "width": 20, + "height": 20, + "color": "#ECDFCC" + } + }, + { + "id": "e12-5-from_flux_L", + "source": "12", + "target": "5", + "sourceHandle": "flux_L", + "targetHandle": null, + "type": "smoothstep", + "data": {}, + "style": { + "strokeWidth": 2, + "stroke": "#ECDFCC" + }, + "markerEnd": { + "type": "arrowclosed", + "width": 20, + "height": 20, + "color": "#ECDFCC" + } + } + ], + "nodeCounter": 13, + "solverParams": { + "dt": "0.02", + "dt_min": "1e-6", + "dt_max": "1.0", + "Solver": "SSPRK22", + "tolerance_fpi": "1e-6", + "iterations_max": "100", + "log": "true", + "simulation_duration": "2", + "extra_params": "{}" + }, + "globalVariables": [] +} \ No newline at end of file diff --git a/src/App.jsx b/src/App.jsx index e0fd06b9..aee0c912 100644 --- a/src/App.jsx +++ b/src/App.jsx @@ -29,6 +29,7 @@ import { makeEdge } from './CustomEdge'; import MultiplierNode from './MultiplierNode'; import { Splitter2Node, Splitter3Node } from './Splitters'; import BubblerNode from './BubblerNode'; +import WallNode from './WallNode'; // Add nodes as a node type for this script const nodeTypes = { @@ -49,6 +50,7 @@ const nodeTypes = { pid: DefaultNode, splitter2: Splitter2Node, splitter3: Splitter3Node, + wall: WallNode, bubbler: BubblerNode, white_noise: SourceNode, pink_noise: SourceNode, @@ -609,6 +611,9 @@ export default function App() { case 'splitter3': nodeData = { ...nodeData, f1: '1/3', f2: '1/3', f3: '1/3' }; break; + case 'wall': + nodeData = { ...nodeData, thickness: '', surface_area: '1', temperature: '', D_0: '1', E_D: '0', n_vertices: '100' }; + break; case 'bubbler': nodeData = { ...nodeData, conversion_efficiency: '0.95', vial_efficiency: '0.9', replacement_times: '' }; break; diff --git a/src/WallNode.jsx b/src/WallNode.jsx new file mode 100644 index 00000000..3827d5ed --- /dev/null +++ b/src/WallNode.jsx @@ -0,0 +1,68 @@ +import React from 'react'; +import { Handle } from '@xyflow/react'; + +export default function WallNode({ data }) { + return ( +
+
{data.label}
+ + {/* Left side handles with labels */} + +
+ c₀ +
+ + +
+ φ₀ +
+ + {/* Right side handles with labels */} + +
+ cₗ +
+ + +
+ φₗ +
+
+ ); +} diff --git a/src/convert_to_python.py b/src/convert_to_python.py index 29b45b75..e1e1aa2b 100644 --- a/src/convert_to_python.py +++ b/src/convert_to_python.py @@ -3,11 +3,7 @@ from inspect import signature from pathsim.blocks import Scope -from .custom_pathsim_blocks import ( - Process, - Splitter, - Bubbler, -) +from .custom_pathsim_blocks import Process, Splitter, Bubbler, FestimWall from .pathsim_utils import ( map_str_to_object, make_blocks, @@ -162,11 +158,29 @@ def make_edge_data(data: dict) -> list[dict]: raise ValueError( f"Invalid source handle '{edge['sourceHandle']}' for {edge}." ) + elif isinstance(block, FestimWall): + if edge["sourceHandle"] == "flux_0": + output_index = 0 + elif edge["sourceHandle"] == "flux_L": + output_index = 1 + else: + raise ValueError( + f"Invalid source handle '{edge['sourceHandle']}' for {edge}." + ) else: output_index = 0 if isinstance(target_block, Scope): input_index = target_block._connections_order.index(edge["id"]) + elif isinstance(target_block, FestimWall): + if edge["targetHandle"] == "c_0": + input_index = 0 + elif edge["targetHandle"] == "c_L": + input_index = 1 + else: + raise ValueError( + f"Invalid target handle '{edge['targetHandle']}' for {edge}." + ) else: input_index = block_to_input_index[target_block] diff --git a/src/custom_pathsim_blocks.py b/src/custom_pathsim_blocks.py index 1165c40e..a9bdf5fc 100644 --- a/src/custom_pathsim_blocks.py +++ b/src/custom_pathsim_blocks.py @@ -188,3 +188,101 @@ def create_reset_events(self) -> list[pathsim.blocks.Schedule]: events.extend(self._create_reset_events_one_vial(vial, reset_times[i])) return events + + +# FESTIM wall + + +class FestimWall(Block): + def __init__( + self, thickness, temperature, D_0, E_D, surface_area=1, n_vertices=100 + ): + super().__init__() + try: + import festim as F + except ImportError: + raise ImportError("festim is needed for FestimWall node.") + + self.thickness = thickness + self.temperature = temperature + self.surface_area = surface_area + self.D_0 = D_0 + self.E_D = E_D + self.n_vertices = n_vertices + self.t = 0.0 + + self.initialise_festim_model() + + def initialise_festim_model(self): + import festim as F + + model = F.HydrogenTransportProblem() + + model.mesh = F.Mesh1D( + vertices=np.linspace(0, self.thickness, num=self.n_vertices) + ) + material = F.Material(D_0=self.D_0, E_D=self.E_D) + + vol = F.VolumeSubdomain1D(id=1, material=material, borders=[0, self.thickness]) + left_surf = F.SurfaceSubdomain1D(id=1, x=0) + right_surf = F.SurfaceSubdomain1D(id=2, x=self.thickness) + + model.subdomains = [vol, left_surf, right_surf] + + H = F.Species("H") + model.species = [H] + + model.boundary_conditions = [ + F.FixedConcentrationBC(left_surf, value=0.0, species=H), + F.FixedConcentrationBC(right_surf, value=0.0, species=H), + ] + + model.temperature = self.temperature + + model.settings = F.Settings( + atol=1e-10, rtol=1e-10, transient=True, final_time=1 + ) + + model.settings.stepsize = F.Stepsize(initial_value=1) + + self.surface_flux_0 = F.SurfaceFlux(field=H, surface=left_surf) + self.surface_flux_L = F.SurfaceFlux(field=H, surface=right_surf) + model.exports = [self.surface_flux_0, self.surface_flux_L] + + model.show_progress_bar = False + + model.initialise() + + self.dt = model.dt + self.c_0 = model.boundary_conditions[0].value_fenics + self.c_L = model.boundary_conditions[1].value_fenics + + self.model = model + + def update_festim_model(self, c_0, c_L, stepsize): + self.c_0.value = c_0 + self.c_L.value = c_L + self.dt.value = stepsize + + self.model.iterate() + + return self.surface_flux_0.data[-1], self.surface_flux_L.data[-1] + + def update(self, t): + # no internal algebraic operator -> early exit + # if self.op_alg is None: + # return 0.0 + + # block inputs + c_0, c_L = self.inputs.to_array() + + if t == 0.0: + flux_0, flux_L = 0, 0 + else: + flux_0, flux_L = self.update_festim_model( + c_0=c_0, c_L=c_L, stepsize=t - self.t + ) + + flux_0 *= self.surface_area + flux_L *= self.surface_area + return self.outputs.update_from_array([flux_0, flux_L]) diff --git a/src/pathsim_utils.py b/src/pathsim_utils.py index a7aaf476..0274f377 100644 --- a/src/pathsim_utils.py +++ b/src/pathsim_utils.py @@ -20,7 +20,14 @@ Schedule, ) from pathsim.blocks.noise import WhiteNoise, PinkNoise -from .custom_pathsim_blocks import Process, Splitter, Splitter2, Splitter3, Bubbler +from .custom_pathsim_blocks import ( + Process, + Splitter, + Splitter2, + Splitter3, + Bubbler, + FestimWall, +) from flask import jsonify import inspect @@ -49,6 +56,7 @@ "function": Function, "delay": Delay, "bubbler": Bubbler, + "wall": FestimWall, "white_noise": WhiteNoise, "pink_noise": PinkNoise, } @@ -419,11 +427,30 @@ def make_connections(nodes, edges, blocks) -> list[Connection]: raise ValueError( f"Invalid source handle '{edge['sourceHandle']}' for {edge}." ) + elif isinstance(block, FestimWall): + if edge["sourceHandle"] == "flux_0": + output_index = 0 + elif edge["sourceHandle"] == "flux_L": + output_index = 1 + else: + raise ValueError( + f"Invalid source handle '{edge['sourceHandle']}' for {edge}." + ) else: output_index = 0 if isinstance(target_block, Scope): input_index = target_block._connections_order.index(edge["id"]) + # TODO we should do the same for all blocks with several input/target ports + elif isinstance(target_block, FestimWall): + if edge["targetHandle"] == "c_0": + input_index = 0 + elif edge["targetHandle"] == "c_L": + input_index = 1 + else: + raise ValueError( + f"Invalid target handle '{edge['targetHandle']}' for {edge}." + ) else: input_index = block_to_input_index[target_block]