Convertion to C code of derivatives in NL Poisson Problem

1

Hi everyone, I am trying to resolve an equation with non linear coefficient, but as soon as I change my initial function u to a non constant function I receive the following error:

f = // Not supported in C:
// Derivative
// Derivative
// im
// re
((10*sin(20*x[0])*sin(20*x[1]) + 10.1)*(re(x[0])*Derivative(re(x[0]), x[0]) + im(x[0])*Derivative(im(x[0]), x[0])) - 200*(fabs(x[0]) + 1)*sin(20*x[1])*cos(20*x[0])*fabs(x[0]))/(pow(fabs(x[0]) + 1, 2)*fabs(x[0]))

here is the code I used for this problem, Iā€™m working on Ubuntu 18.04

""" Fenics Code used to search the solution of Richard's Equation in a steady state (no time dependence)."""

# Packages used
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
from sympy import sin, cos, symbols, lambdify
import sympy as sym

################################################################
# Defining Relativity hydraulic conductivity
# Defining k_s(x), positive function.
def k_s(x,y):
    return 0.1+10*(sin(20*(x))*sin(20*y)+1)
    

# Haverkamp model
# Constants used
A=1
B=1
Gamma=1

def k(x,y,u):
    return k_s(x,y)*(A/(A+((abs(u)/B)**Gamma)))
    
# Use SymPy to compute f from the manufactured solution u
from sympy import re, im, Derivative, symbols, lambdify
import sympy as sym
x, y = sym.symbols('x[0], x[1]')
u = x
f = -sym.diff(k(x,y,u)*sym.diff(u, x), x) -sym.diff(k(x,y,u)*sym.diff(u, y), y)
f = sym.simplify(f)
u_code = sym.printing.ccode(u)
f_code = sym.printing.ccode(f)
print('u =', u_code)
print('f =', f_code)

Thanks for your answers.

2

I usually do manufactured solutions by using UFL to take the derivatives, which avoids having to go through sympy. In your example, it would look something like the following:

from fenics import *

################################################################
# Defining Relativity hydraulic conductivity
# Defining k_s(x), positive function.
def k_s(x,y):
    return 0.1+10*(sin(20*x)*sin(20*y)+1)
    
# Haverkamp model
# Constants used
A=1
B=1
Gamma=1

def k(x,y,u):
    return k_s(x,y)*(A/(A+((abs(u)/B)**Gamma)))
    
mesh = UnitSquareMesh(10,10)
x,y = SpatialCoordinate(mesh)
u = x

# UFL representation of the source term that can be used directly in
# variational problems:
f = -div(k(x,y,u)*grad(u))
2 Likes
3

Thank you so much :slightly_smiling_face: