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

Unify assembler structure #1785

Merged
merged 5 commits into from Nov 22, 2017
Merged

Conversation

gassmoeller
Copy link
Member

Our current assembler system is somewhat messy, because the structure for how an assembler has to look like is not well defined. Essentially any function can be attached to the available signals, and we do not know which functions are already attached (because we can not check which type each connection has), or which quantities an individual assembler needs to work correctly. This has lead to the extremely confusing logic inside of set_assemblers, and to the fact that for every new assembler we need to modify three objects: Connect to the associated signal. Add the parent object to assembler_objects. Modify the correct assembler_properties object to compute the correct quantities before calling the assembler.
If instead we had an interface-plugin system as for the other parts of ASPECT, all this would be handled in one object, the individual assembler plugin. It would have a defined function that is called by Simulator (which replaces the signal), and could (optionally) implement another function that specifies which quantities need to be provided (which replaces the properties object). We could still keep the set_assemblers signal that allows every part of the code to modify the list of existing assembler plugins with the added benefit that we can then remove specific plugins, instead of wiping them all and starting from scratch, which is our only option right now.
This is a very early draft on which changes we would need to implement in order to move in such a direction, but I wanted to start a discussion about it first, before spending more time with it. The important bits can be seen in assembly.h below. For now I moved the StokesPreconditioner assembly out of the StokesAssembler class and made it a separate fully-functional assembler plugin. It still connects its execute function to the old signal, but we could as well call it via the pointer to the base class stored in assembler_objects. The one disadvantage of the new system is that we have less flexibility of providing plugin-specific information into the assembler (see my workaround for pressure_scaling below). However, I would argue that we can provide all necessary information in two ways: Either add it to the scratch object that is already designed to store all necessary information or provide it via the SimulatorAccess class.

I would like to ask for opinion from @tjhei and @bangerth (and all other interested parties). I think the new way would keep the code cleaner and will be easier to understand. Everything else (performance, flexibility) should be similar to the current approach with signals.

@bangerth
Copy link
Contributor

Yes, I like this approach. I think it is fair to say that my idea with the signals has failed and collapsed under its own weight.

@gassmoeller
Copy link
Member Author

Do not be so harsh on your idea 😏. It works very well for simple functions. Our assemblers are just getting too complicated. And we should definitely keep the set_assemblers signal to allow users to modify the selected assemblers.
I will continue working on this after next week.

@tjhei
Copy link
Member

tjhei commented Jun 1, 2017

Yes, this is what I have been thinking about for a while. Should we delay this re-shuffling until after we are done with the whole Newton stuff?

@gassmoeller
Copy link
Member Author

Yes lets wait until after the newton assemblers are merged, I discussed this with @MFraters in #1784.

@gassmoeller gassmoeller force-pushed the extend_assembler_base branch 5 times, most recently from 8084b7d to 935f229 Compare July 3, 2017 15:48
@gassmoeller
Copy link
Member Author

I am reasonably far with the rework of the assembler structure, but I would like to discuss a few interface questions before asking for a formal review, because they might change quite a few things:

  1. The new Assemblers::Interface base class has one execute function that has parameters of type ScratchBase *. This allows every assembler to implement that single function, but it also means every assembler has to dynamic_cast the type of scratch to its appropriate scratch object (i.e. StokesPreconditioner, StokesSystem, or AdvectionSystem). I could instead declare three separate overloads for the execute function in assemblers/interface, with default implementations that throw assertions. This would shorten each Assembler, and it would also implicitly check that each assembler is used only in the assembly of the correct system. I would prefer the three function way, but wanted to ask for opinions before changing all those lines again (it would significantly shorten this PR, and also require less changes to user written assemblers).

  2. The scratch and data parameters are conceptually similar to MaterialModelInputs and MaterialModelOutputs, however contrary to the latter, scratch is currently not handed over as const. This leads to the strange situation that the order of assembler calls can matter, if an assembler modifies the scratch object. Does anyone see a necessity for this behavior, or can I just make scratch const?

  3. Many assemblers need a cell, and face (and some face_no) input, but they are currently not explicitly stored in scratch. They can be reconstructed from the fe_(face_)values object, but it is quite a bit of replicated code. Unless somebody objects I would prefer to store those explicitly in the scratch objects to make assemblers simpler.

@gassmoeller gassmoeller changed the title [WIP] Unify assembler structure Unify assembler structure Jul 10, 2017
@gassmoeller gassmoeller force-pushed the extend_assembler_base branch 6 times, most recently from 8961d04 to 76a5e23 Compare October 23, 2017 19:49
@gassmoeller
Copy link
Member Author

I updated this PR to incorporate the Newton assemblers. This is ready for a review.
Sorry for the long changelog, but it really mostly is a conversion of single functions to assembler objects.

@bangerth
Copy link
Contributor

There's now a merge conflict. Can you resolve it?

@gassmoeller
Copy link
Member Author

The merge conflict is resolved and the tests pass. I am not sure why the astyle tester fails in a file that this PR does not even touch though.

@MFraters
Copy link
Member

I had the same problem in #1953, that a unchanged file failed the astyle test. A rebased 'fixed' the problem. I am not sure what causes this problem.

@gassmoeller
Copy link
Member Author

Thanks Menno, rebasing did the trick. Can we move forward with this PR, and block other changes at the assemblers for the moment? Every time something changes in the assemblers there will be a merge conflict with this PR, and I feel this would really improve the structure of our assembler system.

@bangerth
Copy link
Contributor

This is difficult to review for me because I haven't paid attention to how the assemblers are built these days. The patch is already pretty unreadable :-) (Not because it's poorly written, but because it moves so much code around.)

Do any of the other reviewers have a better overview of how this all works and could do the review?

@jdannberg jdannberg self-assigned this Nov 2, 2017
@jdannberg
Copy link
Contributor

I'll have a look tomorrow.

Copy link
Contributor

@jdannberg jdannberg left a comment

Choose a reason for hiding this comment

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

I'll look at the rest tomorrow.

internal::Assembly::Scratch::StokesSystem<dim> &scratch,
internal::Assembly::CopyData::StokesSystem<dim> &data) const
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
Copy link
Contributor

Choose a reason for hiding this comment

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

Can you fix the number of spaces so that the &s are aligned?

/**
* The number of the face object with respect to the current
* cell on which we operate. If we currently
* operate on a cell, this is not meaningful.
Copy link
Contributor

Choose a reason for hiding this comment

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

What happens if we operate on a cell, I mean, what's the value of face_number in that case? It could make sense to have that in the description here, so that it can be checked in the code.

HeatingModel::HeatingModelOutputs face_heating_model_outputs;
HeatingModel::HeatingModelOutputs neighbor_face_heating_model_outputs;

const typename Simulator<dim>::AdvectionField *advection_field;
Copy link
Contributor

Choose a reason for hiding this comment

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

Would it make sense to write a short comment for the AdvectionField (as for the other parameters)?

@@ -0,0 +1,491 @@
/*
Copyright (C) 2016 by the authors of the ASPECT code.
Copy link
Contributor

Choose a reason for hiding this comment

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

Copyright year?

* in Simulator::set_assemblers() (which is the only place where
* we know their actual type) and destroyed in the destructor of
* the Simulator object.
*/
Copy link
Contributor

Choose a reason for hiding this comment

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

That's the same documentation as for the Interface class...

* that correspond to velocity or pressure degrees (or, in fact
* any other variable outside the block we are currently considering)
*/
std::vector<types::global_dof_index> local_dof_indices;
Copy link
Contributor

Choose a reason for hiding this comment

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

Empty line missing.

* internal::Assemblers::AssemblerBase::create_additional_material_model_outputs()
* functions from each object in Simulator::assembler_objects.
*/
void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &) const;
Copy link
Contributor

Choose a reason for hiding this comment

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

You moved this function to the assemblers interface, correct? But I didn't see this extensive documentation over there. Did you think it wasn't necessary?

Copy link
Member Author

Choose a reason for hiding this comment

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

I moved it over. Copy-Paste error.

"rate is constant.")
}
}
#include "../cookbooks/inner_core_convection/inner_core_assembly.cc"
Copy link
Contributor

Choose a reason for hiding this comment

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

Perfect, thanks! I must have missed that when I went over the benchmarks to remove the duplicate code.





Copy link
Contributor

Choose a reason for hiding this comment

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

That's super nice that all this stuff doesn't have to be in here any more.

* equation for the current cell.
*/
template <int dim>
class AdvectionSystem : public Assemblers::Interface<dim>,
Copy link
Contributor

Choose a reason for hiding this comment

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

One point that confused me about this design, and I guess we can talk about that, is that now, if you call an assembler, you don't know any more if that assembler changes the Stokes or the Advection system. Before, with the signals, they were all in one class that either said advection or Stokes or other or whatever, but now they are just in different files. So if the assembler is just called AdvectionSystem, that's totally fine, but if it's called something weird like PhaseBoundaryAssembler, it's suddenly not clear anymore what it does.
So I guess we can either say we don't care, or we may want to think about some naming convention, or maybe you can think of another way for how to make this clearer (I think, in particular as the set_assemblers function is quite long, it could make sense to have names that are easy to understand).

Copy link
Contributor

Choose a reason for hiding this comment

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

As discussed, one could also use different namespaces for the Stokes/Advection/Other assemblers. Then, the assembler would only have to be called IncompressibleTerms instead of StokesIncompressibleTerms, and would be in the Stokes namespace.

Copy link
Member Author

Choose a reason for hiding this comment

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

I like the namespace idea (disclaimer: it was mine) 😄. If nobody objects, I will regroup the assemblers accordingly.

const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
const double pressure_scaling = scratch.pressure_scaling;

// First loop over all dofs and find those that are in the Stokes system
// save the component (pressure and dim velocities) each belongs to.
Copy link
Contributor

Choose a reason for hiding this comment

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

I discussed this with @gassmoeller, and this is not something for this pull request, but something that relates to #733 and could be addressed later. To avoid duplicating so much code in the different assemblers, we could move at least the lines below where things in the scratch object are filled out of the assembler, and into the local_assemble_... functions. Many of those things are similar for the all Stokes assemblers, or all Stokes preconditioners, and while we would have to fill the things needed by all assemblers in the local_assemble_... functions (so potentially also some things that wouldn't be needed later on), it might still not be slower in the case of multiple assemblers being used together, as we would have to go through the loops only once.
And it would have the advantage that the assemblers would be quite a bit shorter.

* This class assembles the right-hand-side term of the Stokes equation
* that is caused by the compressibility in the mass conservation equation.
* This class approximates this term as
* $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$
Copy link
Contributor

Choose a reason for hiding this comment

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

That is not the correct documentation.

* equation for the current cell.
*/
template <int dim>
class AdvectionSystem : public Assemblers::Interface<dim>,
Copy link
Contributor

Choose a reason for hiding this comment

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

As discussed, one could also use different namespaces for the Stokes/Advection/Other assemblers. Then, the assembler would only have to be called IncompressibleTerms instead of StokesIncompressibleTerms, and would be in the Stokes namespace.

Copy link
Contributor

@jdannberg jdannberg left a comment

Choose a reason for hiding this comment

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

I've seen a lot of this pull request now, and I think I understand the general concept. I like the concept of having classes for the assemblers rather than signals, and I think this simplifies the code a lot, and also makes it easier to debug things. Thanks for doing this!

@tjhei: As we've discussed the structure of the assembly extensively with you, maybe you could have a look at the basic structure of these changes (basically the new assemblers/interface.h) and say if you are okay with the new interface?

std::vector<double> phi_p_c;
std::vector<Tensor<1,dim> > grad_phi_p;

double pressure_scaling;
Copy link
Contributor

Choose a reason for hiding this comment

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

Question: Your assemblers are classes now, and you derive them from SimulatorAccess. So you could also just get the pressure scaling from SimulatorAccess, instead of making it a member variable of the class. Then you wouldn't have to worry about initializing it, or if it's set in local_assemble_stokes.

* This class assembles the right-hand-side term of the Stokes equation
* that is caused by the compressibility in the mass conservation equation.
* This class approximates this term as
* $- \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$
Copy link
Contributor

Choose a reason for hiding this comment

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

You should mention that it uses the depth derivative of the reference density (in \frac{\partial rho}{\partial z}).

Copy link
Contributor

Choose a reason for hiding this comment

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

Oh, and maybe that it's the explicit version.

* It includes this term implicitly in the matrix,
* which is therefore not longer symmetric.
* This class approximates this term as
* $ - \nabla \mathbf{u} - \frac{1}{\rho} * \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$
Copy link
Contributor

Choose a reason for hiding this comment

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

Same thing about the reference density.

* This class assembles the right-hand-side term of the Stokes equation
* that is caused by the compressibility in the mass conservation equation.
* This class approximates this term as
* $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$
Copy link
Contributor

Choose a reason for hiding this comment

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

Maybe mention that the compressibility is used in place of \frac{1}{\rho} * \frac{\partial rho}{\partial p}.

* This class assembles the right-hand-side term of the Stokes equation
* that is caused by the compressibility in the mass conservation equation.
* This class approximates this term as
* $ - \nabla \mathbf{u} = \frac{1}{\rho} * \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$
Copy link
Contributor

Choose a reason for hiding this comment

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

This also doesn't look like the right documentation.

* get_additional_output() returns an instance before adding a new
* one to the additional_outputs vector.
*/
virtual void create_additional_material_model_outputs(MaterialModel::MaterialModelOutputs<dim> &) const {}
Copy link
Contributor

Choose a reason for hiding this comment

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

I don't really understand how this works now. Before, this was a function that was part of the simulator, and now you moved it to the assembler interface, but only some assemblers implement it? For example, for Stokes only the StokesIncompressibleTerms. What if I don't use that assembler?
And you also changed some of the test results that tested if this is working correctly, so it would be great if you could explain this a bit more.

Copy link
Member Author

Choose a reason for hiding this comment

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

If you do not use this assembler then this additional material model outputs object does not get created. But that is fine, each Assembler is only responsible for creating the additional material model outputs that it actually uses during assembly. This also means each assembler has to check if the relevant additional material model outputs already exists (because another assembler already created them).

namespace OtherTerms
template <int dim, class AssemblerType>
void
initialize_simulator(Simulator<dim> &simulator,
Copy link
Contributor

Choose a reason for hiding this comment

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

Can you document that function?

initialize_simulator(*this,assemblers->advection_system);
initialize_simulator(*this,assemblers->advection_system_on_boundary_face);
initialize_simulator(*this,assemblers->advection_system_on_interior_face);
initialize_simulator(*this,assemblers->advection_system_residual);
Copy link
Contributor

Choose a reason for hiding this comment

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

Hm, this is not really nice. Isn't there a way to do this for all of them at once?

Copy link
Member Author

Choose a reason for hiding this comment

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

It turns out there is a way, but for that I need to implement a constructor for all assemblers (this constructor would take the this pointer to the simulator and call the constructor of the SimulatorAccess class that immediately does what initialize_simulator does). But this needs more lines than I have here. Still if you would prefer the constructor I can implement that (it might look cleaner, but require more lines).

Copy link
Contributor

Choose a reason for hiding this comment

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

Hm. What kind of error message do you get if you forget to call this line for your assembler? I am just thinking, if you make a new type of assembler and forget to implement this line, how easy is it to find out what the problem is? If that is obvious, I would be okay with leaving it like it is.

@@ -1001,6 +399,9 @@ namespace aspect

// Prepare the data structures for assembly
scratch.finite_element_values.reinit (cell);
scratch.cell = cell;
scratch.pressure_scaling = material_model->reference_viscosity()
/ geometry_model->length_scale();
Copy link
Contributor

Choose a reason for hiding this comment

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

As I mentioned earlier, it would be much better to use the simulator access for this instead.

@@ -1189,6 +592,9 @@ namespace aspect

// Prepare the data structures for assembly
scratch.finite_element_values.reinit (cell);
scratch.cell = cell;
scratch.pressure_scaling = material_model->reference_viscosity()
/ geometry_model->length_scale();
Copy link
Contributor

Choose a reason for hiding this comment

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

Same point about the pressure scaling.

@gassmoeller
Copy link
Member Author

I think I addressed all comments, except the discussion points I answered above.

Copy link
Member

@tjhei tjhei left a comment

Choose a reason for hiding this comment

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

This looks very good! I left a bunch of minor comments/suggestions. In general we should figure out if all implementations should go outside the class declaration or into the .cc files instead.

s/struct AssemblerLists/class Manager/g
s/AssemblerLists/Manager/g
s/AssemblerBase/Interface/g
s:assembly.h:simulator/assemblers/interface.h:g
Copy link
Member

Choose a reason for hiding this comment

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

we could keep assembly.h (that includes the other files) to avoid this change?!

Copy link
Member Author

Choose a reason for hiding this comment

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

I would vote for removing the file. The interface is broken anyway for anybody who wrote an individual assembler (not many people), and those would be the only ones including assembly.h. interface.h is the new assembly.h.

class ApplyStabilization: public Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
void execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch,
Copy link
Member

Choose a reason for hiding this comment

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

why does this work without "public:" ?

public SimulatorAccess<dim>
{
public:
virtual ~MeltInterface () {};
Copy link
Member

Choose a reason for hiding this comment

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

I don't think you need the destructor here as this is still an abstract base class

public SimulatorAccess<dim>
{
public:
virtual ~NewtonStokesReferenceDensityCompressibilityTerm () {};
Copy link
Member

Choose a reason for hiding this comment

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

remove semicolon. Do we normally have all those in the .cc files? I don't mind having them here...

{
namespace Assembly
{
namespace Scratch
Copy link
Member

Choose a reason for hiding this comment

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

My vote would be to remove this namespace as it doesn't convey additional information (Scratch::ScratchBase vs ScratchBase)

// lines in set_assemblers(), where this operation appears multiple times.
template <int dim, class AssemblerType>
void
initialize_simulator(Simulator<dim> &simulator,
Copy link
Member

Choose a reason for hiding this comment

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

can you make the simulator const here?

*/
namespace OtherTerms
// This function initializes the simulator access for all assemblers
// inside of the assemblers parameter. This is just a shortcut to save some
Copy link
Member

Choose a reason for hiding this comment

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

very cool!

for (unsigned int i=1; i<assemblers->advection_system.size(); ++i)
{
std::vector<double> new_residual = assemblers->advection_system[i]->compute_residual(scratch);
for (unsigned int i=0; i<residual.size(); ++i)
Copy link
Member

Choose a reason for hiding this comment

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

can you please use a different parameter than reusing i?

@@ -8,24 +8,14 @@ Number of active cells: 4 (on 2 levels)
Number of degrees of freedom: 84 (50+9+25)

*** Timestep 0: t=0 seconds
* create_additional_material_model_outputs() called
Copy link
Member

Choose a reason for hiding this comment

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

why is this no longer called?

Copy link
Member Author

Choose a reason for hiding this comment

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

Because previously if an assembler was in the list of assembler_objects create_additional_material_model_outputs() would be called for all equations, while now it is only called for the equations in which this assembler is actually used (in this case the Stokes equations). Since the Assembler can only be used in the Stokes equation (to successfully cast the scratch object), I can not reproduce the old output, but the new one is more reasonable anyway.

@@ -8,32 +8,14 @@ Number of active cells: 4 (on 2 levels)
Number of degrees of freedom: 84 (50+9+25)

*** Timestep 0: t=0 seconds
* create_additional_material_model_outputs() called
Copy link
Member

Choose a reason for hiding this comment

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

same here?

@gassmoeller
Copy link
Member Author

I think I addressed most comments. About the destructor implementations: I agree that we usually move those into the .cc file. It feels like a serious case of wasted lines of code in this case though. Each of the ~20 assemblers only implements one single function that is const, and none has any member variables, so having a destructor is kind of useless anyway. I would vote for keeping it in the .h files.

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.

None yet

5 participants