SymPy script for PFHub Problem 7

Problem7-sympy.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

## Sympy code to generate expressions for PFHub Problem 7 (MMS)

from sympy import Symbol, symbols, simplify
from sympy import Eq, sin, cos, cosh, sinh, tanh, sqrt
from sympy.physics.vector import divergence, gradient, ReferenceFrame, time_derivative
from sympy.printing import ccode, pprint
from sympy.abc import kappa, S, t
import string

# Spatial coordinates: x=R[0], y=R[1], z=R[2]
R = ReferenceFrame('R')

# sinusoid amplitudes
A1, A2 = symbols('A1 A2')
B1, B2 = symbols('B1 B2')
C1, C2 = symbols('C1 C2')

# Define interface offset (alpha)
alpha = 1/4 + A1 * t * sin(B1 * R[0]) \
            + A2 * sin(B2 * R[0] + C2 * t)
print("\nInterface displacement:")
pprint(Eq(symbols('alpha'),
          alpha))

# Define the solution equation (eta)
eta = 1/2 * (1 - tanh((R[1] - alpha) /
                       sqrt(2*kappa)))
print("\nManufactured Solution:")
pprint(Eq(symbols('eta'),
          eta))

# Compute the initial condition
print("\nInitial condition:")
pprint(Eq(symbols('eta0'),
          eta.subs(t, 0)))

# Compute the source term from the equation of motion
S = simplify(time_derivative(eta, R)
             + 4 * eta * (eta - 1) * (eta - 1/2)
             - divergence(kappa * gradient(eta, R), R))
print("\nSource term:")
pprint(Eq(symbols('S'), S))


# === Check Results against @stvdwtt ===
dadx = A1*B1*t*cos(B1*R[0]) + A2*B2*cos(B2*R[0]+C2*t)
dadt = A1*sin(B1*R[0])      + A2*C2*cos(B2*R[0]+C2*t)
d2adx2 = -A1*B1**2*t*sin(B1*R[0]) - A2*B2**2*sin(B2*R[0]+C2*t)
Sdw = 1/(cosh((R[1]-alpha)/sqrt(2*kappa))**2)/(4*sqrt(kappa)) * (-2*sqrt(kappa)*tanh((R[1]-alpha)/sqrt(2*kappa))*(dadx)**2 \
                                                             + sqrt(2)*(dadt - kappa*d2adx2))
print("@tkphd and @stvdwtt agree?")
notZero = simplify(Sdw - S)
pprint("True" if (not notZero) else delta)


print("\nC codes (without types):")
CA = ccode(alpha).replace('R_x','x').replace('R_y','y')
CH = ccode(eta).replace('R_x','x').replace('R_y','y')
C0 = ccode(eta.subs(t, 0)).replace('R_x','x').replace('R_y','y')
CS = ccode(S).replace('R_x','x').replace('R_y','y')

print("\nalpha(x, y, t) {")
pprint(CA)
print("}\n\nmanufacturedSolution(x, y, t) {")
pprint(CH)
print("}\n\ninitialCondition(x, y, t) {")
pprint(C0)
print("}\n\nsourceTerm(x, y, t) {")
pprint(CS)
print("}")


'''
Returns:

Interface displacement:
α = A₁⋅t⋅sin(Rₓ⋅B₁) + A₂⋅sin(Rₓ⋅B₂ + C₂⋅t) + 0.25

Manufactured Solution:
              ⎛√2⋅(R_y - A₁⋅t⋅sin(Rₓ⋅B₁) - A₂⋅sin(Rₓ⋅B₂ + C₂⋅t) - 0.25)⎞      
η = - 0.5⋅tanh⎜────────────────────────────────────────────────────────⎟ + 0.5
              ⎝                          2⋅√κ                          ⎠      

Initial condition:
               ⎛√2⋅(R_y - A₂⋅sin(Rₓ⋅B₂) - 0.25)⎞      
η₀ = - 0.5⋅tanh⎜───────────────────────────────⎟ + 0.5
               ⎝              2⋅√κ             ⎠      

Source term:
     ⎛   ⎛           ⎛     2                     2                  ⎞                                                     2     ⎛√2⋅(-R_y + A₁⋅t⋅sin(Rₓ⋅B₁) + A₂⋅sin(Rₓ⋅B₂ + C₂⋅t) + 0.25)⎞⎞                                                    ⎞ ⎛    2⎛√2⋅(-R_y + A₁⋅t⋅sin(Rₓ⋅B₁) + A₂⋅sin(Rₓ⋅B₂ + C₂⋅t) + 0.25)⎞    ⎞ 
    -⎜√κ⋅⎜0.25⋅√2⋅√κ⋅⎝A₁⋅B₁ ⋅t⋅sin(Rₓ⋅B₁) + A₂⋅B₂ ⋅sin(Rₓ⋅B₂ + C₂⋅t)⎠ + 0.5⋅(A₁⋅B₁⋅t⋅cos(Rₓ⋅B₁) + A₂⋅B₂⋅cos(Rₓ⋅B₂ + C₂⋅t)) ⋅tanh⎜─────────────────────────────────────────────────────────⎟⎟ + 0.25⋅√2⋅(A₁⋅sin(Rₓ⋅B₁) + A₂⋅C₂⋅cos(Rₓ⋅B₂ + C₂⋅t))⎟⋅⎜tanh ⎜─────────────────────────────────────────────────────────⎟ - 1⎟ 
     ⎝   ⎝                                                                                                                      ⎝                           2⋅√κ                          ⎠⎠                                                    ⎠ ⎝     ⎝                           2⋅√κ                          ⎠    ⎠ 
S = ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
                                                                                                                                                              √κ                                                                                                                                                         

@tkphd and @stvdwtt agree?
True

C codes (without types):

alpha(x, y, t) {
A1*t*sin(x*B1) + A2*sin(x*B2 + C2*t) + 0.25
}

manufacturedSolution(x, y, t) {
-0.5*tanh((1.0L/2.0L)*sqrt(2)*(y - A1*t*sin(x*B1) - A2*sin(x*B2 + C2*t) - 0.25)/sqrt(kappa)) + 0.5
}

initialCondition(x, y, t) {
-0.5*tanh((1.0L/2.0L)*sqrt(2)*(y - A2*sin(x*B2) - 0.25)/sqrt(kappa)) + 0.5
}

sourceTerm(x, y, t) {
-(sqrt(kappa)*(0.25*sqrt(2)*sqrt(kappa)*(A1*pow(B1, 2)*t*sin(x*B1) + A2*pow(B2, 2)*sin(x*B2 + C2*t)) + 0.5*pow(A1*B1*t*cos(x*B1) + A2*B2*cos(x*B2 + C2*t), 
2)*tanh((1.0L/2.0L)*sqrt(2)*(-y + A1*t*sin(x*B1) + A2*sin(x*B2 + C2*t) + 0.25)/sqrt(kappa))) + 0.25*sqrt(2)*(A1*sin(x*B1) + A2*C2*cos(x*B2 + C2*t)))*(pow(t
anh((1.0L/2.0L)*sqrt(2)*(-y + A1*t*sin(x*B1) + A2*sin(x*B2 + C2*t) + 0.25)/sqrt(kappa)), 2) - 1)/sqrt(kappa)
}

'''