Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

TDycore polymorphism, part 3: separating concerns #215

Merged
merged 75 commits into from
Dec 16, 2021

Conversation

jeff-cohere
Copy link
Collaborator

@jeff-cohere jeff-cohere commented Nov 15, 2021

Continuing in our journey (see #211 and #214):

This pull request defines interfaces for all important operations needed by our different numerical discretizations. Specifically:

  • TDy is now just a container with all the relevant information that describes a problem
  • The TDy struct contains a context pointer that is used to manage the specific implementations using the function pointers defined in the ops table.
  • The context pointer stores the WY, BDM, or MPFA-O implementation selected by TDySetDiscretization.
  • The calculation of boundary conditions (for pressure, temperature, velocity) and sources/sinks (forcings) have been separated into a Conditions type. They have been "vectorized" in the sense that they now can compute their values on multiple points.
  • All model data has been removed from the MaterialProps and CharacteristicCurves types, which now have reworked, vectorized interfaces. This allows different implementations to store material properties differently, instead of forcing them to have similar layouts. The MaterialProps and CharacteristicCurves interfaces now only compute model data.
  • Moving model data out of the above types required some changes in the IO system, in addition to those made in the specific numerical implementations.
  • There's a new EOS type that calculates properties of water.
  • Some source files have been moved into new directories:
    • fe: the finite element discretizations and related utilities live here
    • materials: EOS, CharacteristicCurves, and MaterialProp live here
    • mpfao: the MPFA-O code lives here (not a new folder, but it has taken on new code)
  • Many functions have been removed from the public interface and made internal.
  • The MPFA-O mesh has been consolidated into a smaller number of source files, and it has been separated from TDy itself (because it's part of the MPFA-O implementation).
  • The Fortran interface has been reworked with vectorized functions for material properties, boundary conditions, and forcings, like its C counterpart. All Fortran-related source files now live in src/f90-mod.
  • Several memory leaks in the MPFA-O code have been plugged.
  • Bonus: source files in the library are now built in a deterministic order, which helps when you're trying to fix a lot of things in various places. :-)

In particular, there's a much stronger sense of an "interface" for TDycore. Instead of any piece doing whatever it wants at any point, we're moving toward a set of interfaces that perform specific tasks in specific contexts in a way that suggests a flow of information that is easier to understand.

So where are we?

Everything works.

  • I had to disable a test and some diagnostic output in transient_snes_mpfaof90 because it uses some low-level functions that arguably don't belong in a demo, like directly setting coefficients and extracting data, which implies knowledge of how things are indexed. Maybe we should discuss this when we have a chance.
  • I've disabled the mesh geometry reading/writing capability, since it makes assumptions about the mesh construction process that no longer hold.
  • The TH demos exhibited a very tricky problem that led to the addition of a set_dm_fields operation, required by all models, that set up the field layout on a DM before it's distributed. The way it was done before was not correct, but for some reason, we hadn't experienced any problems with it till now.

What's next?

The bad news is: not everything will be perfect after this PR goes in. A lot of cleanup remains, but much of it is within the separate numerical implementations, which lets different parts of our team adjust things without affecting others.

The good news: the changes in this PR make it easier for us to move in the direction of the DM-oriented problem definition that we started discussing in #211. Briefly: we attach functions to a DM for computing residuals, Jacobians, etc, and then we hand the DM over to the solvers. This allows us to make TDycore's interface smaller and tighter--the dycore produces a DM that gets used by the solvers directly.

Fixes #133
Fixes #196
Fixes #197
Fixes #22

This is a real mess at the moment, but it contains most of the essential
changes. Now to dig our way back to sanity.
Also pausing to see whether I can get rid of an irritating build system
"feature."
@jeff-cohere
Copy link
Collaborator Author

jeff-cohere commented Nov 29, 2021

Okay, I've fixed some issues I found in the higher-level driver interface, which drives the richards and th demos. Now these drivers run, but they get bad answers. I could go ahead and find out exactly what's going on, but it would be faster if someone who knew more about the driver interface could help me understand the assumptions that it makes.

These demos and the transfer of saved geometry info in the transient_snes_mpfaof90 demo are the only issues remaining, as far as I can tell. @bishtgautam , maybe we should discuss how to proceed.

@jeff-cohere
Copy link
Collaborator Author

The Richards driver tests now pass. The TH ones still fail.

Also, I found the source of the issue with reading in mesh geometry from a binary file: the new mesh construction process has already converted indices to compressed row format by the time the data is read, whereas the old construction process read in the geometry before doing that conversion.

@jeff-cohere
Copy link
Collaborator Author

@bishtgautam and I agreed that we can comment out the mesh geometry reading/writing stuff for now, as we may not actually need it.

The TH tests are the only remaining failures.

Comment on lines +326 to +339
! call TDySetPorosityFunction(tdy,PorosityFunctionPFLOTRAN,ierr);
! CHKERRA(ierr);
!
! do c = 1,ncell
! index(c) = c-1;
! call MaterialPropAlpha_PFLOTRAN(alpha(c))
! call MaterialPropM_PFLOTRAN(m(c))
! enddo
!
! call TDySetCharacteristicCurveVanGenuchtenValuesLocal(tdy,ncell,index,m,alpha,ierr)
! CHKERRA(ierr);
!
! call TDySetCharacteristicCurveMualemValuesLocal(tdy,ncell,index,m,ierr)
! CHKERRA(ierr);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a new issue to add this functionality back in the future?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes--see #216

mass_pre = mass_pre + liquid_mass(g)
enddo
! TODO: Mass conservation diagnostic output is disabled for the moment.
!call TDyGetLiquidMassValuesLocal(tdy,nvalues,liquid_mass,ierr)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a new issue to add this functionality back in the future?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not yet. I think we should discuss how we'd like to do this and the similar stuff below. We need to be careful about what we guarantee to be true when we don't necessarily know the discretization (from outside the library). Alternatively, we could have functions that begin with TDyMPFAO for this, which would tie the discretization to the finite volume one, for example.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My point here is that we need to be able to extract the liquid mass, saturation, etc from TDycore. Since that functionality has been disabled in this PR, it would be good to track the issue. Though, the exact approach of how we go about adding this functionality back is still undecided at the moment. We can make a more generic issue of extracting data from the library. Does that make sense?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep. I've added #217

@@ -468,8 +452,8 @@ program main

end do

call TDyGetSaturationValuesLocal(tdy,nvalues,liquid_sat,ierr)
CHKERRA(ierr);
!call TDyGetSaturationValuesLocal(tdy,nvalues,liquid_sat,ierr)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a new issue to add this functionality back in the future?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See #217

Comment on lines +201 to +203
// These don't work, and we'll likely get rid of them.
//PETSC_INTERN PetscErrorCode TDyMeshReadGeometry(TDyMesh*,const char*);
//PETSC_INTERN PetscErrorCode TDyMeshWriteGeometry(TDyMesh*,const char*);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a new issue to add this functionality back in the future?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It depends on whether we need it. Is this just for saving time in debugging? At this point, we could make a proper effort to optimize the calculation of geometry so we don't need this feature.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See #218

@codecov-commenter
Copy link

codecov-commenter commented Dec 16, 2021

Codecov Report

Merging #215 (fe06c7a) into master (b97463a) will increase coverage by 2.14%.
The diff coverage is 67.55%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #215      +/-   ##
==========================================
+ Coverage   49.48%   51.63%   +2.14%     
==========================================
  Files           4        4              
  Lines         685      765      +80     
==========================================
+ Hits          339      395      +56     
- Misses        346      370      +24     
Impacted Files Coverage Δ
demo/transient/transient.c 0.00% <0.00%> (ø)
demo/steady/steady.c 54.13% <72.20%> (+2.72%) ⬆️
demo/th/th_driver.c 100.00% <100.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b97463a...fe06c7a. Read the comment docs.

@jeff-cohere jeff-cohere merged commit f6d7fd0 into master Dec 16, 2021
@jeff-cohere jeff-cohere deleted the jeff-cohere/impls-refactor branch December 16, 2021 17:33
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
cleanup Refactor code enhancement New feature or request MPFAO
Projects
None yet
3 participants