Here is a possible suggestion:
for example, the current NR subroutines for equilibrium and phase field variable are:
[STIFFNESS_MATRIX, strain_en_undgr, displ, res_0, res, history_vars_new] = newton_raphson(sys, D_RHS_tract, p_field, displ, label_el, history_vars_old, args)
[p_field, res_0, res, history_vars_new] = newton_raphson( sys, irrev, strain_en_undgr, p_field, p_field_old, label_el, penalty_irrev, penalty_recov, history_vars_old, args )
and we are thinking about a generalized NR similar to
[residual, Field_vars_new, history_vars_new, STIFFNESS_MATRIX (~), strain_en_undgr(~) ]= newton_raphson(sys, Mesh, Fieldvars, Field_vars_old, History_vars, Assembly_subroutine_string, sol_params)
where for outputs:
residual
is a cell or structure or array and it has res
and res_0
Field_vars_new
is a cell containing the new displacement in position 1, phase field variable in position 2, and possibly other fields like concentration.
history_vars_new
is the cell containing the new history variables
STIFFNESS_MATRIX
(~)
, strain_en_undgr``(~)
could be called or not in the output depending on the solved problem
As for the inputs:
We insert label_el
in Mesh
(Or maybe use Mesh%label_el
)
Fieldvars
is the cell containing the vectors solutions: displacement (1), phase field variable (2) , and other possible fields. The positions are hard coded in the called assembly subroutine; for example if we solve only for the displacement, the assembly will call disp =Fieldvars{1}, while if we solve for the equilibrium and phase field, the corresponding assembly subroutine will call disp =Fieldvars{1}, pf =Fieldvars{2}
Field_vars_old
:same fields in a similar cell at previous time step
history_vars_new
is the cell containing the new history variables
sol_params
will be structure adapted to the solved problem; For equilibrium, it will have only D_RHS_tract
for equilibrium as field. While for Phase field it will have only irrev
, undgr
, penalty_irrev
, penalty_recov
for equilibrium as fields; we can think of creating a global sol_params_global
where sol_params_global.disp
is called in the NR solving for equilibrium, and sol_params_global.pf
called when solving for the phase field.
It could be nice to have it very soon, since my master thesis student and I are currently trying to add my Diffusion Code into Griphith. @pcarrara what do you think?
1 - After discussion with @jheinzmann , we think that using arrays could be challenging since each field have its own different number of degrees of freedom (for example, in a coupled diffusion problem in three dimensions, we will have 3 degrees of freedom per node for displacement, 2 degrees of freedom per node for Diffusion, and 1 degree of freedom for damage ). Same applies for history variables. 2 - It could be also a nice idea to generalize the NR routine, which means that we need include flags in the Input file to create the assembly subroutine string. However, there will be an issue regarding the inputs of NR as they are different from one problem to another.
@jheinzmann , I have just discussed the issue of projecting "generalized" stress from GP to nodes this morning with Pietro and we still don't know exactly how to handle it. We have 2 ideas that I will try locally and hopefully before next Wednesday.
Thanks for the feedback.
In the near future, new unknowns solutions need to be accounted for in addition to "displ" and "p_field", such as "Concentration" when dealing with Diffusion coupled problem.
I don't know if this needs to introduce a structure variable called "Solution", such that we no longer modify the number of inputs in phase_field.fem.solver.equilibrium.newton_raphson
and phase_field.fem.solver.pf.newton_raphson
; depending on the problem studied (for example elasticity and fracture only), one could call his solution vectors by Solution.displ and Solution.pf.
The Newton Raphson subroutine will look like this
function [STIFFNESS_MATRIX, strain_en_undgr, Solution, res_0, res] = newton_raphson(... sys, D_RHS_tract, Solution, label_el,... args... )
instead of
function [STIFFNESS_MATRIX, strain_en_undgr, displ, res_0, res] = newton_raphson(... sys, D_RHS_tract, p_field, displ, label_el,... args... )
Obviously, this new change means also modifying how the unknowns are called in every fortran subroutine...
Hi, I don't know if this could help but: https://discourse.paraview.org/t/vtk-separate-mesh-and-data-files/10134 it suggests to look at .pvtu extension instead of .vtk, such that you can call the mesh from an independent file. The details are here : https://kitware.github.io/vtk-examples/site/VTKFileFormats/ I didn't take a deep look at it yet. Best, Hamza