Équation de Poisson: Imposer le gradient complet comme condition aux limites via les multiplicateurs de Lagrange

J'ai un problème physique régi par l'équation de Poisson en deux dimensions J'ai des mesures des deux composantes du gradient ∂ u / ∂ x et ∂ u / ∂ y le long d'une partie de la frontière, Γ m , je voudrais donc imposer ∂

$$-\nabla^2 u = f(x,y), \; in \; \Omega$$ $\partial{u}/\partial{x}$$\partial{u}/\partial{y}$$\Gamma_m$ et se propager dans le champ lointain. $$\frac{\partial{u}}{\partial{x_i}}_0 = g_m, \; on \; \Gamma_m$$

La composante du gradient tangentiel, , je peux simplement intégrer puis appliquer à travers une condition de Dirichlet, telle que ∫Γm∂$\frac{\partial{u}}{\partial{x}}_{(t,0)}$ Afin d'imposer simultanément la composante normale, ∂

$$\int_{\Gamma_m}\frac{\partial{u}}{\partial{x}}_{(t,0)} \, ds = u_0$$ , j'ai compris que je devrais passer par des multiplicateurs de Lagrange.$\frac{\partial{u}}{\partial{x}}_{(n,0)}$

Je pense donc que la forme variationnelle est alors J'ai passé beaucoup de temps à essayer de le reconstituer à partir des informations sur des problèmes connexes tels que https://answers.launchpad.net/fenics/+question/212434https://answers.launchpad.net/fenics/+question / 216323

$$\int_\Omega \nabla{u} \cdot\nabla{v}\, dx - \lambda \int_{\Gamma_m} ( \frac{\partial{u}}{\partial{x}}_{(n,0)}-g_m ) v \,ds = \int_\Omega f\, v\, dx$$

mais je ne vois toujours pas où je vais mal. Jusqu'à présent, ma tentative de solution est la suivante:

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "Lagrange", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Create mesh function over cell facets
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

# Mark left boundary facets as subdomain 0
class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and x[0] < DOLFIN_EPS

Gamma_Left = LeftBoundary()
Gamma_Left.mark(boundary_parts, 0)

class FarField(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and ( (x[0] > 1.0-DOLFIN_EPS) \
               or (x[1]<DOLFIN_EPS) or (x[1]> 1.0-DOLFIN_EPS) )

Gamma_FF = FarField()
Gamma_FF.mark(boundary_parts, 1)

# Define boundary condition
u0 = Expression("sin(x[1]*pi)")
bcs = [DirichletBC(V, u0, Gamma_Left)]

# Define variational problem
(u, lmbd) = TrialFunctions(W)
(v, d) = TestFunctions(W)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Constant(0.0)
h = Constant(-4.0)
n = FacetNormal(mesh)

F = inner(grad(u), grad(v))*dx + d*dot(grad(u),n)*ds(0) + lmbd*dot(grad(v),n)*ds(0)-\
   (f*v*dx + g*v*ds(1) + h*d*ds(0) + lmbd*h*ds(0))

a = lhs(F)
L = rhs(F)

# Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
for bc in bcs: bc.apply(A, b)

w = Function(W)
solve(A, w.vector(), b, 'lu')
(u,lmbd) = w.split()

# Plot solution
plot(u, interactive=True)

qui s'exécute mais donne un résultat bruyant ne ressemblant pas du tout à une solution d'une équation de Poisson. Cela semble avoir quelque chose à voir avec les espaces fonctionnels combinés, mais je ne trouve pas l'erreur.
J'apprécierais toute aide ou indication dans la bonne direction - merci beaucoup déjà!
À la vôtre
Markus

Tout d'abord, un point général: vous ne pouvez pas prescrire des conditions aux limites arbitraires pour un opérateur différentiel partiel et vous attendre à ce que l'équation différentielle partielle (qui inclut toujours à la fois l'opérateur et les conditions aux limites) soit bien posée, c'est-à-dire qu'elle admette une solution unique qui dépend en permanence de la données - tout cela est une condition nécessaire pour réellement essayer de calculer quelque chose.

Selon l'opérateur, il existe souvent un certain nombre de conditions valides que vous pouvez imposer (pour avoir un avant-goût, voir la monographie en trois volumes de Lions et Magenes). Cependant, ce que vous essayez de faire (spécifiez le gradient complet, qui équivaut à la fois aux conditions de Dirichlet et de Neumann sur la même (partie de la) frontière pour un PDE elliptique du deuxième ordre) ne fait pas partie d'eux - c'est ce qu'on appelle un problème de Cauchy latéral, et est mal posé: il n'y a aucune garantie qu'une paire donnée de données aux limites admet une solution, et même s'il en existe une, il n'y a pas de stabilité par rapport aux petites perturbations dans les données. (En fait, c'est le problème d'origine mal posé au sens de Hadamard, et l'exemple classique pour lequel vous ne pouvez pas ignorer les conditions aux limites lorsque vous discutez de la bonne pose. Vous pouvez trouver un exemple explicite dans ses conférences sur le problème de Cauchy en différentiel partiel linéaire équations des années 1920.)

$(r,R)\times(a,b)$$x=r$$R$$x\to \infty$$y=a$$y=b$

1. Si vous pouvez imposer des conditions aux limites (Neumann, Robin, Dirichlet - dont vous auriez besoin pour fixer la constante dans l'intégration de la dérivée tangentielle, soit dit en passant), alors il suffit d'utiliser soit les composantes normales de votre gradient comme condition de Neumann (si vous pouvez fixer le mode constant) ou intégrez les composants tangentiels en tant que condition de Dirichlet. Étant donné que les deux conditions correspondent probablement à la même fonction, vous obtenez la même solution dans les deux cas.

2. $y=a$$y=b$$-\Delta u=f$$-\Delta u + \varepsilon \Delta^2 u=f$$\varepsilon>0$$H^2$$u$$u_\varepsilon$$u$$\varepsilon\to0$

   $H^2$

Vous ne pouvez pas vous attendre à ce que la solution à votre problème modifié soit une solution au problème de Poisson parce que vous devez changer le problème d'une manière ou d'une autre pour qu'il soit bien posé.

$$F(u, \lambda) = \int_\Omega \frac{1}{2}|\nabla u|^2 \;\mathrm{d}x - \int_\Omega f u \;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} g u \;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \lambda (u-u_\mathrm{D})\;\mathrm{d}S$$ $(u,\lambda) \in V\times L^2(\Gamma_\mathrm{N})$$V = \{v\in H^1;v|_{\Gamma_\mathrm{D}}=0\}$$\Gamma_\mathrm{D}$$u_\mathrm{D}$$\Gamma_\mathrm{N}$$F(u)$ $$0 = \mathrm{D}F(u)[v] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad \forall v\in V,$$ ${\Gamma_\mathrm{N}}$${\Gamma_\mathrm{D}}$ $$0 = \mathrm{D}F(u,\lambda)[v,\mu] = \int_\Omega \nabla u\cdot\nabla v\;\mathrm{d}x - \int_\Omega fv\;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} gv\;\mathrm{d}S \quad + \int_{\Gamma_\mathrm{N}} \lambda v\;\mathrm{d}S + \int_{\Gamma_\mathrm{N}} \mu (u-u_\mathrm{D})\;\mathrm{d}S \quad \forall (v,\mu)\in V\times L^2(\Gamma_\mathrm{N}),$$ $-\Delta u=f$$\frac{\partial u}{\partial\mathbf{n}} = g-\lambda$$\Gamma_\mathrm{N}$$\Gamma_\mathrm{N}$

$\lambda\ll|g|$

$\Gamma_\mathrm{D}$$v\in V$$\Gamma_\mathrm{D}$

La conclusion est que vous ne pouvez pas vous attendre à ce que le PDE de second ordre admette deux conditions aux limites indépendantes.

---

$F(u,\lambda)$$\mathrm{D}F(u,\lambda)[v,\mu]$derivative()

$F(u,\lambda)$$\lambda\in L^2(\Gamma_\mathrm{N})$$\lambda\in L^2(\Omega)$$\lambda\in\mathbb{R}$

Votre approche ne peut pas fonctionner, certainement à cause de la mise en œuvre et probablement à cause de votre formulation.

L'imposition de conditions Dirichlet dans dolfin définit finalement les DOF ​​correspondants de votre espace de test à zéro.

Ceci est un extrait du manuel fenics :

Chapitre 3.3.9 (fin): L'application d'une condition aux limites de Dirichlet à un système linéaire identifiera tous les degrés de liberté qui devraient être définis à la valeur donnée et modifiera le système linéaire de sorte que sa solution respecte la condition aux limites. Ceci est accompli en mettant à zéro et en insérant 1 sur la diagonale des lignes de la matrice correspondant aux valeurs de Dirichlet, et en insérant la valeur de Dirichlet dans l'entrée correspondante du vecteur de droite [...]

$v$$\Gamma_m$

En résumé, en utilisant la routine par défaut dans dolfin, vous ne pouvez pas appliquer à la fois Dirichlet et d'autres conditions sur la même frontière.

Cependant, avant d'essayer de résoudre ce problème dans votre implémentation, allez trouver les bons espaces de test pour votre formulation mathématique (comme vient de le mentionner @Jan Blechta.)