Source code for opencmp.models.poisson

########################################################################################################################
# Copyright 2021 the authors (see AUTHORS file for full list).                                                         #
#                                                                                                                      #
# This file is part of OpenCMP.                                                                                        #
#                                                                                                                      #
# OpenCMP is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public  #
# License as published by the Free Software Foundation, either version 2.1 of the License, or (at your option) any     #
# later version.                                                                                                       #
#                                                                                                                      #
# OpenCMP is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied        #
# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more  #
# details.                                                                                                             #
#                                                                                                                      #
# You should have received a copy of the GNU Lesser General Public License along with OpenCMP. If not, see             #
# <https://www.gnu.org/licenses/>.                                                                                     #
########################################################################################################################

import ngsolve as ngs
from ..helpers.ngsolve_ import get_special_functions
from ..helpers.dg import jump, grad_avg
from . import Model
from typing import Dict, List, Union, Optional
from ngsolve.comp import FESpace, ProxyFunction
from ngsolve import Parameter, GridFunction, BilinearForm, LinearForm, Preconditioner


[docs]class Poisson(Model): """ A poisson model. """ def _pre_init(self) -> None: # Nothing needs to be done pass def _post_init(self) -> None: # Nothing needs to be done pass def _define_model_components(self) -> Dict[str, Optional[int]]: return {'u': 0} def _define_model_local_error_components(self) -> Dict[str, bool]: return {'u': True} def _define_time_derivative_components(self) -> List[Dict[str, bool]]: return [{'u': True}] def _define_num_weak_forms(self) -> int: return 1 def _define_bc_types(self) -> List[str]: return ['robin', 'dirichlet', 'neumann']
[docs] @staticmethod def allows_explicit_schemes() -> bool: return True
def _construct_fes(self) -> FESpace: fes = getattr(ngs, self.element['u'])(self.mesh, order=self.interp_ord, dirichlet=self.dirichlet_names.get('u', ''), dgjumps=self.DG) compound_fes = ngs.FESpace([fes]) return compound_fes def _set_model_parameters(self) -> None: self.dc = self.model_functions.model_parameters_dict['diffusion_coefficient']['all'] self.f = self.model_functions.model_functions_dict['source']['all']
[docs] def construct_bilinear_time_ODE(self, U: Union[List[ProxyFunction], List[GridFunction]], V: List[ProxyFunction], dt: Parameter = Parameter(1.0), time_step: int = 0) -> List[BilinearForm]: # Split up the trial (weighting) and test functions u = U[0] v = V[0] # Laplacian term a = dt * self.dc[time_step] * ngs.InnerProduct(ngs.Grad(u), ngs.Grad(v)) * ngs.dx # Bulk of Bilinear form if self.DG: # Define the special DG functions. n, _, alpha, I_mat = get_special_functions(self.mesh, self.nu) if self.dirichlet_names.get('u', None) is not None: # Penalty terms for Dirichlet BCs a += dt * self.dc[time_step] * ( - u * n * ngs.Grad(v) # 1/2 of penalty for u=g on 𝚪_D - ngs.Grad(u) * n * v # ∇u^ = ∇u + alpha * u * v # 1/2 of penalty term for u=g on 𝚪_D from ∇u^ ) * self._ds(self.dirichlet_names['u']) # Robin BCs for u for marker in self.BC.get('robin', {}).get('u', {}): r, q = self.BC['robin']['u'][marker][time_step] a += -dt * self.dc[time_step] * r * u * v * self._ds(marker) return [a]
[docs] def construct_bilinear_time_coefficient(self, U: List[ProxyFunction], V: List[ProxyFunction], dt: Parameter, time_step: int) -> List[BilinearForm]: # Split up the trial (weighting) and test functions u = U[0] v = V[0] if self.DG: # Define the special DG functions. n, _, alpha, I_mat = get_special_functions(self.mesh, self.nu) jump_u = jump(u) avg_grad_u = grad_avg(u) jump_v = jump(v) avg_grad_v = grad_avg(v) # Penalty for discontinuities a = dt * self.dc[time_step] * ( - jump_u * n * avg_grad_v # U - avg_grad_u * n * jump_v # 1/2 term for u+=u- on 𝚪_I from ∇u^ + alpha * jump_u * jump_v # 1/2 term for u+=u- on 𝚪_I from ∇u^ ) * ngs.dx(skeleton=True) else: a = ngs.CoefficientFunction(0.0) * ngs.dx return [a]
[docs] def construct_linear(self, V: List[ProxyFunction], gfu_0: Optional[List[GridFunction]], dt: Parameter, time_step: int) -> List[LinearForm]: # V is a list of length 1 for Poisson, get just the first element v = V[0] # Source term L = dt * self.f[time_step] * v * ngs.dx # Dirichlet BCs for u, they only get added if using DG if self.DG: # Define the special DG functions. n, _, alpha, I_mat = get_special_functions(self.mesh, self.nu) for marker in self.BC.get('dirichlet', {}).get('u', {}): g = self.BC['dirichlet']['u'][marker][time_step] L += dt * self.dc[time_step] * ( alpha * g * v # 1/2 of penalty term for u=g on 𝚪_D from ∇u^ - g * n * ngs.Grad(v) # 1/2 of penalty for u=g on 𝚪_D ) * self._ds(marker) # Neumann BCs for u for marker in self.BC.get('neumann', {}).get('u', {}): h = self.BC['neumann']['u'][marker][time_step] L += dt * self.dc[time_step] * h * v * self._ds(marker) # Robin BCs for u for marker in self.BC.get('robin', {}).get('u', {}): r, q = self.BC['robin']['u'][marker][time_step] L += -dt * self.dc[time_step] * r * q * v * self._ds(marker) return [L]
[docs] def construct_imex_explicit(self, V: List[ProxyFunction], gfu_0: Optional[List[GridFunction]], dt: Parameter, time_step: int) -> List[LinearForm]: # The Poisson equation is linear and has no need for IMEX schemes. pass
[docs] def solve_single_step(self, a_lst: List[BilinearForm], L_lst: List[LinearForm], precond_lst: List[Preconditioner], gfu: GridFunction, time_step: int = 0) -> None: self.construct_and_run_solver(a_lst[0], L_lst[0], precond_lst[0], gfu)
[docs] def solve_stationary_iteration(self, a_lst: List[BilinearForm], L_lst: List[LinearForm], precond_lst: List[Preconditioner], gfu: GridFunction) -> None: self.construct_and_run_solver(a_lst[0], L_lst[0], precond_lst[0], gfu) return(0., 0.)
[docs] def accept_stationary_solution(self, gfu: GridFunction): pass