Diffusion [non-local plasticity is still in progress]

[alarshi]

This branch has implementation of the diffusion scheme written by Julianne in a simple material model. A test is included for simple model but more changes for diffused strain weakening is yet to be done.

[jdannberg]

@alarshi, you mentioned you wanted to push some commits to this pull request before I should look at this? Or is this the state where you had the problem with the compositional field that you want to advect being solved correctly only once?

[alarshi]

@jdannberg Thank you for looking into it! I pushed the material model i modified from visco plastic: diffused plastic. It is exactly the same but has only copy outputs and not plastic additional outputs as the named additional outputs. The parameter file has a field called 'seed' which is not solved after the first time step and I am not sure how to fix it. Thanks again! :)

[jdannberg]

I think the problem is that this function goes through every compositional field and interpolates the copy outputs onto all of these fields. I think what you really want to do here is instead of looping over all the fields (below), you want to pass this function the compositional field that you want to interpolate your outputs onto as an argument. Then, in solver_schemes.cc, you could call interpolate_material_output_into_field(c), and then is this function the interpolation would only be done for that specific field.

[jdannberg]

This is the place I was talking about above. All the fields are copied, and not just the ones you want. So here and in line 1804, you want to remove the loop (and c will just be a function argument).

[jdannberg]

Edit: I guess what I really wanted to say was:
So you want to change that here too, using interpolate_material_output_into_field(const unsigned int c).

[jdannberg]

Alright, I made some comments in the code based on what I think might be the problem, but the test file visco_plastic_copy_and_diffuse.prm still crashes for me (although it solves for the seed field now). So let me know if that helped, and if you encounter any other problems.

[alarshi]

Thank you @jdannberg. I will have to work more on the test but this is a huge help!

[alarshi]

@jdannberg : After running some tests from simple setups, it appears to me that without diffusion and just copy and diffuse field, the results are not consistent with corresponding test not using this compositional field method. If we run the parameter file copy-and-diffuse.prm without the diffusion or the diffusivity parameters enabled, i see that the density output is copied into the composition but over the time the density outputs do not match the density outputs if we were not using the copy-and-diffuse method. I do not understand the reason, could you help me why this is happening?

[jdannberg]

@alarshi, I'll have a look at this tomorrow.

[jdannberg]

@alarshi, can you explain a bit more which parameter files exactly you ran that gave you different results? One reason I can imagine why the results would still be different, even if the diffusivity is set to zero, is that a field that evolves using the copy and diffuse method will not be advected, and in the composition-copy-and-diffuse.prm, the compositional field has an influence on the density. So over time, the density will evolve differently, because if you use any other compositional field methods, the field will be advected (and not copied), so the model will have a different density. If you want to see the same density in a model with fields that are advected normally, and fields that use the copy and diffuse method, then your compositional field that you use to track the density has to be passive, without an influence on the density.

Or is that not the problem you are seeing?

[alarshi]

@jdannberg Thank you, this is the exact problem I was trying to understand. I used composition-copy-and-diffuse.prm parameter file, I still need to make sense of the values of the density in the results but it is much clearer now.

[alarshi]

@jdannberg : I am still confused in the results. If we run the current parameter file with copy-and-diffuse field and higher differential density (which is copied as the composition field in our case), say 0.4, I see an increase in density from 1.4 at time_step=0 to 1.56 at time_step=1 for the denser circular composition in the center. With much lesser diffusivity (10^-5 times), I do not see the material diffusion or the hot center decreasing in size (as expected) but the values of the density still increase similar to the case with lower diffusivity. In other words, changing diffusivity did not change the rate at which the background material density is increased. I do not quite understand these results in terms of reaction of the composition with the density.
Could you help in this, whenever you have time?

[jdannberg]

I think the problem is that the model setup is a little confusing, as the model density depends on the compositional field that you use to track the density, so you would not set up a real model like this. Essentially what happens is that the density is copied into a compositional field, then it is diffused, and then the actual density in the model is computed as the reference density + the Density differential * c.
So even if the compositional field c was constant everywhere (and then it would not matter if there was diffusion or not), your density would still change every time step, as

new density = reference density + Density differential * old_density

So in our example, the background density starts out as being 1 (before time step 0), which is the reference density. Assuming a Density differential = 0.4, that gives us a new density = 1 + 0.4 * 1 = 1.4 after time step 0, and a new density = 1 + 0.4 * 1.4 = 1.56 after time step 1 and so on. So the background density is totally independent of the diffusivity.

In a real model, you wouldn't set up the compositional field in this way if you wanted to use it to track the density. You would start out initializing the compositional field with the real density, and then having a reference density of 0 and a Density differential of 1 for the field that tracks the density (so essentially, you would use your new compositional field as the density). Just to demonstrate how that would look like, I've modified the composition-copy-and-diffuse test like that (see attached file), if you run that, it should be much easier to understand.

composition-copy-and-diffuse.prm.txt

[alarshi]

Ah, I understand now, thank you so much!

[jdannberg]

@alarshi
I saw you pushed some commits, just let me know if you get to a point where you want me to have another look at this, or if you want help with rebasing this to ASPECT's current master (you might get lots of conflicts due to my pull request #2677, so I am happy to help with that). I feel it would be really useful to have this functionality in ASPECT soon (and I want to use something similar to your approach in my models too), so just ask if you need help with anything.

[alarshi]

Thank you so much @jdannberg for offering your help. I think we are close but not there yet.
Could you tell me if the outputs copied into copy and diffuse field (in the older version, it is called different now) are used in the next time step? For example, I use copy_out->strain_rate, would the diffused strain_rate be used in the subsequent time steps to calculate other variables like viscosity? I do not see the variation currently, maybe I am doing something incorrect.
If not, could you tell me how I can approach it so that the other variables calculated are dependent on the diffused field?

[jdannberg]

Yes, you can use the outputs copied into the copy and diffuse field in the next time step. They are only changed when you actually do the copying and diffusing (which is done, I think, after the temperature and the other compositional fields are solved). If you want to use this in your material model, you simply have to use in.composition[i][this->introspection().compositional_index_for_name("name_of_your_strain_rate_field")] in place of in.strain_rate[i]. (You can not simply use copy_out->strain_rate, because that would just copy the current strain rate, it is not the field that you have diffused, it's just the material model output that you used to copy the strain rate into a compositional field).

[alarshi]

Understood, thanks! :)

[alarshi]

@jdannberg : We realized recently that the equation we were trying to solve for the unknown, D, is - $\dot{D} - \kappa \nabla \dot{D} = source$ and not $\dot{D} - \kappa \nabla D = source$. So, we cannot use the diffusion solver and will have to think of another way to approach it.
I do not want to delay your work so please let me know and I will clean up the files and push it here. I did some tests on using the diffused plastic strain using the viscoplastic model to observe the effect of increasing diffusion with a mesh resolution and modifying resolution with constant diffusivity. Let me know if it will be of use to you.

[jdannberg]

So do you still have D anywhere in the equation (like in the source term), or only \dot{D}?

If you only have \dot{D}, I think it would be quite easy to still use the diffusion solver, you just have to change some terms (and even if you have D itself in the source term, I think this approach might still work). What we are doing at the moment is something like

D_new/Delta_t - \nabla \kappa \nabla D_new = D_old/Delta_t + source

(where Delta_t would be the time step, and in the assembly we multiply by Delta_t).
If you just have D instead of \dot{D} in the first term, you can write this as

D_new - \nabla \kappa \nabla D_new = source

(so in the assembly in advection.cc, you would have to remove the term with scratch.old_field_values[q] in line 393, and you would have to remove the timestep variable). So if I change the terms that are in advection.cc in your pull request, I think it should look like

data.local_rhs(i) += scratch.phi_field[i] * reaction_term/time_step * JxW;

and

data.local_matrix(i,j) += (diffusivity * (scratch.grad_phi_field[i] * scratch.grad_phi_field[j]) + (scratch.phi_field[i] * scratch.phi_field[j]) ) * JxW;

Depending on the source term (I don't quite remember what has to go in there), the right-hand sight might look slightly different. If you don't need the reaction_term, you might as well remove that, otherwise you would have to make sure the equation is never solved for a zero time step size. I think that the solver should still be able to deal with this type of equation.

That means that your equation is basically instantaneous now (the time step doesn't go in there anymore) and you are solving for \dot{D} instead of D. So what you then have in your compositional field is \dot{D}, and if you need the real D anywhere and not just its update, you would have to compute it as D_new = D_old + \dot{D} * timestep.

Does that sound like it could work for what you want to do?

So I guess that means that you need to solve a different equation compared to what we were planning to do (we really want to solve the original diffusion equation). So I would say, your first priority is probably getting this to work for the problem that you want to solve, and once this is working, we can think about how to finish the pull request. However, it might be easier for me to help if you push your updates to the branch (no need to clean them up beforehand).

In the meantime, I might make a separate pull request with the equation that we want to solve, but using the same general structure of the code (having a new advection method, a new assembler, and using the copy fields) so that once you finish up this pull request it should fit into the general structure, and would just add a new method for solving the diffusion problem. This would also mean that hopefully you won't get too many conflicts when you rebase.

[alarshi]

@jdannberg : Thank you very much!, this would work! :D . This PR is up-to date with my local branch. I will get this working and push the commit once done. Thank you again.

[alarshi]

@jdannberg: Hi Juliane. Since we have a working model of diffusion compatible with the recent ASPECT, should I close this PR?