###################################################################################
# Copyright 2021 National Technology & Engineering Solutions of Sandia, #
# LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the #
# U.S. Government retains certain rights in this software. #
# If you want to use this code, please refer to the README.rst and LICENSE files. #
###################################################################################
import numpy as np
from PyNucleus_base import REAL
from PyNucleus_base.utilsFem import problem
from . functions import complexLambda, wrapRealToComplexFunction, waveFunction, radialIndicator
from . import functionFactory
[docs]
class diffusionProblem(problem):
[docs]
def setDriverArgs(self):
p = self.driver.addGroup('problem')
self.setDriverFlag('domain', 'square', acceptedValues=['interval', 'cube', 'standardSimplex3D', 'fichera', 'gradedSquare', 'gradedCube',
'sphere1', 'sphere2'], group=p)
self.setDriverFlag('problem', 'sin', acceptedValues=['reac-sin', 'diffusivity-sin', 'poly', 'fichera', 'cos'], group=p)
self.setDriverFlag('noRef', argInterpreter=int, group=p)
self.setDriverFlag('element', 'P1', acceptedValues=['P1', 'P2', 'P3'], group=p)
self.setDriverFlag('symmetric', False, group=p)
self.setDriverFlag('reorder', False, group=p)
[docs]
def processCmdline(self, params):
domain = params['domain']
element = params['element']
noRef = params['noRef']
if domain in ('interval', 'unitInterval'):
if noRef is None:
noRef = {'P1': 15, 'P2': 14, 'P3': 13}[element]
elif domain in ('sphere1', ):
if noRef is None:
noRef = {'P1': 10, 'P2': 9, 'P3': 8}[element]
elif domain in ('square', 'unitSquare', 'gradedSquare'):
if noRef is None:
noRef = {'P1': 9, 'P2': 8, 'P3': 7}[element]
elif domain in ('square', 'unitSquare', 'gradedSquare'):
if noRef is None:
noRef = {'P1': 9, 'P2': 8, 'P3': 7}[element]
elif domain == 'graded_disc':
if noRef is None:
noRef = {'P1': 5, 'P2': 4, 'P3': 3}[element]
elif domain in ('sphere2', ):
if noRef is None:
noRef = {'P1': 5, 'P2': 4, 'P3': 3}[element]
elif domain in ('cube', 'gradedCube'):
if noRef is None:
noRef = {'P1': 6, 'P2': 5, 'P3': 4}[element]
elif domain == 'standardSimplex3D':
if noRef is None:
noRef = {'P1': 2}[element]
elif domain == 'fichera':
if noRef is None:
noRef = {'P1': 5, 'P2': 4}[element]
params['noRef'] = noRef
super().processCmdline(params)
[docs]
@problem.generates(['dim', 'manifold_dim', 'diffusivity', 'reaction', 'rhsFun', 'exactSolution', 'L2ex', 'H10ex', 'boundaryCond', 'nontrivialNullspace'])
def processProblem(self, domain, problem, noRef, element, symmetric, reorder):
from . functions import constant, Lambda
from . factories import (meshFactory, rhsFunSin1D,
rhsFunSin2D, rhsFunSin3D, solSin1D, solSin2D, solSin3D, cos2D,
rhsCos2D, rhsFichera, solFichera,)
self.diffusivity = None
self.reaction = None
self.dim = meshFactory.getDim(domain)
self.manifold_dim = meshFactory.getManifoldDim(domain)
self.nontrivialNullspace = False
if domain in ('interval', 'unitInterval'):
if problem == 'sin':
self.rhsFun = rhsFunSin1D
self.exactSolution = solSin1D
self.L2ex = 1/2
self.H10ex = np.pi**2/2
self.boundaryCond = None
elif problem == 'reac-sin':
self.rhsFun = Lambda(lambda x: (np.pi**2.0 + 10.)*np.sin(np.pi*x[0]))
self.exactSolution = solSin1D
self.L2ex = 1/2
self.H10ex = (np.pi**2 + 10.)/2
self.reaction = 10.
self.boundaryCond = None
else:
raise NotImplementedError()
elif domain in ('sphere1', ):
if problem == 'sin':
n = 1
self.exactSolution = functionFactory('Lambda', lambda x: np.sin(n*np.arctan2(x[1], x[0])))
self.rhsFun = n**2 * self.exactSolution
self.L2ex = np.pi
self.H10ex = np.pi*n**2
self.boundaryCond = None
self.nontrivialNullspace = True
else:
raise NotImplementedError()
elif domain in ('square', 'unitSquare', 'gradedSquare'):
if problem == 'sin':
self.rhsFun = rhsFunSin2D
self.exactSolution = solSin2D
self.L2ex = 1/4
self.H10ex = 2*np.pi**2/4
self.boundaryCond = None
elif problem == 'cos':
self.rhsFun = rhsCos2D
self.exactSolution = cos2D
self.L2ex = 1/4
self.H10ex = 2*np.pi**2/4
self.boundaryCond = cos2D
elif problem == 'reac-sin':
self.rhsFun = Lambda(lambda x: (2*np.pi**2.0 + 10.)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))
self.exactSolution = solSin2D
self.L2ex = 1/4
self.H10ex = (2*np.pi**2 + 10.)/4
self.boundaryCond = None
self.reaction = 10.
elif problem == 'diffusivity-sin':
self.diffusivity = Lambda(lambda x: np.exp(np.sin(np.pi*x[0]) *
np.sin(np.pi*x[1])))
self.rhsFun = Lambda(lambda x: -np.pi**2 *
np.exp(np.sin(np.pi*x[0])*np.sin(np.pi*x[1])) *
(np.sin(np.pi*x[0])**2 * np.cos(np.pi*x[1])**2 +
np.cos(np.pi*x[0])**2 * np.sin(np.pi*x[1])**2 -
2*np.sin(np.pi*x[0]) * np.sin(np.pi*x[1])))
self.exactSolution = solSin2D
self.L2ex = 1/4
self.H10ex = np.nan
self.boundaryCond = None
elif problem == 'poly':
self.rhsFun = Lambda(lambda x: 32*x[0]*(1-x[0])+32*x[1]*(1-x[1]))
self.exactSolution = Lambda(lambda x: 16*x[0]*x[1]*(1-x[0])*(1-x[1]))
self.L2ex = 256/900
self.H10ex = 256/45
self.boundaryCond = None
elif problem == 'variable-reac-sin':
self.rhsFun = Lambda(lambda x: (2*np.pi**2.0 + 10.)*np.sin(np.pi*x[0])*np.sin(np.pi*x[1]))
self.exactSolution = solSin2D
self.L2ex = 1/4
self.H10ex = (2*np.pi**2 + 10.)/4
self.boundaryCond = None
self.reaction = Lambda(lambda x: 0. if x[0] < 0.5 else 2000.)
else:
raise NotImplementedError()
elif domain == 'graded_disc':
if problem == 'constant':
self.rhsFun = constant(1.)
self.exactSolution = None
self.L2ex = None
self.H10ex = None
self.boundaryCond = None
else:
raise NotImplementedError()
elif domain in ('sphere2', ):
if problem == 'sin':
from scipy.special import sph_harm
n = 1
m = -1
assert n >= 0
assert abs(m) <= n
ev = n*(n+1)
if m == 0:
self.exactSolution = functionFactory('Lambda', lambda x: sph_harm(m, n, np.arccos(x[2]), np.arctan2(np.sqrt(x[0]**2+x[1]**2), x[2])).real)
elif m < 0:
self.exactSolution = functionFactory('Lambda',
lambda x: (np.sqrt(0.5)*1j*(sph_harm(m, n, np.arccos(x[2]), np.arctan2(np.sqrt(x[0]**2+x[1]**2), x[2]))
- (-1)**m * sph_harm(-m, n, np.arccos(x[2]),
np.arctan2(np.sqrt(x[0]**2+x[1]**2), x[2])))).real)
else:
self.exactSolution = functionFactory('Lambda', lambda x:
(np.sqrt(0.5)*(sph_harm(-m, n, np.arccos(x[2]), np.arctan2(np.sqrt(x[0]**2+x[1]**2), x[2]))
+ (-1)**m * sph_harm(m, n, np.arccos(x[2]),
np.arctan2(np.sqrt(x[0]**2+x[1]**2), x[2])))).real)
def fun(x):
r = np.linalg.norm(x)
theta = np.arccos(x[2] / r)
phi = np.sign(x[1]) * np.arccos(x[0]/np.sqrt(x[0]**2+x[1]**2))
# theta \\in [0, pi]
# phi \\in [-pi, pi]
if m == 0:
return sph_harm(m, n, phi, theta).real
elif m > 0:
return (np.sqrt(0.5)*1j*(sph_harm(m, n, phi, theta) - (-1)**m * sph_harm(-m, n, phi, theta))).real
else:
return np.sqrt(0.5)*(sph_harm(-m, n, phi, theta) + (-1)**m * sph_harm(m, n, phi, theta)).real
self.exactSolution = functionFactory('Lambda', fun)
self.rhsFun = ev * self.exactSolution
self.L2ex = 1.0
self.H10ex = ev
self.boundaryCond = None
self.nontrivialNullspace = True
else:
raise NotImplementedError()
elif domain in ('cube', 'gradedCube'):
if problem == 'sin':
self.rhsFun = rhsFunSin3D
self.exactSolution = solSin3D
self.L2ex = 1/8
self.H10ex = 3*np.pi**2/8
self.boundaryCond = None
elif problem == 'variable-reac-sin':
self.rhsFun = constant(1.)
self.exactSolution = None
self.L2ex = np.nan
self.H10ex = np.nan
self.boundaryCond = None
self.reaction = Lambda(lambda x: 0. if x[0] < 0.5 else 2000.)
else:
raise NotImplementedError()
elif domain == 'standardSimplex3D':
if problem == 'poly':
self.rhsFun = Lambda(lambda x: 2*(x[1]*x[2]+x[0]*x[2]+x[0]*x[1]))
self.L2ex = 1/8
self.H10ex = 3*np.pi**2/8
self.boundaryCond = None
else:
raise NotImplementedError()
elif domain == 'fichera':
if problem == 'fichera':
self.rhsFun = rhsFichera
self.exactSolution = solFichera
self.L2ex = None
# H10ex = 7/8**9.52031/4
self.H10ex = None
self.boundaryCond = solFichera
else:
raise NotImplementedError()
else:
raise NotImplementedError()
[docs]
class helmholtzProblem(problem):
[docs]
def setDriverArgs(self):
p = self.driver.addGroup('problem')
self.setDriverFlag('domain', acceptedValues=['square', 'interval', 'cube'], group=p)
self.setDriverFlag('problem', acceptedValues=['wave', 'greens'], group=p)
self.setDriverFlag('element', 'P1', acceptedValues=['P1'], group=p)
self.setDriverFlag('frequency', 40., group=p)
self.setDriverFlag('symmetric', False, group=p)
self.setDriverFlag('reorder', False, group=p)
[docs]
@problem.generates(['dim', 'noRef', 'solEx', 'rhs', 'boundaryCond'])
def processProblem(self, domain, problem, element, frequency, symmetric, reorder):
from . import meshFactory
self.dim = meshFactory.getDim(domain)
if domain == 'interval':
self.noRef = 7
def n(x):
if x[0] == 0:
return np.array([-1.], dtype=REAL)
elif x[0] == 1:
return np.array([1.], dtype=REAL)
else:
raise NotImplementedError()
if problem == 'wave':
xi = np.array([0.5], dtype=REAL)
self.solEx = complexLambda(lambda x: np.exp(1j*np.vdot(xi, x)))
self.rhs = complexLambda(lambda x: (np.vdot(xi, xi)-self.frequency**2) * self.solEx(x))
self.boundaryCond = complexLambda(lambda x: 1j*(np.vdot(xi, n(x))+self.frequency) * self.solEx(x))
elif problem == 'greens':
self.rhs = wrapRealToComplexFunction(radialIndicator(1e-2, np.array([0.5])))
self.solEx = None
self.boundaryCond = None
else:
raise NotImplementedError(problem)
elif domain == 'square':
self.noRef = 8
def n(x):
if x[1] == 0:
return np.array([0., -1.], dtype=REAL)
elif x[1] == 1.:
return np.array([0., 1.], dtype=REAL)
elif x[0] == 0.:
return np.array([-1., 0.], dtype=REAL)
elif x[0] == 1.:
return np.array([1., 0.], dtype=REAL)
else:
raise NotImplementedError()
if problem == 'wave':
xi = np.array([0.5, 0.25], dtype=REAL)
self.solEx = waveFunction(xi)
self.rhs = (np.vdot(xi, xi)-self.frequency**2) * self.solEx
self.boundaryCond = complexLambda(lambda x: 1j*(np.vdot(xi, n(x))+self.frequency) * self.solEx(x))
elif problem == 'greens':
self.rhs = wrapRealToComplexFunction(radialIndicator(1e-2, np.array([0.5, 0.5])))
self.solEx = None
self.boundaryCond = None
else:
raise NotImplementedError(problem)
elif domain == 'cube':
self.noRef = 6
def n(x):
if x[2] == 0:
return np.array([0., 0., -1.], dtype=REAL)
elif x[2] == 1.:
return np.array([0., 0., 1.], dtype=REAL)
elif x[1] == 0:
return np.array([0., -1., 0.], dtype=REAL)
elif x[1] == 1.:
return np.array([0., 1., 0.], dtype=REAL)
elif x[0] == 0.:
return np.array([-1., 0., 0.], dtype=REAL)
elif x[0] == 1.:
return np.array([1., 0., 0.], dtype=REAL)
else:
raise NotImplementedError()
if problem == 'wave':
xi = np.array([0.75, 0.5, 0.25], dtype=REAL)
self.solEx = waveFunction(xi)
self.rhs = (np.vdot(xi, xi)-self.frequency**2) * self.solEx
self.boundaryCond = complexLambda(lambda x: 1j*(np.vdot(xi, n(x))+self.frequency) * self.solEx(x))
elif problem == 'greens':
self.rhs = wrapRealToComplexFunction(radialIndicator(1e-1, np.array([0.5, 0.5, 0.5])))
self.solEx = None
self.boundaryCond = None
else:
raise NotImplementedError(problem)
else:
raise NotImplementedError(domain)