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

Mass conservation workaround #153

Merged
merged 7 commits into from
Feb 23, 2024

Conversation

emmetfrancis
Copy link
Collaborator

For "surface_to_volume" and "volume_to_surface" reactions, volume species are now interpolated onto the surface function space at the start of each Newton iteration. This adds quite a bit to computation time, but does remedy mass conservation issues. Fixes to integration over Meshviews in dolfin should render this unnecessary. (see https://bitbucket.org/fenics-project/dolfin/issues/1139/integral-over-single-facet-in-meshview)

@emmetfrancis
Copy link
Collaborator Author

FYI @jorgensd @finsberg I updated this PR to fix an issue when curvatures were defined as None and introduced the adaptive time stepping function to the model module.

Comment on lines 1423 to 1430
h = h5py.File(str(filename.with_suffix(".h5")), "r")
all_keys = list(h.keys())
if len(all_keys) > 3:
for key in all_keys[3:]:
cur_dim = int(key[-1])
mf_cur = d.MeshFunction("size_t", mesh, cur_dim)
hdf5.read(mf_cur, f"/{key}")
subdomains.append(mf_cur)
Copy link
Collaborator

Choose a reason for hiding this comment

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

This looks a little bit dangerous. First of all, you are opening the same file with h5py that is already open with dolfin.HDFFile. Second, here we assume that mesh functions are saved in a specific order and that what is saved is a MeshFunction.

Also, bringing in a new dependency (h5py) which is know to be tricky to install.

I am wondering if is is better if the user can provide these keys as an argument when loading the mesh, e.g

def load_mesh(
    filename: Union[pathlib.Path, str],
    comm: MPI.Intracomm = MPI.COMM_WORLD,
    mesh: Optional[d.Mesh] = None,
    extra_keys: Optional[List[str]] = None,
) -> LoadedMesh:
    if extra_keys is None:
        extra_keys = []

and then just reading the meshfunctions from those keys.

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, this would be much better, I'll make this change today, thanks!

smart/model.py Outdated
@@ -1158,13 +1164,15 @@ def initialize_discrete_variational_problem_and_solver(self):
# if use snes
if self.config.solver["use_snes"]:
logger.debug("Using SNES solver", extra=dict(format_type="log"))
# Does passing the entire model object to smartSNESProblem slow down execution?
Copy link
Collaborator

Choose a reason for hiding this comment

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

I think one potential problem here is that you are passing in self as an object and setting it as an attribute onself.problem, meaning that self.problem.model = self (=self.problem.model.problem.model and so on). This is called a cycle, and such an object is never cleaned up which might lead to memory leaks.

@emmetfrancis emmetfrancis changed the base branch from development to new-meshview February 22, 2024 02:00
@emmetfrancis emmetfrancis changed the base branch from new-meshview to development February 22, 2024 22:34
@@ -686,34 +686,35 @@
" \"initial_dt\": DT_INITIAL,\n",
" \"adjust_dt\": adjust_dt,\n",
" \"time_precision\": 6,\n",
" \"attempt_timestep_restart_on_divergence\": True\n",
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Note that allowing timestep restart on divergence was required for convergence in this example. We might want to look at this more closely down the line, something about using the 2D mesh with axisymmetry in this case made it more difficult to get this to converge.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Do we need to remove assemble_mixed here? i know that the form isn’t Mixed-dimensional, But it shouldn’t change the result and would introduce less confusion for a user as to when using what assemble command

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Good point, changing that.

Copy link
Collaborator

@jorgensd jorgensd left a comment

Choose a reason for hiding this comment

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

Only a minor comment. Looks good!

@emmetfrancis emmetfrancis merged commit 2f13d2e into development Feb 23, 2024
7 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants