Ecuación de Poisson: imponer el gradiente completo como condición de frontera a través de multiplicadores de Lagrange

11

Tengo un problema físico regido por la ecuación de Poisson en dos dimensiones Tengo mediciones de los dos componentes de gradienteu /x yu /y a lo largo de alguna parte del límite, Γ m , por lo que me gustaría imponer u

2u=f(x,y),inΩ
u/xu/yΓm propagarse al campo lejano.
uxi0=gm,onΓm

El componente de gradiente tangencial, , solo puedo integrar y luego aplicar a través de una condición de Dirichlet, de modo que Γmuux(t,0) Para imponer simultáneamente el componente normal,u

Γmux(t,0)ds=u0
, deduje que tendría que ir a través de multiplicadores de Lagrange.ux(n,0)

Entonces creo que la forma variacional es entonces Pasé mucho tiempo tratando de reconstruirlo a partir de la información sobre problemas relacionados, como https://answers.launchpad.net/fenics/+question/212434https://answers.launchpad.net/fenics/+question / 216323

ΩuvdxλΓm(ux(n,0)gm)vds=Ωfvdx

pero todavía no puedo ver a dónde me estoy equivocando. Mi intento de solución hasta ahora es:

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)

que se ejecuta pero da un resultado ruidoso que no se parece en nada a una solución para una ecuación de Poisson. Parece tener algo que ver con los espacios de funciones combinadas, pero no puedo encontrar el error.
Agradecería cualquier ayuda o consejos en la dirección correcta. ¡Muchas gracias ya!
Saludos
Markus

Markus
fuente
Permítanme aclarar esto: ¿tienen datos tanto de Dirichlet como de Neumann, pero solo en parte del límite?
Christian Clason
1
Como he entendido el OP, es el gradiente que se da en el límite. Los datos de Dirichlet se utilizan para imponer la derivada tangencial. Pensé que era extraño imponer tanto a Dirichlet como a Neumann en una parte del límite, pero tal vez en esta situación particular sea consistente. Entonces, el problema es más bien cómo aplicar datos de gradiente en el límite (a través de multiplicadores).
enero
Es cierto que eso daría datos consistentes, pero aún tiene el problema de la falta de estabilidad y el hecho de que tiene condiciones de límite en solo una parte del límite.
Christian Clason
Ok, déjame darte más información sobre el problema físico específico que estoy tratando de resolver. Tengo un campo magnético estático que razonablemente puedo suponer que es simétrico rotacionalmente, por lo tanto, 2D. Mido componentes radiales y axiales del vector de densidad de campo magnético a lo largo de una línea, bastante cerca del eje de rotación y me gustaría ver este campo magnético a una distancia considerable de este eje de rotación. La combinación de Dirichlet y Neumann BC fue solo mi idea de abordar el problema como describió Jan elocuentemente: datos de gradiente impuestos en el límite.
Markus
1
OK, eso cambia las cosas significativamente. ¿Entonces tiene un dominio ilimitado e información derivada en toda la parte "finita" del límite?
Christian Clason

Respuestas:

8

Primero, un punto general: no puede prescribir condiciones de límite arbitrarias para un operador diferencial parcial y esperar que la ecuación diferencial parcial (que siempre incluye tanto las condiciones de operador como las de límite) esté bien planteada, es decir, admite una solución única que depende continuamente de datos: todo lo cual es una condición necesaria para realmente tratar de calcular algo.

Dependiendo del operador, a menudo hay una serie de condiciones válidas que puede imponer (para probar, vea la monografía de tres volúmenes de Lions and Magenes). Sin embargo, lo que está tratando de hacer (especifique el gradiente completo, que es equivalente a las condiciones de Dirichlet y Neumann en el mismo límite (parte del) para un PDE elíptico de segundo orden) no está entre ellos, esto se conoce como un problema de Cauchy lateral, y está mal planteado: no hay garantía de que un par de datos límite admitan una solución, e incluso si existe, no hay estabilidad con respecto a pequeñas perturbaciones en los datos. (De hecho, este es el problema original mal planteado en el sentido de Hadamard, y el ejemplo clásico de por qué no se pueden ignorar las condiciones límite cuando se habla de una buena postura. Puede encontrar un ejemplo explícito en sus Conferencias sobre el problema de Cauchy en diferencial lineal parcial ecuaciones de la década de 1920).

(r,R)×(a,b)x=rRxy=ay=b

  1. Si puede imponer condiciones de contorno (Neumann, Robin, Dirichlet, que necesitaría para fijar la constante en la integración de la derivada tangencial, por cierto), entonces es suficiente usar los componentes normales de su gradiente como condición de Neumann (si puede arreglar el modo constante) o integrar los componentes tangenciales como una condición de Dirichlet. Dado que ambas condiciones corresponden presumiblemente a la misma función, obtendrá la misma solución de cualquier manera.

  2. y=ay=bΔu=fΔu+εΔ2u=fε>0H2uuεuε0

    H2

Christian Clason
fuente
Para la implementación por elementos mixtos en FEniCS ver biharmonicdemo. Esto probablemente no tenga el término de Laplace, pero supongo que se puede agregar fácilmente.
Jan Blechta
Hola Christian, gracias por tu sugerencia! Tenía la impresión de que la ecuación de Poisson era benigna en lo que respecta a la estabilidad numérica, gracias por señalarlo. Lo leeré como sugirió. No estoy seguro si esto cambia las cosas sustancialmente, pero como se menciona en el comentario más adelante, la combinación Dirichlet-Neumann es quizás engañosa. 'Todo' que estoy buscando es una forma de imponer datos de gradiente en el límite.
Markus
2
La ecuación de Poisson es benigna, pero esa no es la ecuación que está tratando de resolver :) (Las condiciones de contorno son una parte integral de la ecuación; el operador por sí solo es insuficiente para decidir la estabilidad.)
Christian Clason
Muy bien, eso me da algo para masticar. Gracias a todos por su tiempo, consejos y paciencia, y mis disculpas por caer en la trampa XY ...
Markus
4

No puede esperar que la solución a su problema alterado sea una solución al problema de Poisson porque necesita cambiar el problema de alguna manera para que esté bien planteado.

F(u,λ)=Ω12|u|2dxΩfudxΓNgudS+ΓNλ(uuD)dS
(u,λ)V×L2(ΓN)V={vH1;v|ΓD=0}ΓDuDΓNF(u)
0=DF(u)[v]=ΩuvdxΩfvdxΓNgvdSvV,
ΓNΓD
0=DF(u,λ)[v,μ]=ΩuvdxΩfvdxΓNgvdS+ΓNλvdS+ΓNμ(uuD)dS(v,μ)V×L2(ΓN),
Δu=fun=gλΓNΓN

λ|g|

ΓDvVΓD

La conclusión es que no puede esperar que el PDE de segundo orden admita dos condiciones límite independientes.


F(u,λ)DF(u,λ)[v,μ]derivative()

F(u,λ)λL2(ΓN)λL2(Ω)λR

Jan Blechta
fuente
2ufuH1
Lamentablemente, mis matemáticas no están a la altura y no estoy seguro de las implicaciones matemáticas de los espacios de Banach, pero me cuesta ver por qué la ecuación no es una solución a una ecuación de Poisson cuando el término multiplicador de Lagrange desaparece. Desde el punto de vista físico, una solución (al problema práctico que describí, no me refiero a la solución en el sentido matemático) debe existir hasta donde puedo ver
Markus
Entonces, es más bien un problema inverso, encontrar la condición de límite para el campo lejano que, junto con la condición de Dirichlet que puede imponer, produce el gradiente normal observado en el límite en el que mide.
Markus
3

Su enfoque no puede funcionar, definitivamente debido a la implementación y probablemente debido a su formulación.

La imposición de condiciones de Dirichlet en dolfin, eventualmente establece los DOF ​​correspondientes de su espacio de prueba a cero.

Este es un extracto del manual de fenics :

Capítulo 3.3.9 (final): la aplicación de una condición de límite de Dirichlet a un sistema lineal identificará todos los grados de libertad que deben establecerse en el valor dado y modificará el sistema lineal de modo que su solución respete la condición de límite. Esto se logra poniendo a cero e insertando 1 en la diagonal de las filas de la matriz correspondiente a los valores de Dirichlet, e insertando el valor de Dirichlet en la entrada correspondiente del vector del lado derecho [...]

vΓm

En resumen, al usar la rutina predeterminada en dolfin no puede aplicar tanto Dirichlet como otras condiciones en el mismo límite.

Sin embargo, antes de intentar solucionar esto en su implementación, busque los espacios de prueba correctos para su formulación matemática (como acaba de mencionar @Jan Blechta).

ene
fuente
Comprendo su punto: creo que mi formulación podría no reflejar exactamente lo que he implementado: mis disculpas. El principio variacional no es más que un recuerdo borroso y estoy tratando de entenderlo nuevamente. He leído el manual hacia adelante y hacia atrás junto con algunos ejemplos de código FEniCS que implementa multiplicadores de Lagrange. Pensé que el problema que plantea es exactamente la razón por la que usaría una segunda función de prueba 'd'.
Markus
Estoy de acuerdo con @JanBlechta. Al principio, debes encontrar el espacio correcto para el multiplicador, que no es trivial. Tal vez los textos sobre optimización de restricciones PDE, donde uno usa multiplicadores para incorporar condiciones secundarias, brinden algunas ideas útiles. En este documento , se utiliza un multiplicador para tener en cuenta las condiciones de Dirichlet dependientes del tiempo.
enero