This doesnāt seem to me to be a well posed problem. The solution u is overconstrained on Gamma_2, because it is asked to be both equal to V (a Dirichlet condition, which would be acceptable on its own), and such that the average of a partial derivative is equal to one.
The point here is that V is an unknown quantity, which must be solved for using the Integral constraint involving the gamma term.
OK, I didnāt get that from reading the text the first time. I donāt have specific experience for this kind of constraint, but you may want to start by searching the forum with ānullspaceā as keyword, and start reading recent posts about that.
I found the following post that seems to solve a similar type of problem
https://fenicsproject.org/olddocs/dolfin/2018.1.0/python/demos/singular-poisson/demo_singular-poisson.py.html
Unfortunately, the solution is in fenics instead of fenicsx, so I donāt know how to convert it. (I am a complete beginner and trying to teach myself this stuff.) So here is a comment I would like to make about the problem that is outlined, called the āsingular Poisson problemā. If one were to choose any single point on the boundary and impose a Dirichlet condition at that point, the problem would no longer be singular. Now impose the condition that the average of u over the entire domain must be zero, and use that condition to solve for the unknown value in the Dirichlet condition. This is completely analogous to what I am trying to do in my problem. Apparently this can be handled in fenics, but how do we do it in fenicsx?
In the tutorial, they use something called a Krylov solver to impose the average condition. My method above seems more straightforward. The question is whether or not fenicsx is set up to handle my approach.
I just noticed that the singular Poisson problem was implemented in dolfinx in
Not sure if I can figure out how to modify this approach to solve the above problem with a current condition, but I will spend some time thinking about it and report my questions or answers.
Having seen your recent posts in Issues Solving Pure Neumann Problems in dolfinx - #7 by Daniel.Baker, Iād like to remark that one difference between the problem here and the one in the other post is that here you actually have two constraints:
- the integral constraint on Gamma_2
- the fact that the solution u must have the same value on all DOFs on Gamma_2
Having said that, I have no experience in how to include these constraints in dolfinx, and hence I canāt help any further.
I agree entirely with your observation. Having said that, it is perhaps worth being more explicit about the second constraint. Letās choose an arbitrary point x_0 in Gamma_2. Then the second constraint says that the difference u(x)-u(x_0)=0 for every point x in Gamma_2. Thus, if there are N DOFs for u in Gamma_2, there will be N-1 constraints of the second type and the integral constraint is also a single constraint. So the total number of constraints is N, the number of DOFs for u in Gamma2. Thus the problem is well defined.
In Issues Solving Pure Neumann Problems in dolfinx one can also choose an arbitrary point x_0 on the boundary of the region and stipulate a Dirichlet condition at x_0. Everywhere else is Neumann boundary conditions. The value of the Dirichlet condition is then determined by the integral constraint on u over the entire domain. In this regard, the two problems are very similar.
@Daniel.Baker and @dokken may want to continue here the discussion from Issues Solving Pure Neumann Problems in dolfinx - #8 by dokken
I will try to reply to comments from @dokken that seem to be relevant to the problem I posed at the top of this thread. I have looked through his suggestions, and I think that, in principle, the information he provided could be used to solve the above problem, if I only understood how to implement it in code. So here is the reference from @dokken that I am interested in
Specifically, he refers to
In my case, let u_0, u_1, ā¦,u_N be the variables on the boundary Gamma2. These constraints would be
u_i=u_0 for all 0<i<N
There is one more inhomogeneous constraint I need, which, in variational form, could be written as
So I think that all of this would fit into the format that it described. My job is to understand how this code is implemented, for example, in the tutorial
https://jsdokken.com/dolfinx_mpc/python/demos/demo_stokes.html
It would greatly help me to understand what is in this tutorial, if someone could document for me the mathematical PDEās and boundary conditions that are being solved in it. In the tutorial itself, all I see is code, which so far I have been unable to understand. Thanks ahead of time for your help!
This was my bad. Full documentation of the equations is cited in the demo_stokes in the link
https://arxiv.org/pdf/2001.10639.pdf
It might help to put this reference at the top of the demo. The following reference is also helpful
https://hal.science/hal-04077986/document
A presentation of DOLFINx_mpc as a framework was done at FEniCS 21: https://fenics2021.com/slides/dokken.pdf
and at Sorbonne/DāAlembert institute fall of 2023: https://www.youtube.com/watch?v=4LSWOdyiGH8&ab_channel=InstitutJeanLeRonddāAlembert
I have dolfinx up and running. How do I find dolfinx_mpc so I can install it? Sorry to have to keep asking these questions. I should be starting to know this stuff by now.
Got it! Iām up and running.
Depending on how you installed dolfinx, choose the appropriate way:
or
I have some extra questions about dolfinx_mpc.
- Will it handle inhomogeneous constraints? For example, an inhomogeneous periodic boundary condition, where the difference between the left and right boundary is a known function?
- If the known function in question (1) depends on a single unknown constant, can one introduce an additional constraint to determine the value of this constant? For example if the two boundaries are at x=0 and x=1, suppose we want to introduce the condition
u(x=1)=u(x=0)+K
where K is to be determined by \int_Omega u dx =1 (again, an inhomogeneous constraint)?
No, neither of those are currently supported.
-
should be quite easy to add in, But i simply havenāt had the time for it atm.
-
that one i dont know how one would easily implement as it is non-linear, all degrees of freedom are connected through the integral constraint (+ having to solve 1. first)
Regarding 2 above, the system I want to solve is the Laplace equation, so I think determining K via the integral constraint is a linear problem. That said, without the capability to do 1 in particular, I donāt think that the problem I originally posed at the top of this thread is solvable with the current resources in fenicsx. Do you agree? This surprises me because everything is linear with inhomogeneous linear constraints. Was this solvable in the earlier fenics?
Legacy fenics had support for real spaces, Which would solver the int u dx=1 and u(x_i)=u(x_j) if the meshes align such that x_i and x_j has a dof at that location. It would silently fail of this mapping did not exist.
Dolfinx_mpc extended this to handle u_k = sum_i=1,iā k^M alpha_i u_i
which covers a broad range of constraints (But not all constraints).
Note that dolfinx_mpc is open source, and Im happy for others to contribute to it. I have alot of other projects on my plate these days, so adding new functionality to dolfinx_mpc is not on the top of my list of priorities.
You have already been very gracious with your time, and I appreciate your help. As a new user, my focus has been to understand what the current capabilities of FEniCSx are and to see if I can use it for problems of interest to me. If I ever get to the point where I understand the code that is already up and running, I would also be interested in modifying it to meet my own needs and those of others. Right now that point seems to be very far off in the future. Thanks again for your help!
I am coming back to this problem now that I understand FEniCSx a little better. To implement a current condition as above, one must introduce a new variable V, which is the unknown constant voltage imposed on the boundary Gamma2. This voltage V must be determined by solving the constraint that the integral of the normal gradient of u over Gamma2 is required to be equal to 1, as noted in the condition initially posted in the problem. It seems to me that the best way to do this is to rewrite a Newton Solver where I can expand the original matrix A by an extra row and column to include the current contraint as a way to determine V.
There are some good descriptions of how to calculate the matrix A and the vector b, but I am not sure what sort of operations I would need to expand this matrix by adding an additional row and column to solve for V. Are there some examples out there that someone could point me too? So far, I have been looking at
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html
as well as Section 5.3.2 of the book
If needed, I could supply more details about how I want to do this, but perhaps what I have written above is enough for someone to point me to the right material. Thanks in advance!
Not sure if I can be of any help, but your problem looks like what I am doing in certain ways. I solve 2 coupled PDEs, one of them is Laplace equation for the voltage, where I specify a given current (the flux of voltage in some way) passing through 2 boundaries like in your case. As voltage is defined up to an arbitrary constant, I do not impose any Dirichlet b.c. for the voltage on any surface, I let the solver do its job and if itās unable.to converge, I can deal with specifying the null space. The thing is that the function āāvoltāā or electrostatics potential is determined automatically just given the input/output current.
If youāre interested, have a look at the codes I posted in my recent post history.