Error with DG discretization and EMPM binding model

CADET Troubleshooting
1

Expected behavior.

Using the DG build of CADET-Core and CADET-Process, SMB built with CarouselBuilder using ExtendedMobilePhaseModulator binding model and DiscontinuousGalerkin discretization runs successfully with float beta values >=0.

Actual behavior

Simulation fails when any non-int value of beta is passed (e.g. 0.1, 0.5, 1.3, etc.). Int values for beta appear to work even when explicitly passed as floats (e.g. 0.0, 1.0, 2.0 etc.). I’ve tested with multiple unit operation models (LRM, LRMP, GRM) and the error seems to occur with all of them.

Simulations typically fail with the following error:

Traceback (most recent call last):
  File "<frozen runpy>", line 198, in _run_module_as_main
  File "<frozen runpy>", line 88, in _run_code
  File "/home/cadet/mre.py", line 123, in <module>
    run_sim()
  File "/home/cadet/mre.py", line 119, in run_sim
    results = simulator.simulate(_model)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/simulator/simulator.py", line 275, in simulate
    results = self.simulate_n_cycles(
              ^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/log.py", line 164, in wrapper
    result = function(*args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/log.py", line 225, in wrapper
    results = function(*args, **kwargs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/log.py", line 200, in wrapper
    raise e
  File "/venv/lib/python3.12/site-packages/CADETProcess/log.py", line 192, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/simulator/simulator.py", line 334, in simulate_n_cycles
    return self._run(process, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/simulator/cadetAdapter.py", line 154, in wrapper
    results = func(self, process, *args, **kwargs)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/venv/lib/python3.12/site-packages/CADETProcess/simulator/cadetAdapter.py", line 292, in _run
    raise CADETProcessError(
CADETProcess.CADETProcessError.CADETProcessError: CADET Error: Simulation failed with double free or corruption (out)

I’ve also observed other variations when parameters are tweaked:

CADET Error: Simulation failed with munmap_chunk(): invalid pointer

CADET Error: Simulation failed with free(): invalid size

How to produce bug (including a minimal reproducible example)

My development environment is in a devcontainer based on Ubuntu 24.04. I’ve also included my dockerfile below for reference as it includes my configuration for the DG build of CADET-Core.

MRE code:

from CADETProcess.processModel import (
    ComponentSystem,
    Inlet,
    Outlet,
    Cstr,
    LumpedRateModelWithoutPores,
    LRMDiscretizationDG,
    ExtendedMobilePhaseModulator
    )
from CADETProcess.modelBuilder import CarouselBuilder, SerialZone
from CADETProcess.simulator import Cadet

def run_sim():
    feed_conc = [1, 5]
    eluent_conc = [0.001, 0.001]

    component_system = ComponentSystem(["a", "b"])

    binding_model = ExtendedMobilePhaseModulator(component_system)
    binding_model.name = "empm"
    binding_model.is_kinetic = True
    binding_model.adsorption_rate = [10, 1]
    binding_model.desorption_rate = [3, 1]
    binding_model.capacity = [10, 1]
    binding_model.ion_exchange_characteristic = [0.1, 0]
    binding_model.hydrophobicity = [0, 0]
    binding_model.component_mode = [2, 0]

    column = LumpedRateModelWithoutPores(component_system, name="column")
    column.binding_model = binding_model
    column.total_porosity = 0.4
    column.length = 40 / 100
    column.diameter = 15 / 1000
    column.axial_dispersion = 5e-5
    column.c = eluent_conc
    column.q = [0.001, 0.001]
    column.discretization = LRMDiscretizationDG()

    eluent = Inlet(component_system, name='eluent')
    eluent.c = eluent_conc
    eluent.flow_rate = 4e-7

    cstr_init_vol = 0.03

    eluent_cstr = Cstr(component_system, name="eluent_cstr")
    eluent_cstr.init_liquid_volume = cstr_init_vol

    feed = Inlet(component_system, name='feed')
    feed.c = feed_conc
    feed.flow_rate = 4e-7

    feed_cstr = Cstr(component_system, name="feed_cstr")
    feed_cstr.init_liquid_volume = cstr_init_vol

    eluate = Outlet(component_system, name='eluate')
    eluate_cstr = Cstr(component_system, name="eluate_cstr")
    eluate_cstr.init_liquid_volume = cstr_init_vol

    depleted = Outlet(component_system, name='depleted')
    depleted_cstr = Cstr(component_system, name="depleted_cstr")
    depleted_cstr.init_liquid_volume = cstr_init_vol

    adsorb_zone = SerialZone(component_system, "adsorb_zone", 1)
    desorb_zone = SerialZone(component_system, "desorb_zone", 1)

    def model():
        builder = CarouselBuilder(component_system, "smb")
        builder.column = column

        adsorb_train = (feed, feed_cstr, adsorb_zone, depleted_cstr, depleted)
        desorb_train = (eluent, eluent_cstr, desorb_zone, eluate_cstr, eluate)

        # Add main connections
        for train in (adsorb_train, desorb_train):
            for unit in train:
                builder.add_unit(unit)
            
            # Make connections
            for this_unit, next_unit in zip(train, train[1:]):
                builder.add_connection(this_unit, next_unit)

        # Recycle connections
        builder.add_connection(eluate_cstr, feed_cstr)
        recycle_to_feed = 0.1
        builder.set_output_state(eluate_cstr, [1-recycle_to_feed, recycle_to_feed])
        
        builder.add_connection(depleted_cstr, eluent_cstr)
        depleted_to_eluent = 0.1
        builder.set_output_state(depleted_cstr, [1-depleted_to_eluent, depleted_to_eluent])

        builder.switch_time = 100

        return builder.build_process()

    _model = model()

    simulator = Cadet("/opt/CADET/install")
    simulator.n_cycles = 3
    simulator.timeout = 5 * 60
    results = simulator.simulate(_model)

    return results

run_sim()

Dockerfile:

# Base image
FROM mcr.microsoft.com/devcontainers/base:ubuntu-24.04

# Prep to compile CADET-Core
RUN apt-get update && apt-get install -y \
    build-essential \
    cmake \
    git \
    libhdf5-dev \
    libsuperlu-dev \
    libeigen3-dev \
    intel-mkl \
    python3.12 \
    python3.12-dev \
    python3.12-venv \
    python3-pip \
    && rm -rf /var/lib/apt/lists/*

RUN export MKLROOT=/opt/intel/mkl

RUN update-alternatives --install /usr/bin/python python /usr/bin/python3.12 1

# Create venv
RUN python -m venv /venv \
        && /venv/bin/python -m pip install --upgrade pip setuptools wheel \
        && chown -R vscode:vscode /venv
        
# Install requirements
COPY requirements.txt /tmp/requirements.txt
RUN /venv/bin/pip install --no-cache-dir -r /tmp/requirements.txt

# Pull and compile Cadet-Core
RUN git clone https://github.com/cadet/CADET-Core /opt/CADET
WORKDIR /opt/CADET
RUN mkdir build && mkdir install && cd build && \
    cmake -DCMAKE_INSTALL_PREFIX="../install" \
    -DENABLE_DG=ON \
    -DCMAKE_BUILD_TYPE=Release \
    -DCADET_BUILD_EXAMPLES=OFF .. ; \
    make -j$(nproc) && \
    make install

# Set python path
ENV PATH="/venv/bin:$PATH"

# Set up working dir and copy contents
RUN mkdir -p /home/cadet
WORKDIR /home/cadet

USER vscode

File produced by conda env export > environment.yml

requirements.txt (950 Bytes)

2

Hi Nick,

thank you for reporting this apparent bug.

Did you try running the simulation with a FV discretization? Maybe the DG discretization causes negative concentration values. In this case, c^\beta is not always defined when \beta is not an integer. The error message should definitely be improved.

For context: DG is more accurate than FV for sufficiently smooth concentration profiles, but can produce oscillations causing negative values when the solution is not smooth enough. Depending on the binding model, this can cause the simulation to crash, even when the negative values are small.

Best regards,
Jan

3

Hi Jan,

Thanks for the follow up. The simulation does run successfully with FV using the same parameters, and when I set the eluent and feed concentrations closer together (within ~1-2 orders of magnitude) it also runs successfully with DG. I had assumed that an SMB-style system with recycle loops by nature would result in concentration profiles which are smooth enough, but I think you’re right that this might not be the case for my system which explains the crashing.

Sounds like I’ll stick with FV then, thanks again for the support!

Nick

4

The concentration profiles could be smooth enough so that DG would yield good results in the sense of an absolute metric. But it might still give small negative values, which doesnt happen for the FV code, and which can crash the simulation if negative values are mathematically disallowed by the binding model.

We are aware of that problem and a stabilizing mechanism is planned for the DG code.

2 Likes
5

I debugged the simulation and could confirm that there are small negative values with DG, but not with FV. I double checked this because the error message CADET Error: Simulation failed with double free or corruption (out) is misleading. Ive added the issue Abort simulation and raise error if NaN is encountered to our next dev-call.

6

Thanks for confirming Jan. An improved error message would be quite helpful in this case.