Broken FEEC projections on polar domains#576
Broken FEEC projections on polar domains#576alisa-kirkinskaia wants to merge 82 commits intodevelfrom
Conversation
…into polar_splines_alisa
- Remove polar_mapping and use_logical_sol options - Change the analytical solution - Avoid Laplacian() for right hand side - Add disk radius as parse argument
Coverage summary from CodacySee diff coverage on Codacy
Coverage variation details
Coverage variation is the difference between the coverage for the head and common ancestor commits of the pull request branch: Diff coverage details
Diff coverage is the percentage of lines that are covered by tests out of the coverable lines that the pull request added or modified: See your quality gate settings Change summary preferences |
yguclu
left a comment
yguclu
left a comment
There was a problem hiding this comment.
Thanks for this nice work!
I have started reviewing this PR and I have a couple of questions about test_conga_projections.py
| # Checking projection property P0(P0(phi)) = P0(phi) | ||
| assert np.allclose(P0.dot(y)[:, :], y[:, :], atol=1e-12, rtol=1e-12) | ||
|
|
||
| sp_P0 = P0.tosparse() |
There was a problem hiding this comment.
The tosparse method of a LinearOperator object does not give any guarantee about the sparse format of the returned object, which only needs to be of type scipy.sparse.spmatrix (or even better, scipy.sparse.sparray).
A few lines below you access the attributes row, col, and data, which implies that sp_P0 is of type scipy.sparse.coo_matrix. In order to make sure that the current line is always correct, we need to modify it as follows:
| sp_P0 = P0.tosparse() | |
| sp_P0 = P0.tosparse().tocoo() |
Since in this case P0.tosparse() already creates a COO matrix, the tocoo() call returns the same object with negligible computational overhead (see here).
FYI, you can check object identity using the is operator,
sp_P0_any = P0.tosparse() # COO format in this case
sp_P0_coo = sp_P0_any.tocoo() # no copy is made
assert sp_P0_any is sp_P0_coo # True| # gather sparse matrix entries (rows, columns, data) on root process | ||
| payload = (sp_P0.row, sp_P0.col, sp_P0.data) | ||
| parts = mpi_comm.gather(payload, root=0) | ||
|
|
||
| if mpi_comm.rank == 0: | ||
|
|
||
| # Concatenate the result of gather (equivalent to summation of local matrices) | ||
| rows = np.concatenate([p[0] for p in parts]) | ||
| cols = np.concatenate([p[1] for p in parts]) | ||
| data = np.concatenate([p[2] for p in parts]) |
There was a problem hiding this comment.
Assuming that mpi_comm contains N+1 processes, here we are taking these steps:
- Gather on process
0the three-tuples(sp_P0.row, sp_P0.col, sp_P0.data)to obtain
parts = [(r0, c0, d0), (r1, c1, d1), ..., (rN, cN, dN)]- Separate
partsback into(rows, cols, data):
rows = [*r0, *r1, ..., *rN]
cols = [*c0, *c1, ..., *cN]
data = [*d0, *d1, ..., *dN]Instead, could we not gather separately rows, cols, and data on process 0? This way we should obtain a NumPy array for each of these without the need to call np.concatenate(). A basic version of the code should read
| # gather sparse matrix entries (rows, columns, data) on root process | |
| payload = (sp_P0.row, sp_P0.col, sp_P0.data) | |
| parts = mpi_comm.gather(payload, root=0) | |
| if mpi_comm.rank == 0: | |
| # Concatenate the result of gather (equivalent to summation of local matrices) | |
| rows = np.concatenate([p[0] for p in parts]) | |
| cols = np.concatenate([p[1] for p in parts]) | |
| data = np.concatenate([p[2] for p in parts]) | |
| # gather sparse matrix entries (rows, columns, data) on root process | |
| rows = mpi_comm.gather(sp_P0.row, root=0) | |
| cols = mpi_comm.gather(sp_P0.col, root=0) | |
| data = mpi_comm.gather(sp_P0.data, root=0) | |
| if mpi_comm.rank == 0: |
but I believe that Gather (capitalized) should be used for increased efficiency.
| P0.dot(x, out=y) | ||
|
|
||
| # Checking projection property P0(P0(phi)) = P0(phi) | ||
| assert np.allclose(P0.dot(y)[:, :], y[:, :], atol=1e-12, rtol=1e-12) |
There was a problem hiding this comment.
I have noticed that you use a hard-coded tolerance of 1e-12 at many places in this file. How about defining the constants RTOL and ATOL just once at the beginning of this file, and then using them everywhere? This way it will be very easy to modify one of the two tolerances in the future, if need be
| assert np.allclose(y_sp, y.toarray(), atol=1e-12, rtol=1e-12) | ||
|
|
||
| if __name__ == '__main__': | ||
| test_PolarProjection_V2(1, [8, 10], [2, 2], False) No newline at end of file |
There was a problem hiding this comment.
Apparently the newline character is missing at the end of this file
| except NameError: | ||
| import matplotlib.pyplot as plt | ||
| plt.show() | ||
| plt.show() No newline at end of file |
There was a problem hiding this comment.
Please revert this whitespace-only change. It removes the required newline character at the end of file
2 participants
This PR introduces broken FEEC projections for C0 and C1 sequences (conga_projections.py) in 2D. The formulas for the projections can be found here on pages 19-23.
Added projections
C0PolarProjection_V0,C0PolarProjection_V1,C0PolarProjection_V2onto C0 sequence, and analogously onto C1 sequence.dotmethods of the projections are implemented for parallel execution for arbitrary domain decompositions in the angular and/or radial dimensions.Added unit tests for the projections checking:
tosparsemethods against reference matrix operatorsdotandtosparsemethodsPoisson example
Example of run:
mpirun -n 2 python poisson_2d.py -S -n 16 24 -d 2 2 -t disk -D 0.2 -m 'C0conga'Exact solution, approximate solution and error plot:
Example of run:
mpirun -n 2 python maxwell_2d.py -S -n 16 20 -d 2 2 -T 1 -D 0.2 -s 1Plot of exact solution and approximate solution at final time T = 1: