search

festim-dev / FESTIM

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

[BUG] Can't use natural log in expressions #813

Open RemDelaporteMathurin opened 1 month ago

RemDelaporteMathurin commented 1 month ago

Describe the bug

@KulaginVladimir and I recently noticed that we can't use natural log in expressions (probably temperature, sources, boundary conditions) due to this issue

To Reproduce Run

import festim as F
import sympy as sp
import numpy as np

my_model = F.Simulation()

my_model.mesh = F.MeshFromVertices(np.linspace(0, 1))
my_model.materials = F.Material(1, 1, 0)

my_model.T = F.Temperature(sp.log(F.t))
my_model.settings = F.Settings(1e10, 1e-10)
my_model.dt = F.Stepsize(1)
my_model.initialise()

It produces

Moving new file over differing existing file:
src: /home/remidm/V-V-report/jitfailure-dolfin_expression_832dccb795d293387fe65678894ff9aa/error.log.b1a76b4574e94437bc5aa74cc38833e9
dst: /home/remidm/V-V-report/jitfailure-dolfin_expression_832dccb795d293387fe65678894ff9aa/error.log
backup: /home/remidm/V-V-report/jitfailure-dolfin_expression_832dccb795d293387fe65678894ff9aa/error.log.old
Backup file exists, overwriting.
------------------- Start compiler output ------------------------
/tmp/tmpmu0esvf9/dolfin_expression_832dccb795d293387fe65678894ff9aa.cpp: In member function 'virtual void dolfin::dolfin_expression_832dccb795d293387fe65678894ff9aa::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const':
/tmp/tmpmu0esvf9/dolfin_expression_832dccb795d293387fe65678894ff9aa.cpp:62:26: error: too few arguments to function 'void dolfin::log(int, std::string, ...)'
   62 |           values[0] = log(t);
      |                       ~~~^~~
In file included from /home/remidm/miniconda3/envs/vv-festim-report-env/include/dolfin/common/Array.h:32,
                 from /home/remidm/miniconda3/envs/vv-festim-report-env/include/dolfin/function/Expression.h:27,
                 from /tmp/tmpmu0esvf9/dolfin_expression_832dccb795d293387fe65678894ff9aa.cpp:13:
/home/remidm/miniconda3/envs/vv-festim-report-env/include/dolfin/log/log.h:103:8: note: declared here
  103 |   void log(int debug_level, std::string msg, ...);
      |        ^~~

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/remidm/V-V-report/jitfailure-dolfin_expression_832dccb795d293387fe65678894ff9aa
Traceback (most recent call last):
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dolfin/jit/jit.py", line 165, in compile_class
    module, signature = dijitso_jit(cpp_data, module_name, params,
                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dolfin/jit/jit.py", line 103, in dijitso_jit
    return dijitso.jit(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dijitso/jit.py", line 216, in jit
    raise DijitsoError("Dijitso JIT compilation failed, see '%s' for details"
dijitso.jit.DijitsoError: Dijitso JIT compilation failed, see '/home/remidm/V-V-report/jitfailure-dolfin_expression_832dccb795d293387fe65678894ff9aa' for details

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/remidm/V-V-report/report/validation/thermodesorption_spectra/dark/test2.py", line 13, in <module>
    my_model.initialise()
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/festim/generic_simulation.py", line 282, in initialise
    self.T.create_functions(self.mesh)
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/festim/temperature/temperature.py", line 39, in create_functions
    self.expression = f.Expression(sp.printing.ccode(self.value), t=0, degree=2)
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dolfin/function/expression.py", line 400, in __init__
    self._cpp_object = jit.compile_expression(cpp_code, params)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dolfin/function/jit.py", line 158, in compile_expression
    expression = compile_class(cpp_data, mpi_comm=mpi_comm)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/vv-festim-report-env/lib/python3.11/site-packages/dolfin/jit/jit.py", line 170, in compile_class
    raise RuntimeError("Unable to compile C++ code with dijitso")
RuntimeError: Unable to compile C++ code with dijitso

A solution is provided in the FEniCS discourse ticket, simply replace the string "log" by "std::log".

RemDelaporteMathurin commented 1 month ago

A workaround for now is to do:

import festim as F
import sympy as sp
import numpy as np

my_model = F.Simulation()

my_model.mesh = F.MeshFromVertices(np.linspace(0, 1))
my_model.materials = F.Material(1, 1, 0)

log_bis = sp.Function("std::log")
my_model.T = F.Temperature(log_bis(F.t))

my_model.settings = F.Settings(1e10, 1e-10)
my_model.dt = F.Stepsize(1)
my_model.initialise()
RemDelaporteMathurin commented 1 month ago

Actually the code does not work: 

import festim as F
import sympy as sp
import numpy as np

my_model = F.Simulation()

my_model.mesh = F.MeshFromVertices(np.linspace(0, 1))
my_model.materials = F.Material(1, 1, 0)

log_bis = sp.Function("std::log")
my_model.T = F.Temperature(log_bis(F.t))

my_model.settings = F.Settings(1e10, 1e-10, final_time=10)
my_model.dt = F.Stepsize(1)
my_model.initialise()

produces 

Traceback (most recent call last):
  File "/home/remidm/FESTIM-SurfaceKinetics-Validation/tmp_folder_not_working_cases/D_EUROFER/test.py", line 15, in <module>
    my_model.initialise()
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/festim/generic_simulation.py", line 297, in initialise
    raise AttributeError(
AttributeError: final_time argument must be provided to settings in transient simulations
(festim-surface-kinetics-vv-env) remidm@lap-mathurin-01:~/FESTIM-SurfaceKinetics-Validation/tmp_folder_not_working_cases/D_EUROFER$ python test.py 
Traceback (most recent call last):
  File "/home/remidm/FESTIM-SurfaceKinetics-Validation/tmp_folder_not_working_cases/D_EUROFER/test.py", line 15, in <module>
    my_model.initialise()
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/festim/generic_simulation.py", line 330, in initialise
    self.T.create_functions(self.mesh)
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/festim/temperature/temperature.py", line 41, in create_functions
    self.expression = f.Expression(sp.printing.ccode(self.value), t=0, degree=2)
                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/sympy/printing/codeprinter.py", line 752, in ccode
    return c_code_printers[standard.lower()](settings).doprint(expr, assign_to)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/sympy/printing/codeprinter.py", line 172, in doprint
    lines = self._print(expr).splitlines()
            ^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/sympy/printing/printer.py", line 331, in _print
    return printmethod(expr, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/sympy/printing/codeprinter.py", line 461, in _print_Function
    return self._print_not_supported(expr)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remidm/miniconda3/envs/festim-surface-kinetics-vv-env/lib/python3.11/site-packages/sympy/printing/codeprinter.py", line 582, in _print_not_supported
    raise PrintMethodNotImplementedError("Unsupported by %s: %s" % (str(type(self)), str(type(expr))) + \
sympy.printing.codeprinter.PrintMethodNotImplementedError: Unsupported by <class 'sympy.printing.c.C99CodePrinter'>: std::log
Set the printer option 'strict' to False in order to generate partially printed code.
RemDelaporteMathurin commented 1 month ago

The only way I found to circumvent this issue is by monkey patching the sympy C printer:

import festim as F
import sympy as sp
import numpy as np

my_model = F.Simulation()

my_model.mesh = F.MeshFromVertices(np.linspace(0, 1))
my_model.materials = F.Material(1, 1, 0)

log_bis = sp.Function("std::log")

# Monkey patch the C99CodePrinter class
from sympy.printing.c import C99CodePrinter

original_print_function = C99CodePrinter._print_Function

def _custom_print_Function(self, expr):
    if expr.func == log_bis:
        return f"std::log({self._print(expr.args[0])})"
    return original_print_function(self, expr)

C99CodePrinter._print_Function = _custom_print_Function

my_model.T = F.Temperature(log_bis(F.t))

my_model.settings = F.Settings(1e10, 1e-10, final_time=10)
my_model.dt = F.Stepsize(1)
my_model.initialise()
RemDelaporteMathurin commented 1 month ago

A fix would be to replace this line: https://github.com/festim-dev/FESTIM/blob/47d7dc4831fe51b7bf2fe197e973d6ba5d4f6987/festim/temperature/temperature.py#L41

by 

ccode = sp.printing.ccode(self.value).replace("log", "std::log")
self.expression = f.Expression(ccode , t=0, degree=2)

we'd have to propagate this fix in all the codebase where sp.printing.ccode is used though.

Maybe we could make a festim.ccode function to simplify things:

def ccode(string):
    return sp.printing.ccode(self.value).replace("log", "std::log")

But maybe not really needed.

OR

Monkey patch sp.printing.ccode to replace do the same thing.