Algorithm for P and Q weights

[mdumett]

Is your feature request related to a problem? Please describe.
A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]
Even it is not required for common PDE solution, having P and Q weights will simplify the computation of integrals with a mimetic approach and will help for theoretical proofs. For example, the computation of mass, momenta, energy, etc., will leverage on these weights.

Describe the solution you'd like
A clear and concise description of what you want to happen.
We should implement 1D, 2D, 3D versions of them in C++, each on its own file.
In addition, MATLab/Octave 2D versions of them have wrong sizes, and the 3D versions of them are missing in MATLAB/Octave.

Describe alternatives you've considered
A clear and concise description of any alternative solutions or features you've considered.

I have just noticed that the orders of the MATLAB/Octave 2D weightsP and weightsQ matrices have wrong sizes.
This should be also corrected in this issue.

In addition, Order 8 of the algorithm requires a deeper explanation based on a Linear Programming Problem Simplex Algorithm, which will be suggested in a different issue. The following is for orders of accuracy k = 2,4,6.

a) C++ 1D weights P(order k, x-dir cells m, x-dir cell size dx):
G = grad(k, m, dx);
b = [-1; zeros(m, 1); 1];
P = G'\b;

b) C++ 2D weights P(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy)
Use 2D formula of section 8.4, page 26, of the attached document for P_{xy}^{(k)}

c) C++ 3D weights P(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy, z-dir cells o, z-dir cell size dz)
Use 3D formula of section 8.4, page 26, of the attached document for P_{xyz}^{(k)}

d) MATLAB/Octave 2D weights P(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy)
Use 2D formula of section 8.4, page 26, of the attached document for P_{xy}^{(k)}

e) MATLAB/Octave 3D weights P(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy, z-dir cells o, z-dir cell size dz)
Use 3D formula of section 8.4, page 26, of the attached document for P_{xyz}^{(k)}

f) C++ 1D weights Q(order k, x-dir cells m, x-dir cell size dx), (for orders 2,4,6):
D = div(k, m, dx);
b = [-1; zeros(m-1, 1); 1]; % RHS
Q = [1; D(2:end-1, :)'\b; 1];

g) C++ 2D weights Q(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy)
Use 2D formula of section 8.2, page 24, of the attached document for Q_{xy}^{(k)}

h) C++ 3D weights Q(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy, z-dir cells o, z-dir cell size dz)
Use 3D formula of section 8.2, page 24, of the attached document for Q_{xyz}^{(k)}

i) MATLAB/Octave 2D weights Q(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy)
Use 2D formula of section 8.2, page 24, of the attached document for Q_{xy}^{(k)}

j) MATLAB/Octave 3D weights Q(order k, x-dir cells m, x-dir cell size dx, y-dir cells n, y-dir cell size dy, z-dir cells o, z-dir cell size dz)
Use 3D formula of section 8.2, page 24, of the attached document for Q_{xyz}^{(k)}

CSRCR2023-01.pdf

Additional context
Add any other context or screenshots about the feature request here.

[joehellmersNOAA]

@mdumett and @jcastillo001 I've got this working in my fork. I just need to refine my test a bit then I can create a pull request.

After than I'll create an example or two.

FYI, in order to use Armadillo/SuperLU I needed to add a utility function to solve using QR decomposition, since the system to be solved did not result in a square matrix.

[mdumett]

Joe, it is true, the system for finding P and Q is not a square system. In addition, that system does not enforce positivity of the weights. It just happens up to order six but not for order eight. A way to overcome the fact that positivity is not enforced is to construct a linear programming problem (LPP) where the linear system for P and Q are introduced as constraints together with a positive constraints for P and Q weights (i.e., p_i > epsilon > 0, q_j > epsilon > 0). This LLP can be easily solved by a Simplex algorithm. I have an implementation of this latter algorithm.

[jcastillo001]

see atached paper Prof. Jose E. Castillo, Ph.D. Director/Professor Computational Science Research Center San Diego State University Google Scholar <https://scholar.google.com/citations?user=6qGWVfkAAAAJ&hl=en> www.csrc.sdsu.edu

…

On Tue, Dec 16, 2025 at 10:58 AM Joe Hellmers ***@***.***> wrote: *joehellmersNOAA* left a comment (csrc-sdsu/mole#224) <#224 (comment)> @mdumett <https://github.com/mdumett> and @jcastillo001 <https://github.com/jcastillo001> I've got this working in my fork. I just need to refine my test a bit then I can create a pull request. After than I'll create an example or two. FYI, in order to use Armadillo/SuperLU I needed to add a utility function to solve using QR decomposition, since the system to be solved did not result in a square matrix. — Reply to this email directly, view it on GitHub <#224 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ATWZMCP6JYDWWWPKJ23PDY34CAM2PAVCNFSM6AAAAACJFRB4ZWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTMNRQHE4DENRUGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[joehellmersNOAA]

Thanks @mdumett. The C++ version of MOLE's divergence operator doesn't support more than order 6 currently.

test6: mole/src/cpp/divergence.cpp:21: Divergence::Divergence(arma::u16, arma::u32, Real): Assertion k > 1 && k < 7 failed.

I'm not sure why this limitation is present. Maybe @jcorbino or @jbrzensk knows.

It looks like the gradient operator DOES support 8th order.

[jbrzensk]

We stopped at 6 because no one uses 8. Most schemes are second- or fourth-order. There are incredibly few things that need sixth-order operators. Order 8 has stability issues for most things of interest.

[mdumett]

The Matlab version of gradient, divergence, interpolations, and Laplacian support up to eight order. There are no stability issues in those operators. What happens is that the algorithm proposed for getting P and Q coefficients do not enforce positivity of them. By just pure luck up to six order are positive. In Matlab we could add the coefficients obtained by the Linear Programming approach mentioned. Then those coefficients can be added to the C++ Fortran, Python and Julia versions.

[aboada]

We should implement 1D, 2D, 3D versions of them in C++, each on its own file.

@mdumett , why do we need to put the 1D, 2D, 3D versions in different files? That will go against the design pattern we have in main/src/cpp.

[mdumett]

When I propose 1D, 2D and 3D versions in different files I was not aware of the C++ implementation. Of course, by consistency they should be in the same file. But that is not the main issue here. The main issue is to have coefficients up to eight order to be consistent with the Matlab side.

@jbrzensk: By the way there is an eight order implementation of the gradient in C++.

[aboada]

When I propose 1D, 2D and 3D versions in different files I was not aware of the C++ implementation. Of course, by consistency they should be in the same file. But that is not the main issue here. The main issue is to have coefficients up to eight order to be consistent with the Matlab side.

@jbrzensk: By the way there is an eight order implementation of the gradient in C++.

My point is that the solution is misleading. Could we update the solution to avoid any undesired changes?

[joehellmersNOAA]

@aboada if you look at the branch dev/cpp-weights you'll see my latest code. There are two new classes, one for P and one for Q weights.

[aboada]

@aboada if you look at the branch dev/cpp-weights you'll see my latest code. There are two new classes, one for P and one for Q weights.

Right, thanks. So far, you are only supporting 1D, right?

[mdumett]

There are formulas for 2D and 3D P and Q but they have not been implemented in any of the versions yet.

[joehellmersNOAA]

@aboada The branch I referenced earlier was too far behind, so I created a new one. dev/cpp-weights2.

By the way expanding from the 1D to 2 and 3D weights doesn't look to bad based on what we do in MATLAB/Octave. (which is only 2D for the P weights, the 2D Q weight doesn't look complete at all).

Nonetheless, I think we should do the 1D only for now, and work on the others later.