Adaptive Mesh Refinement #110

hughcars · 2023-10-16T15:41:05Z

Introduces a Solve-Estimate-Mark-Refinement adaptive mesh refinement loop for all problem types aside from Transient.
Allows for conformal rebalancing through a gather scatter model. Supports nonconformal adaptation through patches applied to mfem for higher order Nedelec tets, interior boundaries and coarsening.
Does not currently support WavePorts for nonconformal meshes.
Employs a Dorfler marking scheme (a.k.a bulk chasing) whereby a minimum number of elements are marked that achieves a given fraction of the total global error.

Resolves #2

sebastiangrimberg

I have not yet looked at the larger changes in basesolver.cpp and dorfler.hpp/cpp or geodata.hpp/cpp, but here is a first pass.

Related, but does this comment need addressing for the nonconformal case?

palace/palace/fem/coefficient.cpp

Line 81 in 2b66e53

mesh.GetFaceInfos(i, &info1, &info2); // XX TODO: Nonconforming support

cmake/ExternalMFEM.cmake

sebastiangrimberg · 2023-10-17T19:44:17Z

docs/src/config/model.md

@@ -54,7 +54,19 @@ scale based on the bounding box of the computational domain.
            "Radius": float
        },
        ...
-    ]
+    ],
+    "Adaptation":


I think these should be part of the top level "Refinement" object rather than adding a new sub-object. Also, this is a lot of options, surely some are only for advanced usage? I would remove "WriteSerialMesh" at least, for example. Likewise I'm not sure when SaveStep would not just be 1 so I'd move that to advanced too.

Part of the rationale for the substruct was that I had been thinking of Adaptation as a subtype of Refinement. If the idea is that Refinement becomes Adaptation however, then this could change. A lot of these options would then need to be prefixed with Adaptive or similar. I thought the substruct was more encapsulating though, hence this choice.

The main advantage of moving to Advanced options is not having to write a doc string for an option, but given I have I'm not sure the extra information is worth removing. Independent of that issue though:

I could definitely see "SaveStep" equal to 0, as that amounts to just overwriting the postpro directory every time a new solve completes, without saving off the substeps. If a user didn't care about storing the intermediate results (or cared about disk space), they might set that option. This could be changed to a boolean instead, of storing results or not, but I originally picked "SaveStep" for consistency with other saving off options.

"MaximumImbalance" could be argued as advanced, as adjusting that from 1 ought to be very unusual.

"WriteSerialMesh" is similar in spirit to the "SaveStep" option, if a user didn't care for the extra storage. It's pretty self explanatory though, so again can be made advanced. This option is mainly helpful for cases where a user thinks they might conceivably want to restart on a larger machine, or generate an adapted mesh for use in a parametric sweep.

"UseCoarsening" could also be argued as advanced, given the default to false, and that any reasonable choice of "UpdateFraction" makes it practically unneeded. This argument can be extended to deleting coarsening from the main SEMR loop entirely, but I thought it would be likely that users might want to explore it as an option so opted to keep it.

I see that you already updated this in the last commit, but two comments:

If it would serve better to rename the top level Refinement to something else, I could be on board with that. config["Model"]["Adaptation"] or "MeshAdaptation" could work too, and the existing "UniformLevels" and other parameters wouldn't seem totally out of place in that section. I also do not think it is necessarily confusing to users to have config["Refinement"]["Tol"] or ["MaxIts"] (as currently, but without the "Adapt" prefix). There shouldn't be any conflict with any other tolerances associated with refinement. But if you feel it is necessary I guess we can keep it.

The main benefit to moving items to the advanced section is not really avoiding doc strings (they should have probably all have doc strings but this makes them harder to change quickly as the code changes, and so for now they are undocumented). It is to be clear to users what items you intend for them to specify and which ones they really never should unless doing development work. So, items in the advanced section are not included in the demo config file above the docstrings.

sebastiangrimberg · 2023-10-17T19:47:29Z

palace/utils/timer.hpp

@@ -54,6 +56,8 @@ class Timer
      "Estimation",
      "  Construction",
      "  Solve",
+      "Adaptation",


The "Estimation" and "Adaptation" headers should be combined into one "AMR" section I think.

The first reason I opted for separate is pragmatic. If running without adaptation, it doesn't appear in the block timer and the estimate section ends up indenting weirdly:

Solve Preconditioner Coarse Solve Estimation Construction Solve

Adaptation only shows up in the timing if the max its is greater than one, which makes the indenting of the estimation consistent. The second is more semantic, Adaptation here means the Mark and Refine stages of SEMR. If estimate were placed as a sub timer, then that implies solve should be too, as all four steps are substeps of that loop. Solve being a sub timer feels entirely incorrect, so transitively estimate as a sub timer does too.

Marking and refining themselves ended up being so quick I didn't bother breaking their timers out, but Rebalancing ended up being noticeable, hence the sub timer.

I would always have the top level object show up:

... Operator Construction Solve Preconditioner Coarse Solve Mesh Refinement // <- Or just Refinement? Or Adaptation? Estimation Estimation Solve // <- Lump in construction, probably don't need separate categories? Adaptation Rebalancing Postprocessing ...

The idea being that we wouldn't really be doing any estimation cost if we weren't going to try to do something with it (adapt the mesh), so grouping them together as the "net cost of AMR" seems like the right way to do it.

palace/main.cpp

palace/utils/configfile.cpp

palace/drivers/basesolver.hpp

palace/utils/iodata.hpp

palace/utils/geodata.cpp

test/examples/testcase.jl

scripts/schema/config/model.json

palace/utils/iodata.hpp

palace/utils/dorfler.cpp

palace/utils/geodata.cpp

sebastiangrimberg · 2023-10-19T22:33:15Z

palace/utils/geodata.cpp

+      mesh::DimensionalizeMesh(*smesh, iodata.GetMeshScaleFactor());
+      smesh->Print(serial);  // Do not need to nondimensionalize the temporary
+    }
+    Mpi::Barrier(comm);


Do we need a barrier here?

The barrier here is to match with the BlockTimer bt(Timer::IO) before the branch on Mpi::Root(comm). I copied the pattern in ReadMesh where the timer is instantiated for all ranks rather than just root, and kept coordinated by the barrier. Strictly speaking the other ranks don't need to be waiting here, then the timer could be moved inside the if, and only time the root processor, and the IO timing would then be desynchronized across processors.

Whoops sorry about that that, my mistake. Not sure how I missed the Timer::IO addition already as well. You're correct that everywhere there is IO done from the root rank we should end with a barrier so it is timed correctly.

palace/utils/geodata.cpp

palace/drivers/basesolver.cpp

sebastiangrimberg · 2023-10-19T22:55:44Z

palace/drivers/basesolver.cpp

+                 final_elem_count - initial_elem_count, initial_elem_count,
+                 final_elem_count);
+    }
+    else if (use_coarsening)


Feel free to correct me if this is a dumb question or thought here, but it seems like coarsening should be more useful than only after the final solve when ntdof >= dof_limit. Every time we solve, we get a new estimate which could in principle expose elements for coarsening, right? Now, maybe we only want to try to coarsen after a few rounds of refinement, for which we could have a config file option, but I don't see the value in coarsening only after you've hit your limit of the maximum size.

sebastiangrimberg

Second (more thorough) pass done. Everything is looking really great, this is awesome @hughcars! Mostly all minor comments sprinkled throughout.

- Primary method added to BaseSolver that implements a SolveEstimateMarkRefine loop. - Marking is performed using a parallelized Dorfler marking scheme, aka bulk chasing, wherein a minimum number of elements are marked for refinement in order to refine a fraction of the total error. This marking scheme is particularly efficient in combination with higher order schemes wherein the error can be highly concentrated in a tiny fraction of the total number of elements. - The loop has capability for handling nonconformal meshes, but in this commit that is disabled pending additional work to address issues within MFEM. As such, the options within the Adaptation portion of the config file are present, but the options to specify them have been commented out. - This also introduces an ability to perform rebalancing of conformal meshes through a gather/scatter model wherein the mesh is serialized to the root node, before being repartitioned. - Transient solves are not supported, as the ability to coarsen (only possible with nonconformal adaptation presently) is not present.

…re coarsening

- Removing MFEM exceptions - Use long long int for true dof count - SolveEstimateMarkRefine returns void - Move adaptation configuration options to Refinement level - Make dimensionalization of the mesh in PostOperator consistent

…ng in varying marking with number of processors

- Rebase on main - Fix square error aggregation patch - Rename GetMeshScaleFactor -> GetLengthScaleFactor - Move saving off adapt subdirectory to start of Mark Refine Solve Estimate loop to avoid duplication - Add communicator to Dorfler marking arguments - Remove unneeded forward declare and add in used inclue - Refactor communication in Dorfler to send packed MPI buffers - Move Rebalance into geodata, and make RebalanceConformal a geodata.cpp local function - Make IterationPostDir a basesolver.cpp local function - Fix a missing timer barrier - Revert reordering in json schema

sebastiangrimberg · 2023-10-24T17:04:49Z

palace/drivers/basesolver.cpp

    BlockTimer bt_io(Timer::IO);
+    Mpi::Barrier(comm);


I think this isn't necessary? If anything it should be before the timer, but I think it makes sense for the timer to include all time waiting.

I had put it there as I was worried about the root starting to move the post processing directory while another processor is still writing into it. Given that, I viewed that blocking as part of the IO, hence placing it within the barrier.

sebastiangrimberg · 2023-10-24T17:05:16Z

palace/drivers/basesolver.cpp

  }
-  return dir;
+  return fs::path(dir);


Is there an advantage to fs::path instead of just string as used everywhere else?

For this function it could be a string, the one place it gets called, in the SavePostProcess lambda, it would then be converted to a path immediately. Using fs::path allows for using the directory_iterator and copy within the lambda, in addition to the / operator. So this function could return a string, but fs::path would get called on the result immediately, and given the function is only used in this one place, it seemed neater to have it return the finished product. Either is fine with me though.

- Fix bug where coarsening across processor boundaries might segfault - Fix bug where dorfler marking would segfault for empty estimate array (occurs for empty partitions or sometimes coarsening)

…format

…shold in anticipation of transient AMR

sebastiangrimberg · 2023-10-26T22:33:05Z

palace/utils/configfile.cpp

@@ -287,6 +287,23 @@ void RefinementData::SetUp(json &model)
  uniform_ref_levels = refinement->value("UniformLevels", uniform_ref_levels);
  MFEM_VERIFY(uniform_ref_levels >= 0,
              "Number of uniform mesh refinement levels must be non-negative!");
+  adapt_tolerance = refinement->value("AdaptTol", adapt_tolerance);


Some thoughts on config file parameter naming:

I would prefer "Tol" and "MaxIts" to directly mirror the usage elsewhere. If we need to clarify what the tolerance and maximum iterations refer to, then the header of "Refinement" should be updated. But I think any user will see this and there isn't any confusion associated with other refinement tolerances/maximum iterations.

Similarly, "DOFLimit" to "MaxSize" with appropriate description in the docs?

- Changed DOFLimit -> MaxSize, AdaptTol -> Tol, AdaptMaxIts -> MaxIts - Added MaxSize < 1 implying no size limit during an adaptation. - Caught a few missed deletions of coarsening.

…timate-mark-refine-dev

hughcars added enhancement New feature or request draft Work in progress labels Oct 16, 2023

hughcars force-pushed the hughcars/amr-solve-estimate-mark-refine-dev branch 4 times, most recently from a632258 to 2b66e53 Compare October 17, 2023 17:13

hughcars marked this pull request as ready for review October 17, 2023 17:39

hughcars requested a review from sebastiangrimberg October 17, 2023 17:39

hughcars removed the draft Work in progress label Oct 17, 2023

sebastiangrimberg suggested changes Oct 17, 2023

View reviewed changes

hughcars force-pushed the hughcars/amr-solve-estimate-mark-refine-dev branch from 021117b to 14bb828 Compare October 18, 2023 20:21