##
## Example demonstrating the method of
## manufactured solutions (MMS) using SymPy.
## The PDE used is the 1D diffusion equation
## with a spatially variable diffusivity, e.g.
##   d/dx ( kappa(x) dT/dx ) = f, for x \in [-1,1]
## The solution (T), and diffusivity (kappa)
## are chosen. The RHS (f) is manufactured
## using symbolic differentiation
##

from sympy import *
from sympy import var, Plot

x = Symbol('x')

# Chosen solution
T     = 8*sin(x) * x**4 + 20

# Chosen coefficient
kappa = ( 4*tanh(10*x) + 6 ) * exp(0.8*x)

# Compute, fMS = Lu
#              = d/dx(kappa(x).dT/dx)
dT_dx       = diff( T, x )
kappa_dT_dx = kappa * dT_dx

fMS         = diff( kappa_dT_dx, x )

print('[MMS chosen solution]')
print('  T = ' + ccode(T) + '\n')

print('[MMS chosen coefficients]')
print('  kappa =' + ccode(kappa) + '\n')

print('[MMS right hand side]')
print('  f_MS = ' + ccode(fMS) + '\n')

print('[MMS boundary conditions]')

print('  Dirichlet:')
T_eval = T.subs({x:'-1.0'})
result = T_eval.evalf()
print('    T(x)|(x=-1.0) = ' + str(result))

T_eval = T.subs({x:'1.0'})
result = T_eval.evalf()
print('    T(x)|(x=+1.0) = ' + str(result))

print('  Neumann:')
kgradT_eval = kappa_dT_dx.subs({x:'-1.0'})
result      = kgradT_eval.evalf()
print('    kappa(x).dT(x)/dx|(x=-1.0) = ' + str(result))

kgradT_eval = kappa_dT_dx.subs({x:'1.0'})
result      = kgradT_eval.evalf()
print('    kappa(x).dT(x)/dx|(x=+1.0) = ' +  str(result))