llnl / exaconstit Goto Github PK

Add the necessary additions to ExaConstit to allow it to build with HIP

While fixing issues related to the partial assembly on parallel, I noticed that the current formulation of the l1-jacobi does not work as expected. I need to rework the AssembleDiagonal formulation to return the correct l1-jacobi. In the meantime, a preconditioner is not being used with the partial assembly formulation.

As part of the merged in mfem v4.3 changes, we now have a dependency on RAJA >= v0.13. Although, I will note that this only appears to be a requirement as far as far as compiling on the GPU. Our tests still pass with RAJA v0.11.0 on the CI, so the changes appear localized to just GPU builds.

In the mfem v4.3 changes several mesh ctor methods were deprecated. We should should move to the new methods.

The build system needs to be updated so that the user doesn't have to worry about the various dependencies within MFEM. They should just be able to point to the mfem, ExaCMech, and RAJA builds and be good to go. The rest should be handled for them. It'll also make things easier for us as developers.

We should still require the user to specify if we're using CUDA or OpenMP, and if CUDA what is the machines architecture. We could also think about just going ahead and automatically setting the MPI and Fortran cmake flags internally to true as well, so they don't need to do that.

Helping a new user get everything installed and it was brought up that having an example GPU build script would be useful.

Also, it appears that our fork of MFEM that we use needs to be updated again. I ran into some cmake issues while helping this user out.

The scripts are a bit outdated with our move to mfem v4.3 which requires RAJA v0.13.0.

Also, I've found some issues using the host only version on a WSL2 machine in regards to what appears to be the ENABLE_SHARED_LIBS cmake define passed into ECMech... It's apparently catching itself in the BLT update. I'll open an issue in the ECMech repo about this.
edit: This flag issue appears to be a fault on my end. ECMech is using the correct cmake define in the repo. I was just using the deprecated version which must have worked with my older versions of cmake...

So, we'll want to update things to make sure things work there as well as we might have a few users using the script there.

Currently we don't have any guards checking to see if a property file actually exists and then fail if it does not exist on the system.

We should add those in there probably at the same time when we check the validity of all the options.

For some reason, it appears that if Caliper is found before at least RAJA then cmake crashes and complains about something related to camp. It's not clear to me why this is occurring at all, but it's clear that something weird is going on with cmake. However, if we find caliper last then everything works fine...

I ran into issues where if I had rcarson/mfem/exaconstit-dev fork being used with conjunction with a ExaConstit/mfem_forall_dev branch that the build was barfing when trying to compile the 1st file. If I didn't have ExaConstit configured to being built in debug mode it worked fine.

I need to run this issue down and see if it's a problem with my linux desktop build or if it's happening on other platforms as well. I've never had this happen before in the past when making use of debug builds with both systems at least on a LC machines...

While working with a number of workflow tools, it is starting to seem like it might be useful to make ExaConstit into a callable library. I don't imagine this requiring too many changes. The big thing would be that we would want to take in the same arguments we do now as well as an MPI_Comm argument.

Long term this action might also allow us to drive things from Python and feed data back there as well. However, we would need to break things down quite a bit more to get to that state and do some major refactoring.

Just to put some context into things, we're currently exploring https://emews.github.io/ as a ways to parameterize models on machines like Summit and Frontier.

It would be useful to make use of BLT's google test inclusions to test and make sure we don't see any regressions when things are added.

The PA_test is pretty much there it just needs to be swapped over to make use of google test. Next, a single xtal tension example that's 5x5x5 should allow us to test our voce and mts models and make sure no major changes are introduced. We should test our full and PA formulations here. During all cases, the debug device option should be used to make sure our code memory wise works on both GPUs and CPUs.

The new optimization framework can potentially generate a ton of files due to all the symbolic links to common files for each run conducted.

As part of the post-processing step on a successful run, we should delete all of the symbolic link files and leave only the bare minimum files such as the as the toml file, properties file, and simulation log and error fileS. The average stress and auto dt file are saved off in the checkpoint files. The optimized parameter sets are also saved in the checkpoint file, so if we really need to we could eliminate the property file as well. So, we can recover that data later on if need be. Overall, this would bring us down from something like 10-20 files to only needing 4 files. A drastic reduction in file count if you have a ton of genes being used in the optimization routine.

ExaCMech has a few new changes that break its current interface within ExaConstit. I'll need to update the interface for those. Also, we can get rid of hard coding the number of state variables are needed and actually just have them determined just from the ExaCMech classes themselves.

Also, the history variables have now been updated, so they include a flow strength term as well.

Hi,

I am trying to build ExaConstit and ran into an issue with the dependency on ExaCmech. ExaConstit suggests using tag 0.2.1 for Snls (dependency of ExaCmech) while ExaCmech is asking to build the "templated" branch of Snls. I am running into issues with both options, but before digging deeper wanted to know what version of Snls I should be using for ultimately building ExaConstit.

Also, I was building on Windows (Visual Studio). Is that a nonstarter for ExaConstit?

Thank you.

G Kumar

In order to simplify the addition of this, I'm working on adding a material-specific QuadratureFunction and QuadratureSpace into MFEM. After this is in there, we can start refactoring various parts of the code to work with these classes and finally get multi-material support with the initial caveat that they're the same class of models aka ECMech only or UMAT only and no mixing of the two. It's possible this might complicate things as well once ECMech gets native J2 and Barlat yld2004-18p models added to it. Once those are added, we'll probably also want to revisit and decide if a further restriction on only allowing crystal models to operate with crystal models is wanted or if we'll allow mixing of the two. I do see a benefit of mixing the two models, if you imagine having the J2 model surround the crystal material or have it on like the boundaries... Although, we'll get into some very fun issues related to load balancing with that mixing of things...

One thing that will need to be modified outside the physics classes is the options class/parsing as things will be more complicated. We'll need to allow multiply supplied crystal types, models, and property/state files. We'll have to add in a phase map file for the serial mesh. This phase map will get mapped to the mesh attributes after the initial mesh attributes/grain IDs are used to instantiate everything. Afterwards, we'll map the phase IDs onto the mesh attributes. If only MFEM supported multiple attributes for mesh element attributes this would be easier... One good thing is we can save the grain IDs later on to the data collections so we'll have them there.

The current post-processing stuff is a bit brittle and could be handled better. New models that don't fit the mold that's already within the code could easily break things.

The post-processing stuff really should be it's own class which holds all of the gridfunctions that we might want to use, and their respective FiniteElementSpaces kinda of like what's in the Axom::MfemSidreDataCollection for when they do reconstruction of stuff.

Although, it might also be worth thinking about whether or not we want this class to also be general enough to handle the restart functionality as well. I'm leaning towards keeping them separate just because that functionality could get more complicated if we decide to support restart capabilities through something other than Axom's sidre.

When I run make -j 4 in the ExaConstit build folder after running cmake, I get several errors which seem to be calls to functions with no corresponding definition in the source code. For example :

/path to exaconstit/ExaConstit/scripts/meshing/mesh_generator.cpp: In function ‘void setBdrConditions(mfem::Mesh*)’: /path to exaconstit/ExaConstit/scripts/meshing/mesh_generator.cpp:188:19: error: ‘class mfem::Mesh’ has no member named ‘SetBdrAttribute’; did you mean ‘GetBdrAttribute’? mesh->SetBdrAttribute(i, 1); // bottom

I did a reecursive grep for the function SetBdrAttribute and it does indeed seem that it is not defined. I used the version of mfem specified in the ExaConstit installation notes : the exaconstit-dev branch of MFEM.

While the post-processing tools inlined into ExaConstit are worked on, we should be able to offer python tools that perform the same set of operations such as those outlined in https://doi.org/10.1088/1361-651X/aa6dc5 . We could also look into showing off other tools as well to analyze the data.

The current code doesn't allow any of the boundary conditions to change over time. In order to run more interesting simulations, we need to be able to provide users the ability to modify the boundary conditions over various time steps.

Most of the code to allow for this should already be there. However, we'll need to do a few things first:

• Create a good file format that makes it easy for users to modify the boundary conditions at a given time step
• Additional toml options are needed to signal that modified BCS are being used
• Methods written to modify BCS need to work with partial assembly structure

Overall, the design for this should be for our users to write out the evolving BCS. Next, we also want the new code design to be developer friendly. If we add more advance BCS (for example traction BCS) later on for our users the design should be adaptable to allow for this.

I noticed while making a test to check that my partial assembly implementation is correct that the UMAT object deconstructor errors out when a QuadratureFunction is deconstructed. I need to look into what is causing these issues but it appears to be something within one of these variables:

 mfem::QuadratureFunction loc0_sf_grad;
 // The incremental deformation gradients.
 mfem::QuadratureFunction incr_def_grad;
 // The end step deformation gradients.
 mfem::QuadratureFunction end_def_grad;

It's most likely an issue in one of the following functions:

init_loc_sf_grads(fes);
init_incr_end_def_grad();

since I found that the object properly deconstructed when I'd commented these lines out.

Once I get the partial assembly and GPU implementation working, I'll get back to figuring out what the underlying issue is.

The integrator is added twice to the internal NonlinearForm which would cause performance issues and slight numerical issues.

Essentially, the system being solved in this case is:

(K_mat + K_mat) x = (f_int + f_int) rather than K_mat x = f_int.

I'll push a fix for this shortly onto the dev branch.

Now that the element assembly is supported within MFEM, it opens up another method we can use to get on the GPU. I don't necessarily this should be too difficult, but I will probably need to rework some of the basics from the AssembleElementVector(Grad) functions to be more inline with the type of formulations seen in the mfem::EABilinearFormExtension classes.

Hopefully, we should find that this is more performant than the current PA method that we have in the code. Since we do have some info about the matrix, it should be possible to also use more preconditioners with this method.

We really should allow this to be changed from the toml. It'll allow us to have multiple running tests at a time within python. Also, it allow for batch jobs to have unique file names for their stress values. I'll add this support within the #26 PR. Since, it already contains changes related to the necessary folders.

Been meaning to do this for awhile

These classes now exist within the master branch of MFEM as well as our own branch of MFEM. So, we now longer need them within ExaConstit. Currently, I've got a temporary fix to deal with the name class between the MFEM QuadratureFunctionCoefficient class and our own by just renaming ours.

The removal will require some changes to our code, since the MFEM classes only hold constant references to the QuadratureFunction. My plan is to no longer use them within the model class and instead move them to being only used within the SystemDriver class during the post-processing step. The ExaModel class will now just hold pointers to our quadrature functions. It honestly makes more sense this way, since we never made use of the Coefficient classes within the ExaModel class.

I've been running some MTS DD type models for an FCC material. I'm finding that our global solver either using the plan newton raphson or newton raphson with a line search method is having trouble converging. It seems as if we could really benefit from that addition for these tougher problems.

Once, I've ported that solver we should probably donate it to the MFEM project as their community would also benefit from it as well.

The current kernels are all based on the most general non-tensor element formulation. When tensor type elements are used (hexahedron), we can reduce some overhead in our kernels and potentially speed-up our overall assembly kernel.

I had a user reach out about needing constant strain rate for some of their exploratory studies. It seems like I could probably follow something like what was discussed in the comment thread at #26 (comment) .

It'll probably take a little bit of work to get a POC going, but I don't think it should be too awful with the changing BCs functionality now within the code.

As a first pass, the easiest version will be something like a forward eulerian version of things where we calculate the velocity conditions based on the previous time step nodal locations. After we get that looking, we can look at more expensive versions like ones that do an iterative approach to calculate at the end time step or mid-time step values.

Additionally, I plan on allowing this took work with meshes that don't have a min at origin. Also, I plan on requiring the user to define what degree of freedoms are free and which ones are constrained to the velocity gradient condition.

As part of this set of calculations, I figured I would also calculate the Eulerian strain for the user as we can easily get out the deformation gradient, and so that value isn't hard to get at all. I figure this would make it much easier for them when calculating post-processing values as they'll automatically have the "true" strain values. If they want the lagrangian strain tensor, they can just get that from the deformation gradient...

A typical light-up script would require an elements elastic strain tensor, lattice orientation, and element centroid to work. The element centroid and strain tensor are not currently outputted, but they should be easy enough to add. Although, the elastic strains from ECMech will require just a little bit of extra work as we'll need to add the volumetric contribution to the stored 5-vec form of the deviatoric elastic strain stored in the history variables. We also save off the relative volumetric expansion which will allow us to get at the volumetric strain as well.

For some reason when compiling things on the GPU with gcc/7.3.1 and cuda/10.1.243, the following set of lines:

As far as I can tell this is due to the templates on these functions which is perfectly legal c++ code and that the compiler can't reason about how to capture them. Therefore, I'm going have to revert the templated nature of the functions and just have them passed in as const bools in the function. This fix will require a bug in the point release number to 0.4.1 and the same to a tag release. I still have no idea why this is only captured for these MFEM_FORALL loops and not any earlier ones.

I was testing out a different set of boundary conditions from my usual symmetric BCs where the loading is in the z direction.

The below case failed:

[BCs]
    #essential BC ids for the whole boundary
    essential_ids = [1, 2, 3, 6]
    #component combo (x,y,z = -1, x = 1, y = 2, z = 3, xy = 4, yz = 5, xz = 6, free = 0)
    essential_comps = [3, 1, 2, 0, 0, 2]
    #Vector of vals to be applied for each attribute
    #The length of this should be #ids * dim of problem
    essential_vals = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,  0.0, 0.001, 0.0]

However this case succeeded:

[BCs]
    #essential BC ids for the whole boundary
    essential_ids = [1, 2, 3, 4, 5, 6]
    #component combo (x,y,z = -1, x = 1, y = 2, z = 3, xy = 4, yz = 5, xz = 6, free = 0)
    essential_comps = [3, 1, 2, 0, 0, 2]
    #Vector of vals to be applied for each attribute
    #The length of this should be #ids * dim of problem
    essential_vals = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.001, 0.0]

The error thrown out was this:
terminate called after throwing an instance of 'std::out_of_range'
what(): _Map_base::at

Clearly the map was unhappy with having non-sequential indices provided. I'll need to look around and see how the std::map is being used.

Dear developpers of ExaConstit,

I am interested in your project for my own developments of a mechanical solver based on MFEM.
@rcarson3 already answered several of my questions on MFEM and told me of the ExaConstit project as a source of inspiration. I must thank him for this suggestion as it illustrates nicely how to build the elementary inner forces and stiffness matrices.

However, I am always a bit lost when the behaviour integration is called.

If I am not mistaken, in ExaConstit, it is called in the SetUp method of the NonlinearMechOperator class which inherits from the NonlinearFrom class. This SetUp is called in the Mult method which is part of the NonlinearForm interface. The comments in the Mult function seems to state that this call is more or less unatural...

I have a few questions about this.

I would expect to integrate my behaviour in a class of the type NonlinearFormIntegrator, as this is done in the hyperlelastic examples in MFEM. In this case, I need to implement the AssembleElementVector (inner forces) and AssembleElementGrad (stiffness matrix) (Let alone partial assembly for the moment). But in this case, I don't know where to integrate the behaviour (unless I do some assumptions like AssembleElementVector is always call before AssembleElementGrad). It this the reason why NonlinearMechOperator inherits from NonlinearFrom ? If not, can I only implement a NonlinearFormIntegrator ?

Thanks for any help,

Regards,

Thomas

For large simulations, the ability to restart the simulations is vital.

We need to be able to save out the mesh, velocity field, current time step, end-time step state variables, and end-time step stress state. All the other data can be recovered from our initial data files. In order to do this we'll need to help expand some of the MFEM data collections for this purpose. Currently, the various data collections do not support saving out or loading QuadratureFunction data which is needed for the state variable and stress states.

Conduit and ADIOS2 data collections might be our best bet to get these fields implemented.

We should be able to piggy off a bit of the work being done in mfem/mfem#1991 to get things set-up within here. I don't think we'll need all of the TPLs for our github actions. However, we'll definitely need to make use of RAJA, ExaCMech, Metis v5, Hypre v2.18.2, and MFEM. Additionally, we should be able to just store the built libraries from PR to PR. So, it just has to build ExaConstit and run our test suite.

During my post processing updates, I noticed that the Jacobian calculations is using the beginning of time step coordinates rather than the end of time step coordinates of the mesh. I saw this in the deformation gradient calculations when it was giving me the previous time step values. I'm guessing this error was introduced at some point when doing all of the refactoring for the GPU stuff. Since, it was originally working correctly with the UMAT stuff based on the deformation gradients that were returned back then. Nominally, I wouldn't expect this to make a drastic difference if you have small time steps or going out to small strains. However, it'll definitely cause issues as your strains increase or you're taking large time steps.

I'll push a fix shortly and then will do a new release once I've quantified how a big a difference this makes for large strains. I'll probably do 2 commits for this. One with the fix and then one with the test rebases, so it's easier to see how things have changed.

Currently, data that is saved off to be post-processed is saved off as L2-type elements that have an order of 0 (essentially elemental data). I've recently been experimenting with saving the data off as (p-1) order where p is the order of your field / H1 type elements. I've noticed that this can lead to some interesting oscillations with whatever variable your plotting outside of the solutions variables. An example of this is the below figure that has the von Mises stress plotted for a single crystal at 1% strain. The subfigure on the left makes use of linear hex elements and the right makes use of quadratic hex elements, so the von Mises fields would be constant and linear respectively in each subfigure. It appears to me that the elements need to take in information either about their local neighbors or they need to have their field interpolated across a given grain with L2 projection type method. Depending on the type of data analyst you plan on doing the answer to this question might change. I haven't really thought it through too much other than it will introduce some order of error into the data analyst, since it's all interpolated from the data that lies at the quadrature points.

It's not critical to fix now, but it is something that we'll want to address in the future. I'll need to sit down with some of the MFEM folks and see if we can't come up with a nice solution to this problem.

One last point I think should also be thought about is that orientation type data might be a bit trickier to deal with depending on what orientation space we're in.

In order to use a standard version of MFEM, we need to port a number of features over and swap some of our original implementations for PA/EA assembly methods over to the equivalent MFEM versions.

• Port PA nonlinear extension type class to MFEM's version
• Add EA nonlinear extension type class to MFEM
• Work on splitting up mfem/mfem#722 into smaller PRs
• Create alternative versions of NonlinearForm::SetEssentialBCPartial and GridFunction::ProjectBdrCoefficient that don't depend on intrinsics within ExaConstit aka the boundary condition manager
• Maybe some new code to deal with changing BCs for our solves as MFEM proper has gotten rid of our old workaround...

ExaCMech now has support for HCP materials when using the Kocks Mecking DD MTS models. We should add support for those crystal types within ExaConstit at least for the MTS models.