Adding pure periodic boundary conditions and also a generic boundary condition #99
Conversation
These are gradient, divergence, and interpolation periodic operators. They satisfy the following properties: grad = - div^T and interp from nodes to centers is the transpose of the interp from centers to nodes
addBC allowed periodic BC imposing the conditions explicitly. Per files solve periodic BC without imposing any conditions. The domain does not include boundaries. Gral combines traditional operators and handling periodic boundary conditions without specifying explicitly them
jbrzensk
left a comment
jbrzensk
left a comment
There was a problem hiding this comment.
Thanks @mdumett for these additions.
I had a few items with the license info, and location, but I have suggestions for the overall structure.
We do not want to have users interacting with this many functions. It is confusing. I am confused. I think these should be called internally from the addBC functions you made earlier.
Those functions already take parameters for boundary conditions, periodic included correct? If the user requests periodic, then the addBC function should call the appropriate addGralBC1D operator, to generate the periodic version.
Same for the divPer and gradPer. If the user requests periodic everywhere, then the addBC should call divPer or gradPer, and return them to the user as D or G. The user does not need to see the naming used.
I also think, with all of your excellent additions, we need to have subfolders for the operators. We may have to accept there will be $PATH arguments, but there are a lot of files for humans to parse now. I propose a three tiered structure of:
src/MATLAB/
│── +Boundary_Conditions/ % Boundary Conditions
│ ├── addBC1D.m % Main function to apply BCs (user-facing)
│ ├── addBC2D.m %
│ ├── BC_Helpers/ % All internal functions the user does not need to see
│ │ ├── addGralBC1D.m
│ │ ├── addGralBC2D.m
│ │ └── others....m
│── +Interpolators/ % Interpolators
│ ├── intperCenterToNode.m % Main function to apply BCs (user-facing)
│ ├── interpNodeToCenter.m %
│ ├── Interp_Helpers/ % All internal functions the user does not need to see
│ │ ├── interpolCentersToFacesD1DPer.m
│ │ ├── interpolCentersToFacesD2DPer.m
│ │ └── others....m
│── +Operators/ % Mimetic Operators
│ ├── div.m % overloaded functions (user-facing)
│ ├── grad.m % figure out what to call based on number of arguments
│ ├── Operator_Helpers/ % All internal functions the user does not need to see
│ │ ├── DI1.m
│ │ ├── DI2.m
│ │ ├── Div2D.m
│ │ ├── Div3D.m
│ │ └── curvGrad.m
Now, I get that you are not responsible for this, but I think the boundary conditions (which you are now king of 😜 )can start us off. I ask for @jcorbino and @cpaolini @pritkc @JananiPSrinivasan for comments about the structure proposed. I threw this out there without too much thinking.
src/matlab/div2DPer.m
Outdated
| % dx : Step size along x-axis | ||
| % n : Number of cells along y-axis | ||
| % dy : Step size along y-axis | ||
|
|
There was a problem hiding this comment.
We put the license information at the end of the description.
% ----------------------------------------------------------------------------
% SPDX-License-Identifier: GPL-3.0-or-later
% © 2008-2024 San Diego State University Research Foundation (SDSURF).
% See LICENSE file or https://www.gnu.org/licenses/gpl-3.0.html for details.
% ----------------------------------------------------------------------------
src/matlab/div2DPer.m
Outdated
|
|
||
| function D = div2DPer(k, m, dx, n, dy) | ||
| % Returns a two-dimensional mimetic divergence operator | ||
| % when the boundary condition is periodic |
There was a problem hiding this comment.
Maybe add t a qualifier for clarity:
when all of the boundaries are periodic
There was a problem hiding this comment.
Thank you for the review. I will add more comments explaining the specifics of each file. I agree with the fact of adding more options for boundary conditions could be confusing for the user and that most of the files should be hidden.
The only difference with your reasoning is that addGralBC is the version that is most general and has a better implementation than the others. So, addGralBC should be kept as the only public code for handling BC and all the others (including robinBC should be concealed). This also will require modifying all current examples to use addGralBC). If there is an agreement on this, I volunteer to do all these changes in my next version of the code (as soon as possible). Please, let me know.
There was a problem hiding this comment.
In terms of naming, one could rename addBC as addNonPerBC (and call it when BC along certain direction is non-periodic, as addGralBC does) and rename addGralBC as addBC.
There was a problem hiding this comment.
Your addGralBC1D should be named addBC1D. It does all that addBC1D does and more. Is there an issue if your current "addGralBC" boundary condition files replace the existing "addBC" files?
I see the point of the other files, and there needs to be reorganization. The examples are needed as test cases for the periodic operators.
There was a problem hiding this comment.
Currently, the addGralBC calls both addBC and the pure periodic case. Just renaming addGralBC as addBC will overwrite a large part of the code. I can reorganize following the introduction of directories approach you mentioned, and editing the examples to just use addGralBC (which will be named addBC).
There was a problem hiding this comment.
addGradlBC1D only calls addGralBC1Dlhs and addGralBC1Drhs. It seems like it is the same code, just with an initial check to see if it is periodic in addGralBC1Dlhs and addGralBC1Drhs.
If you think of it as a change of addBC1D, you are only changing a few lines, not rewriting tons of code.
|
You are right. addGralBC is independent of addBC. It is just manner of renaming these operators. The dependence is on the gradGral, divGral and lapGral operators. They depend on the original grad, div and lap as well as on the pure periodic operators. The same is true for the interpolators. I can do all modifications and review the examples since addBC treats explicitly periodic BC and addGralBC only implicitly. If you agree, I can proceed doing all changes. |
|
For clarity: All boundary condition operators will be in a subfolder called You keep the examples because they also act as tests for the Gral functions. (Also, why 'Gral'? ) You are keeping the divPer, gradPer, and lapPer as is until we come up with a better name and way of calling the functions. They are needed, and there is no way around it. I would like them renamed better, maybe "div1DPeriodic". Same for the interpolators. They will go the same way as the boundary conditions, where a few will be used to call the more complicated cases. interpolFacesToCentersG1DPeriodic.m Examples will be combined in a future PR. There are just too many to parse now. Interpolators will be combined in a future PR. Smae reason. I think to really make the boundary condition functions universal, will be for addB1D to detect what is passed, either the grad, div, or lap, based on the size, and to modify the operator appropriately. |
|
Agree with all but except with "... to really make the boundary condition functions universal, will be for addB1D to detect what is passed, either the grad, div, or lap, based on the size, and to modify the operator appropriately." All BC functionalities (addGralBC, addBC, and pure periodic case) consider general second-order PDE operators (i.e., the PDE can be composed of Laplacians, gradients, divergences, curls, etc., provided the whole PDE operator is of second-order) and hence the BC functions are not associated to a Laplacian or a gradient or a divergence per se. |
|
Wait, these addBCs only work with Laplacian minetic operators? |
|
No. addGralBC and addBC work for any PDE of second-order. In general the PDE can have several terms, each of them having any operator provided the order of the PDE is two. |
|
There is an assert that the A matrix be square
…On Fri, Mar 21, 2025 at 9:03 AM mdumett ***@***.***> wrote:
No. addGralBC and addBC work for any PDE of second-order. In general the
PDE can have several terms, each of them having any operator provided the
order of the PDE is two.
—
Reply to this email directly, view it on GitHub
<#99 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AA7CSHKPFWR5YEQMFTN6VCD2VQZ65AVCNFSM6AAAAABZGWAN26VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDONBTHAYTEMJZGQ>
.
You are receiving this because you were assigned.Message ID:
***@***.***>
[image: mdumett]*mdumett* left a comment (csrc-sdsu/mole#99)
<#99 (comment)>
No. addGralBC and addBC work for any PDE of second-order. In general the
PDE can have several terms, each of them having any operator provided the
order of the PDE is two.
—
Reply to this email directly, view it on GitHub
<#99 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/AA7CSHKPFWR5YEQMFTN6VCD2VQZ65AVCNFSM6AAAAABZGWAN26VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDONBTHAYTEMJZGQ>
.
You are receiving this because you were assigned.Message ID:
***@***.***>
|
|
That is true. It is expected that the PDE is written in terms of cell centers both input and output (hence the square matrix). |
|
Are the periodic boundary conditions required? Or is it required by this pull request? They may be better suited to a separate pull request. To save sanity in testing everything. |
|
The periodic boundary conditions are required. I suggests to add the general boundary conditions (addGralBC) renamed as addBC. There will be only nine BC files. In addition, the additional versions of the gradient, divergence and curl have to be added in this PR (PR1). In the next PR (PR2), we can deprecate robinBC and all others BC operators but addBC, and rewrite all examples to use exclusively addBC. Later, in another PR (PR3), we can organize implement the new directory organization. Another sequence of work could be PR1, PR3, PR2. Let me know what you think. |
|
My mistake! I agree with the PR order.
|
Some quick suggestions: The same applies to |
|
@jbrzensk & @aboada : Mimetic differences use staggered grid. Periodic interpolators are needed together with periodic gradient and divergence. According to your item list, |
Sounds good to me. Thanks! |
|
@jbrzensk and @aboada: With respect to the name of the new function that adds boundary conditions to the PDE and rhs operators (addBC), I realized that all the boundary condition operators that are considered for MOLE are just for scalar fields. Vector fields have different Dirichlet and Neumann BCs (and possibly Robin as well as others). For example, if u is a vector field, n the unit exterior normal, and t the unit tangent, u.n = 0 is Dirichlet for vector fields and is not considered currently in MOLE. In addition, [u]_t = 0 (no relative motion between a viscous fluid and a rigid wall), [S n]_t, + b [u]_t = 0, with S = mu (div_x u + div_x^T u), the viscous stress tensor, as well as others are frequent BCs for vector fields and MOLE does not contain any of those. I suggest using the name addScalarFieldBC instead of addBC for the new functionality. What do you think? |
|
So the boundary conditions you have made up are for scalars only? If so, you should call your boundary condition functions "addScalarBC", which will take the place of "addBC", with an eye on the future where you will make a "addVectorBC" function. |
Gral is changed by Scalar, Per by Periodic and added License Identifier. However, there are divGral, gradGral and lapGral operator that still need to be resolved
|
@jbrzensk @aboada: As you can see in the addBCs3 branch, I have changed the name addGralBc by addScalarBC, the name Per by Periodic and added the License-Identifier in all source files. Also modified the examples related to these changes. However, there remain files divGral, gradGral and lapGral. The reason is that the Gral versions call the Periodic operators or the old ones (div, grad, lap) according to the BC type in each coordinate. To handle this the Gral versions have a couple of extra arguments associated to BCs. We need to decide how to incorporate this into MOLE carefully since it may create a legacy issue. It might make incompatible users old codes. I propose to rename the old versions as divNonPeriodic, gradNonPeriodic, lapNonPeriodic and the Gral versions as div, grad, lap. If the last two arguments of the Gral versions are empty, it will call the NonPeriodic versions (i.e, the old ones). If they are non-empty it will choose between Periodic or NonPeriodic operators according to the BCs in each coordinate. I believe this is a simple way to solve this issue. Please, let me know what you think. |
|
This seems needlesssly overly complicated. What does "Gral" mean? I'm going to offer suggestions to simplify this large code change. Would it be easier to modify the existing div, grad, etc, to have the option of taking two extra arguments? Then call divGralXX to decide how to do periodic? That would allow existing code to still work, and allow people to use the familiar OR, just call your grad3DPeriodic is a repeat of the command gradGral3D(k,m,dx,n,dy,o,dz,[0,0],[0,0],[0,0]), I don't think we need the individual grad3DPeriodic, div3DPeriodic, etc, if they are just wrappers to the equivalent grad and div calls. (Also
Here is an example of the first thing I was talking about. function D = div(k, m, dx, dc ,nc)
% Returns a m+2 by m+1 one-dimensional mimetic divergence operator
%
% Parameters:
% k : Order of accuracy
% m : Number of cells
% dx : Step size
% ( optional ) dc : (a0) vector for left/right BC
% ( optional ) nc : (b0) vector for left/right coefficient
% Assertions:
assert(k >= 2, 'k >= 2');
assert(mod(k, 2) == 0, 'k % 2 = 0');
assert(m >= 2*k+1, ['m >= ' num2str(2*k+1) ' for k = ' num2str(k)]);
assert( nargin == 3 || nargin == 5 , 'Need 3 or 5 arguments');
if ( nargin == 5 )
D = divGral1D(k,m,dx,dc,nc);
else
D = sparse(m+2, m+1);
switch k
case 2
for i = 2:m+1
D(i, i-1:i) = [-1 1];
end
case 4
A = [-11/12 17/24 3/8 -5/24 1/24];
D(2, 1:5) = A;
D(m+1, m-3:end) = -fliplr(A);
for i = 3:m
D(i, i-2:i+1) = [1/24 -9/8 9/8 -1/24];
end
case 6
A = [-1627/1920 211/640 59/48 -235/192 91/128 -443/1920 31/960; ...
31/960 -687/640 129/128 19/192 -3/32 21/640 -3/640];
D(2:3, 1:7) = A;
D(m:m+1, m-5:end) = -rot90(A,2);
for i = 4:m-1
D(i, i-3:i+2) = [-3/640 25/384 -75/64 75/64 -25/384 3/640];
end
case 8
A = [-1423/1792 -491/7168 7753/3072 -18509/5120 3535/1024 -2279/1024 953/1024 -1637/7168 2689/107520; ...
2689/107520 -36527/35840 4259/5120 6497/15360 -475/1024 1541/5120 -639/5120 1087/35840 -59/17920; ...
-59/17920 1175/21504 -1165/1024 1135/1024 25/3072 -251/5120 25/1024 -45/7168 5/7168];
D(2:4, 1:9) = A;
D(m-1:m+1, m-7:end) = -rot90(A,2);
for i = 5:m-2
D(i, i-4:i+3) = [5/7168 -49/5120 245/3072 -1225/1024 1225/1024 -245/3072 49/5120 -5/7168];
end
end
D = (1/dx).*D;
end
end |
|
@jbrzensk @aboada: There is a lack of understanding of the proposed changes. There are two different parts in them. The operators and the boundary conditions. Operators: I refer to grad as an example but what is mentioned applies to div too The new grad matrix representation depends on whether the BC along each axis is periodic or not. The changes proposed to add a vargin argument for grad. Whether this argument its empty or not, in each direction, will be used to call the periodic or the old ones. The user will not have to change anything in their codes, since it is a vargin argument. I believe, this is in essence what you propose. Boundary conditions: The implicit code will be used for the Plasma Fusion project at least. I hope this helps to understand the new code. Please, let me know if it is clear now. Thank you for pointing out some places were I forgot to substitute Per by Periodic. I would like to continue modifying the examples that use addBC to use addScalarBC instead and then reorganize the directories. |
|
"Operators: I refer to grad as an example but what is mentioned applies to div too What is the difference between these two commands? G = grad3DPeriodic(4,10,1,11,1,12,1);
G = gradGral3D(4,10,1,11,1,12,1,[0;0;0;0;0;0],[0;0;0;0;0;0]);Matlab says they are the same: >> Gp = grad3DPeriodic(4,10,1,11,1,12,1);
>> G = gradGral3D(4,10,1,11,1,12,1,[0;0;0;0;0;0],[0;0;0;0;0;0]);
>> Gp = full(Gp);
>> G = full(G);
>> diff = sum( abs( Gp(:) - G(:) ) );
>> diff
diff =
0This is where my comment came from. These appear to be the same. If they are, we do not need another function to generate the 2D and 3D versions of a periodic operator. The 1D version IS needed, but the others are superfluous. "Boundary conditions: The "AddScalarBC" seems to be the same as "addBC", except with the additional property that "addScalarBC" can "separate out periodic versus non-periodic boundary conditions". This implies this is exactly the same as the existing addBC functions, with an added internal feature. If so, why not just add the feature to the existing addBC functions? You do not need 9 NEW functions. |
|
Operators: Boundary conditions: |
I would extract the else statement in a separate function; the cognitive complexity is too high. Can we make these functions private (like divGral1D)? Does it make sense? It's been a while since I have coded in Matlab but I think we can benefit from reducing the scope (while still having modularity). |
Typos in grad2DPeriodic and grad3DPeriodic fixed bfore removing those files. Just a sanity check. In addition, removing repeated rows in addBC2D and addScalar2D to speed up these codes. addBC2D will be removed. Again, having a better code before removing it just in case we need to get back to it in the future.
div, grad, and lap extended for non-periodic and periodic bc keeping legacy code. All examples in the MOLE were tested successfully
|
This code has extended div, grad, and lap to include periodic BCs. With respect to the conflicts:
This code can safely be merged with MOLE. The only thing remaining is the organization of the different directories. |
jbrzensk
left a comment
jbrzensk
left a comment
There was a problem hiding this comment.
I noted the comments should show that the last arguments are optional, and that there should be error checking for the number of arguments.
I checked the code out for Octave compatibility, as well, and that the examples still work.
The word optional for arguments are not necessary have been included. Explicit constraints to the number of arguments of div, grad, and lap were added.
|
Optional arguments have been made explicit. Constraints to the number of arguments of div, grad, and lap are included. |
Fixing error message typos.
|
Thank you for finding the typos. I have fixed them. |
|
@mdumett Now you need to fix the conflicts before this branch can be merged. |
|
Even though examples do not use any longer addBC, maybe it is worth to have that working code in the library. It imposes periodic boundary conditions explicitly and not implicitly as the new code. It might be useful for developers in the future, even though users will not have access to it. What you think? I will still need to resolve the operator conflicts. |
|
I do not think it is necessary. If someone really wants to see how it could be done, they could go back through the github history and see how it was done. |
have updated the pure periodic boundary condition code (Per) with a new one called addGralBC which actually combines the capabilities of addBC and the pure periodic boundary conditions. The matrix sizes are reduced because of the chopping of boundaries along directions with periodic boundary conditions.
If you recall,
Examples of how to use it are added.
So, now we have two ways of imposing boundary conditions: explicitly (using addBC) and implicitly (using addGralBC, which calls Per for the directions where there is a periodic boundary condition).
addGralBC and Per are in my addBCs3 branch. Please, review them to be able to push it as soon as possible.
This ends all boundary treatment.