Field on node and its evaluator

[RuixiongHu]

Hello all,

This is Ruixiong Hu, grad student from RPI, I'm taking over the 3D Manufacturing (AMP) project based on Albany. This is a pure 3D thermal conduction problem. I've been working on understanding the code for a while and could understand the general structure and most of tasks regarding parameter and its evaluator.

Currently I'm working on storing a historical max Temperature field on node point. Since from our observation parameters are all based and stored on quadrature point, we decided to create another solution field Tmax, to store this historical value. So I created a new DOF Tmax in the same way of Temperature. This is just a field without residual, which is slightly different than the template problems so I have some questions about how to implement the idea:

1. Does it have to be accompanied with a residual field? Since this Tmax will not be calculated by residual from governing equation, but merely compare the max value of the primary solution Temperature. It shows error of missing Tmax residual evaluator .

2. I would like to know where and how to write a evaluator to achieve this max(T,Tmax) function. Usually parameters are calculated inside evaluators by single equation, and DOF field is calculated through residual. But how can I calculate this DOF field through a evaluator type of file with one function, instead of getting it by solving another common residual term.

I'm a relative new user of Albany and also programming so I would really appreciate any type of help. Please let me know if anything above is wrong or if there's any other easier way already existed in Albany.

Thanks and have a good one
Ruixiong

[bartgol]

Hi Ruixiong,

If I understand correctly, you have a field that you compute at the quadrature points, but you would like to have one version of it at the nodes, and save it, correct?

I think you have a few choices:

1. Save the field as a P0 field (i.e., one value per element): you compute the cell centered field with a QuadPointsToCellInterpolation evaluator.
2. You can save a quadrature point field, actually. It is not that pretty when you open it in paraview, since it gets stored as a vector field (1st component is the 1st quadrature point, 2nd component is the 2nd quadrature point, and so on). This is probably the worst solution.
3. Look at how you compute the field on the QuadPoints, and recreate the same chain of evaluation at the nodes. This is doable only if you are not computing derivatives (otherwise you're stuck at quad points). We do this for some fields in LandIce. E.g., the evaluator BasalFrictionCoefficient is created in StokesFOBase twice, once to compute the output at quad points, and once to compute it at the nodes (for I/O purposes).
4. If you really want the field at the nodes, and your field involves derivatives, you need to solve P x_n = x_q, where P is the projection matrix from nodes to quad points. If the number of quad points is larger than the number of nodes, then this can be solved with the normal equations, and it's even solvable locally, meaning that two bordering elements should be able to recover the same values for the shared nodes, even if the problem is solved locally, element by element. If the number of qp is lower than the number of elements, you need to do some sort of pseudo-inverse, but that will not yield the same value for nodes shared between elements. As a result, the last element to be processed will determine the value that is saved for that node.

In all these cases, you will also need a SaveStateField evaluator.

[RuixiongHu]

Thanks for the detailed explaination!

I think the 3rd option might be the one I wanna check out.

Just to be more specific, I'm looking for the historical max value of my primary solution. Since temperature may go up and drop periodically, but I want to record the peak value of all time, and use this node based value later in some evaluators. So this involves no residual term and cannot be solved with exactly the same governing equation as in Option #3. So I was wondering how to implement this max( ) function.

And for other options, I need different values on different vertices, since I would like to interpolate those nodal values on to edge to do calculation later.

Please let me know if anything new come up. And meanwhile I'll check the StokesFOBase

Thanks and have a nice day!
Ruixiong

[bartgol]

If that's what you need, then you can implement your evaluator easily. Say you call it HistoricalMaxValue, then you would need two things:

1. in the postRegistrationSetup, initialize the values of the MDField to std::numeric_limits<Real>::min()
2. in the evaluateFields method, simply do fieldMax(cell,node) = max(field(cell,node),fieldMax(cell,node)

Note: I'm assuming you're not really interested in the max of the field during Jacobian/Tangent/DistParamDeriv evaluation types, so you should probably have an empty implementation, and only have an implementation for the Residual case. Much like SaveStateField does.

P.s: FYI, StokesFOBase is not the easiest class to read, to be fair. I did it since we were starting to get a lot of code duplication, hard to maintain. But it may be not so easy to read.

[bartgol]

Also, if/when you think your question has been answered, don't forget to close the issue. Thanks.

[RuixiongHu]

Hi thanks a lot for the suggestion!

I tried to implement it in this way, something's still little bit wrong, so can I ask a few questions? Here's how I did it:

In problem.hpp:

• Initialize dof_names Tmax. (Tmax residual not declared since it will give error of not having evaluators for this residual)
• Did not include a seperate Tmax block registerStateVartiable and SaveStateField, since I think this is to initialize Tmax onto integration points.
• Create ParameterList and registerEvaluator for Tmax at the end of problem.hpp file, like the regualr residual parameter list. (In this place not calling SaveStateField since it will somehow give error of no 'Field Name' in this parameter list. Not calling registerStateVariable either (since it's a field instead of a state variable) )

In Tmax.hpp:

• Initialize Tmax.hpp file as how other evaluators do for Parameters. Did not follow the way of residual.hpp (since I think this is not a residual implementation).

In Tmax_Def.hpp:

• Get T_, Tmax_ as dl->node_scalar
• Follow the way of other parameter evaluators (grab old values, do postRegistrationSetup)
• Then loop over nodes to do max( T_(cell,node), Tmax_old(cell,node) )

Above is generally what I added so far, it is compiled but not running successfully.
It does not enter the Tmax_ calculation since I did an output in evaluateField.
And it does not show initializing state: Tmax (this is maybe due to I'm not calling it stateVariable), but it clearly shows it allocated Tmax to be node_scalar dependent field, so may I ask a few questions?

1. I tried to skip some part as stated above based both on my understanding and output error , is there anything that I should add back?
2. Can I know how to declare the MDField to std::numeric_limits<Real>::min() in postRegistrationSetup? Currently I'm doing postRegistrationSetup(typename Traits::SetupData d, PHX::FieldManagere<Traits>& fm)
3. Is this the correct way to skip the residual evaluator but calculate the dof field the same way as other parameters?

I do apologize to ask questions in this detail, but this is all I could figure out after repetitively examining the code and similar examples.

I'll really appreciate the help, thanks!
Ruixiong

[ikalash]

@RuixiongHu : option 4 described by @bartgol sounds like the ProjectIPtoNodalField response in LCM. Basically what this response does is take fields at quadrature points (stresses, equivalent plastic strain, etc.) and projects them onto the nodes. This is done by constructing a mass-like matrix and solving a linear system. The attached paper details the method. You may want to have a look to see if you can adapt the response to what you want to do.

LieGroupLieAlg.pdf

[RuixiongHu]

Thanks for the information @ikalash I I'll learn about this and try it if I need further implementation. Currently I think the last idea by @bartgol should be enough to address my current need, and I'm trying to get that done and then see what comes later.

Thanks!

[bartgol]

@RuixiongHu I think Albany lingo can be a bit confusing. A few comments:

1. If you want to save/load a field to/from the mesh, you do need to add the field as a state (with registerStateVariable). A state here is not referring to the dyanmical system sense (i.e., x in xdot = f(x,t)). A state can be any field that you compute. And it does not need to be on the quadrature points. Any field can be a state.

2. PHX does not evaluate a field (ie., run its evaluator) unless it is asked to, or the field is directly or indirectly required by a requested field. If you look at any problem, you will see that we request the evaluation of the scatter residual operation. This, in a cascade fashion, forces the evaluation of all the evaluators that are needed to compute the residual. However, saving a field is not necessary for the evaluation of the residual, Hence, you need to explicitly ask for its evaluation. To do this, you can do something like this

ev = rcp(new PHAL::SaveStateField(....));
if (ev->evaluatedFields().size()>0) {
   fm.template requireField<EvalT>(*ev->evaluatedFields()[0]);
}

This is telling PHX that, if the evaluator evaluates something (which is true only if EvalT=PHAL::Residual), then it should make sure that that evaluator is evaluated.

3. to init a field to a specific value (like std::numeric_limits<ST>::min()), the correct place is indeed inside postRegistrationSetup, as you already inferred.

So in short, this is the recipe:

• register the nodal field you want to save
• create the SaveStateField with the parameter list that is returned by the registerStateVariable call
• tell PHX to evaluate the save operation, with the syntax I wrote above
• make sure there's an evaluator that computes the nodal field (or else PHX will complain about a missing evaluator)
• enjoy life.

Extra Note: if you are using a STK mesh, then you also need to declare the field in the Required Fields Info sublist in the discretization section. Take a look to the tests/small/LandIce/Hydrology/input_steady.yaml to see how you should declare your field, and mark it as an output field. Again, this is needed only for STK meshes (so far).

[RuixiongHu]

@bartgol Thanks for the explanation!

I've implemented all the above suggestion. It's running but still falls into the same problem, I observed that:

It enters all other evaluators perfectly, but never output the message that it enters the new evaluator, until it shows segmentation fault.

I'm sure that the code is running perfectly before I add this new things in. So here's my thought :

1. The new evaluator still has something wrong for it to be called iteratively like other parameter evaluators, is this because residual does not requires this field so it never gets called? I used these lines at the end of registering the state variable to be node_scalar:

	ev = Teuchos::rcp(new PHAL::SaveStateField<EvalT,PHAL::AlbanyTraits>(*p));
	if (ev->evaluatedFields().size()>0){
	  fm0.template requireField<EvalT>(*ev->evaluatedFields()[0]);
	}
	fm0.template registerEvaluator<EvalT>(ev);

2. For the seg fault, if this really happens after the new evaluators are called. Then the only mistake I would imagine is that inside the new evaluator files, when I was trying to load Temperature (primary solution) to be node scalar. Temperature_ p.get<std::string>("Temperature Name"), dl->node_scalar),
   But I only see this Temperature to be expelicitly registered as qp_scalar.

However, for my understanding, even though calculation is done based on integration point, but Temperature as a DOF will automatically have value on node point. Is this correct? If not (meaning I have to explicitly declare DOF on node like before), while registering primary solution to be also on node, how could I pass the correct calculation result onto it, given that it does not know this nodal field called Temperature is exactly the DOF I'm solving, since therer's already a Temperature on qp point that is enough to complete a simulation.

I'll really appreciate your answer, thanks!
Ruixiong

[bartgol]

	ev = Teuchos::rcp(new PHAL::SaveStateField<EvalT,PHAL::AlbanyTraits>(*p));
	if (ev->evaluatedFields().size()>0){
	  fm0.template requireField<EvalT>(*ev->evaluatedFields()[0]);
	}
	fm0.template registerEvaluator<EvalT>(ev);

Maybe it's nothing but...have you tried register the evaluator before the if statement? I'm not sure if requesting a field that has not yet an evaluator that evaluates it leads to the correct result. I would register the evaluator first, then request for its field.

[bartgol]

Regarding point 2, if your solution field is called 'Temperature', then there should be already an evaluator for that.

As a check, you should add Phalanx Graph Visualization Detail: 4 to the Problem sublist in your input file. This will make Phalanx save a DAG of the evaluation chain to file. You can then process this dag with the dot command (Albany outputs the command you need to run), which then produces a nice figure that helps you to visually check if all the dependencies are met. They should be, since PHX is not erroring out. However, the graph should show you whether the evaluator that saves the field is indeed created or not.

[RuixiongHu]

Thanks for the reply!

Yes I did register the evaluator before the if statement, and this registerEvaluator is for SaveStateField(), I did registerEvaluator twice because all other blocks are calling this after calling SaveStateField(), or otherwise it will show error of not having evaluated field Tmax(dummy) inside this evaluator.

Here is the full block:

 	RCP<ParameterList> p = rcp(new ParameterList("Tmax"));

 	//Input
 	p->set<string>("Temperature Name","Temperature");

 	//Output
 	p->set<string>("Tmax Name","Tmax");

	// Register the name of the evaluator
 	ev =rcp(new TDM::Tmax<EvalT,AlbanyTraits>(*p,dl_));
 	fm0.template registerEvaluator<EvalT>(ev);

	// Register state variable
	p = stateMgr.registerStateVariable("Tmax", dl_->node_scalar,dl_->dummy, eb_name, "scalar", 0.0, true);

	// Save state field
	ev = Teuchos::rcp(new PHAL::SaveStateField<EvalT,PHAL::AlbanyTraits>(*p));
	if (ev->evaluatedFields().size()>0){
	  fm0.template requireField<EvalT>(*ev->evaluatedFields()[0]);
	}
	fm0.template registerEvaluator<EvalT>(ev);

If everything seems correct I'll use the Graph Visualization function. Thanks for all your detailed help!

Have a good one

[bartgol]

Ah, it just occurred to me, that SaveStateField by default assumes the field to be cell centered. Try to add

p->set<bool>("Nodal State", true);

right after registering the state (and before creating the evaluator). It's a poor implementation on my side (I added that), since the PHX layout (with its tags) should already be enough to determine what kind of field it is.

[RuixiongHu]

Thanks!

I've added this, but it shows me error:

Error! Save nodal states available only for stk meshes.

And while I was searching I figured out another example doing this nodal field, it also always has a
p->set("Field Layout", dl->node_scalar);
Do I need this for my implementation?

Thanks!
Ruixiong

[bartgol]

Ah! What kind of discretization are you using?

[RuixiongHu]

I think it's PUMI? The discretization software is Simmetrix and the partitiond mesh is in smb format.

[bartgol]

In that case, I think I'm not going to be of great help. Probably the best people to assist you would be @ibaned or @bgranzow , since they wrote most of the pumi/apf stuff. Do you guys have any idea how one should proceed to save a generic node field to the mesh?

[RuixiongHu]

Hello, hope you guys had a great weekend!

During this time, except for trying to add PUMI support for the nodal field, I also tried to use a brief STK mesh to test the evaluators I added so far, and I've figured out some questions that I would really appreciate some answers.

1. The STK mesh I used is from test/small/LCM/Necking3DSTKAdapt, and in 'Discretization' block I copied the corresponding block and changed 'Exodus Solution Name' and 'Exodus Residual Name'. And the example has no 'required fields' term.

I would like to know if doing this is reasonable to test the evaluators I added so far?

I also see the examples in LandIce/Hydrology/input_steady.yaml. But in that case the discretization method is Gmsh, and input terms are somewhat different. So I'm bit confused which one to use as examples.

2. While running it ,at start of step 0, it gives me error as below, is there anything I can do by this error?

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Start of Continuation Step 0 : Attempting to converge initial guess at initial parameter values.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Phalanx writing graphviz file for graph of Residual fill (detail =4)
 Process using 'dot -Tpng -O phalanx_graph' 
 

p=0: *** Caught standard std::exception of type 'std::runtime_error' :

 
/lore/hur3/Albany_314/src/evaluators/state/PHAL_SaveStateField_Def.hpp:231:
 
 Throw number = 1
 
 Throw test that evaluated to true: scalar_field==0
 
 Error! Field not found.
 
*** Teuchos::StackedTimer::report() - Remainder for a block will be ***
*** incorrect if a timer in the block does not exist on every rank  ***
*** of the MPI Communicator.                                        ***
Albany Stacked Timer: 0.270191 [1]
Albany: Total Time: 0.269382 [1]
|   Albany: Setup Time: 0.249354 - 92.5652% [1]
|   Thyra::DefaultModelEvaluatorWithSolveFactory<double>::evalModel(...): 0.00522776 - 1.94065% [1]
|   |   Albany: **Total Fill Time**: 0.00500972 - 95.8292% [1]
|   |   |   > Albany Fill: Residual: 0.00492828 - 98.3744% [1]
|   |   |   |   Phalanx::SortAndOrderEvaluators: 1.4732e-05 - 0.298928% [1]
|   |   |   |   Phalanx: Evaluator 0: [<Residual>] Gather Solution<Residual>: 5.0653e-05 - 1.0278% [1]
|   |   |   |   Phalanx: Evaluator 14: [<Residual>] Save Field Tmax to State TmaxResidual: 0.00196367 - 39.8448% [1]
|   |   |   |   Remainder: 0.00289923 - 58.8284%
|   |   |   Remainder: 8.1436e-05 - 1.62556%
|   |   Remainder: 0.000218037 - 4.17076%
|   Remainder: 0.0148002 - 5.49411%

3. As advised by Luca I've created a DAG Graph Visualization and installed Graphviz in my machine, by executing the line dot -Tpng -O phalanx_graph, it tells me "png not recognized". However among the options it provides, there's no one could open up the graph created. Is there anything I should do to open this graph?

I would really appreciate any type of reply,
thanks!
Ruixiong

[bartgol]

Hi Ruixiong,

in reverse order:

3. check the name of the generated dag file. I think it should be phalanx_graph_0. The phalanx_graph part is just the root, Albany adds an underscore plus the iteration number.
4. the error means that you did not add the file to the mesh. You need to declare a field in the 'Required Fields' section, and set 'Field Usage' to 'Output'. The hydrology test should show you how. The fact that it is a gmsh mesh instead of STK3D is not a big deal. They are both stk discretizations, they only differ on how they create the mesh (gmsh loads from gmsh file, stk3d builds a structured mesh from scratch).
5. I have never used 'exodus solution name' or 'exodus residual name'. I think those simply allow you to give a different name to the solution/residual than the one albany would give by default.

[RuixiongHu]

Hi Luca,

Thanks a lot for you reply.

• For the DAG file part, I'm actually executing the correct file name, so I guess this should not be the error. After I do a dot -Tpng -O phalanx_graph_0 It just tells me :
  Format: "png" not recognized. Use one of: canon cmap cmapx cmapx_np dot dot_json eps fig gv imap imap_np ismap json json0 mp pic plain plain-ext pov ps ps2 svg svgz tk vml vmlz xdot xdot1.2 xdot1.4 xdot_json

• As for the 'Required Field' part, I included all the state variable that calls SaveStateField() in the problem.hpp file, and fixed most of the error. But finally it still gives me the same 'Field not found' error.
  In the 'Required Field' block I'm not adding terms that didn't call SaveStateField() and registerStateVariable, such as Residual, even though they appears as intermediate terms between iterations. Is this correct?
  What I'm guess is in ' Required Field Info' I did not add 'Time' term, which is also registered and saved as state variable, but somewhat differently, saved as a 'workset_scalar'. But when I added this, it tells me this "the field 'Time' was not found in the mesh".
  So I was wondering if this term is the cause, or there could be others

Thanks and hope to hear from you
RuixiongHu

[bartgol]

Re: Dag. That's odd. What version of dot are you using? I have 2.30 and it works.

Re: field not found in the mesh. The 'Required Field Info' section is only for element node fields (scalar, vector, or tensor), so don't add Time to that section. I don't think you need to add the residual to that list, so that's ok. Can you please post here the discretization section of your input file?

[RuixiongHu]

Thanks for the reply!

Ah, should I use webdot for that dot command? or graphviz, I thought it's graphviz 2.4, but I see there's a webdot 2.3 under it.

Here's thee discretization section

` Discretization:

Method: Exodus
Workset Size: 50
Exodus Input File Name: smooth_tension_coarse_tet4.exo
Exodus Output File Name: smooth_tension_coarse_tet4_out.exo
Use Serial Mesh: true
Required Fields Info:
  Number Of Fields: 7
  Field 0:
Field Name: Temperature
    Field Type: Elem Scalar
    Field Usage: Output
  Field 1:
Field Name: Tmax
    Field Type: Elem Scalar
    Field Usage: Output
  Field 2:
Field Name: Phi1
    Field Type: Elem Scalar
    Field Usage: Output
  Field 3:
Field Name: Phi2
    Field Type: Elem Scalar
    Field Usage: Output
  Field 4:
Field Name: Psi1
    Field Type: Elem Scalar
    Field Usage: Output
  Field 5:
Field Name: Psi2
    Field Type: Elem Scalar
    Field Usage: Output
  Field 6:
Field Name: Laser Source
    Field Type: Elem Scalar
    Field Usage: Output

`

I wasn't using exactly the hydrology mesh is because that's a 2D mesh. This is a mesh in /test/small/Necking3DSTKAdapt

All the terms here, Temperature and Tmax aree DOF, while others are state variable. Type wise, only Tmax is explicitly initialized as node_scalar, others are all initialized as qp_scalar.

These are pretty much all the fields in the simulation. I checked in the problem.hpp file which evaluator called registerStateVariable(), and also those are only terms explicitly initialized by stateManager at starting of the run, and it's not complaining about any wrong term in the required section that is not defined in the simulation. There are some parameter calculated inside the simulation such as Conductivity, but they're never registered or saved so I did not include them.

I used Field type : elem scalar because otherwise it says field is in rank 2, my required is rank1.

Thanks a lot!

[bartgol]

graphviz is what I use. 2.4 should be plenty good. I don't understand why it gives you that error.

Meanwhile, a possible clarification: you cannot save to the mesh a qp_scalar field. Well, technically you can, but it's messy and not much useful for visualization. If you specify 'Elem Scalar', then the field must have layout cell_scalar2 in Albany; in other words, it should be one value per element.

More precisely, these are the field layouts that corresponds to the 'Field Type' specs in the input file:

• Elem Scalar: cell_scalar2, that is, MDALayout
• Elem Vector: cell_vector, that is, MDALayout<Cell,Dim>
• Node Scalar: node_scalar, that is, MDALayout<Cell,Node>
• Node Vector: node_vector, that is, MDALayout<Cell,Node,Dim>

There is also support for layered fields, but that's not relevant in your case, I think.

So to summarize: if you specify, say, an Elem Scalar output field, you must have an evaluator that evaluates a field with cell_scalar2 data layout.