Topic: Sets in TransientTransport solver

Dear developers,
I encountered inconsistency in TransientTransport solver when using sets. Let us take tmpatch21-8.in as a reference and change a solver to TransientTransport.

I can assign the same Newton BC either as "quadaxisym1ht 1 nodes 4 1 3 4 2 boundaryloads 2 2 4" or through set as "Set 3 elementboundaries 2 1 4". I obtain two slightly different results in node 1 (5.68352780e+00 vs. 5.68563857e+00) which I expect to be the same. Attached is the modified tmpatch21-8-tmp.in (un/commenting two lines gives two BC assignments).

Thank you for clarification/fix. Vit

Re: Sets in TransientTransport solver

I will look into this. I have a few possible explanations but i shoud take a closer look first.

Re: Sets in TransientTransport solver

The solution to this difference is very simple; the default error tolerance isn't working well for this particular test, or rather, the output tolerance is much stricter than the NR solver permits. There is just so little happening in the simulation that NRSolver is considered converged already, with 0 iterations, for the default rtolf or 1e-3. But, for many steps, these small errors accumulate.

Setting:
rtolf 1e-9
or
miniter 1

for the TransientTransport solver will produce the same output.

Edit: Just one thing to add, I would expect different results (depending on time-stepping) for nonlinear problems. This is intentional, as I think the nonlineartransport solver does this incorrectly (it computes time-stepping on force level, instead of on the unknown itself, and this relationship is not linear).

4 (edited by smilauer 14-09-2016 16:48:50)

Re: Sets in TransientTransport solver

Dear Mikael,
thank you for the explanation and quick fix. However, the above problem should be linear hence I never suspected iterations to be a problem. For some reason, introducing Newton's boundary condition with linear constitutive material causes iterations in TransientTransport solver (it was linear in a previous NlTransientTransportProblem solver). The unknown heat flow q=h(T-T_inf) depends on unknown T, h*T went in NlTransientTransportProblem to the left hand side and made the problem linear. It seems to me that T is now taken from previous iteration which renders the problem nonlinear or maybe the residuum is not evaluated correctly.

Attached is an input file, demonstrating this nonlinearity. If you change the solver by uncommenting 4th line, you see no iterations. Could you please have a look on that and made the problem linear again? It would speed up computation significantly.

Best regards. Vit

Post's attachments

tmquad12-error.in 1 kb, 3 downloads since 2016-09-14 

You don't have the permssions to download the attachments of this post.

Re: Sets in TransientTransport solver

Short answer: I checked the file and I only see 1 linear step (NRSolver prints the residual before the step, and after the step, to confirm that it was indeed solved, but only one linear solve is computed).

Long long answer:
NRSolver allows for 0 iterations, as long as the tolerance is satisfied. Here, the errors were small, but over many time steps, it accumulated to an appreciable difference. We should reconsider the default tolerance 1e-3, to maybe 1e-6 (it's normalized, so it should always be around this value, regardless of units) to avoid these problems.

The output, for a linear problem, would be

1. Compute the residual, print it
2. Tolerance above , compute do  u += K^(-1) * R
3. Compute the new residual, print it.
4. Tolerance is now very small, no need to compute another step.

This is precisely what I see when i run tmquad12-error (NRSolver shows "2" iterations, having converged on the second one, meaning one linear solve).

It doesn't matter at all if you move the constant part to the right hand side. It still ends up in the exact same residual  R = F_external - F_internal.
Splitting it into internal vs external forces only has a small influence in how we normalize tolerances.



Where i think nltransienttransport solver is incorrect, is it's approach to nonlinear problems.
The method "assembleAlgorithmicPartOfRhs" is a hardcoded assumption of linearity.
It's also entirely unnecessary. These forces should enter the system through the internal forces by just integration the flux over the element, not caring whether the unknowns are from boundary conditions or not.
For sufficiently small time-steps, everything gets asymptotically linear, so you can hide this error.


So, the general case (as done by TransientTransport):
1. Compute the unknown/velocities according to the time discretization, not caring where the values are from   (b.c., initial guess, unknowns)
2. Compute the flux and integrate it.
is a lot simpler, and always correct.

6

Re: Sets in TransientTransport solver

I am to reopen this post. I have tried to run the above mentioned "tmquad12-error.in"  with the last git version (ab25d255741c803bd694de9a74c934a6edc0de64) and I have obtained 13 NR iterations in the first step and 12 iterations in the second. The norm of residual is going from initial value 1.0 to 6.e-4 in the final one.

So there is something wrong. First, the problem is linear and there is no reason for iterations. So the only conclusion is that the system is not correctly linearized, tracing back to Newton BC contribution to the LHS. This BC should have contribution to LHS. I agree that in general such a  relation could be nonlinear, but it could be linearized (as other terms) and LHS contribution should be taken into account during the assembly.

Re: Sets in TransientTransport solver

If you are getting 13 iterations, then yes, something is definitely wrong. It's been quite a while since I merged with oofem.org  (been quite busy), so it might be something missing in the official version.

With my github version I only get 1 iteration (as is expected for a linear problem) and as far as I can see, I do take into account the tangent contributions from the Newton BC.

I will try to get time to merge in the changes this week (I need to test it out so nothing breaks first, since it's been quite a long time).

Re: Sets in TransientTransport solver

Sorry for there being so much of a delay on this. (I've recently switched jobs and it's a lot to do right now). I have merged in changes from oofem.org and have begun testing. Unless I discover any problem i'll push this weekend.

9

Re: Sets in TransientTransport solver

Mikael,
I did as well (created a separate branch for it called transient_tm), however, I made some changes, that I want to discuss before merging:
I have seen in transientTransport problem you use two different time steps one related to time at the end of time increment and related to intrinsic time. My suggestion, already implemented in the transient_tm branch is to use the already existing capability of timeStep to have intrinsic time as well. So my proposed modification consists in introducing new value mode type called VM_TotalIntrinsic, allowing to ask for total value evaluated in intrinsic (assembly) time, the VM_Total is now reserved to total values evaluated at target (final) time of time step, as has been implemented in
giveUnknownComponent for example.

Borek

Re: Sets in TransientTransport solver

Hi Borek,
I've considered something similar, sounds good to me.

I'm hoping to move the actual time-discretization logic into the Field class eventually.
I have merged in the changes with oofem.org, but I haven't pushed to oofem.org yet, since there were some issues to fix first (see http://www.oofem.org/forum/viewtopic.php?id=1673 )

Re: Sets in TransientTransport solver

I had a closer look to the old code here from the recent merge

https://github.com/Micket/oofem/commit/ … 2c2d539c2c

All the differences in time-stepping was intentional. I think these modifications are wrong.

+            //stepWhenIcApply.reset( new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 0, -giveDeltaT(nFirst), giveDeltaT(nFirst), 0) );//previous version for [-dt, 0]
+            stepWhenIcApply.reset( new TimeStep(giveNumberOfTimeStepWhenIcApply(), this, 0, 0., giveDeltaT(nFirst), 0) );//now go from [0, dt]

The first real time step, which starts after the time step when IC is applied should go from [0, dt], meaning the "stepWhenICApply" must be before this, ending at t = 0 and starting at t < 0.

------------------

EDIT: Ignore this and just look at the code examples in the next post. I think it's clearer to just write the code than to explain it badly


All the temporary stuff I had in DofDistributedPrimaryField and TransientTransport was also intentional. The changes in these files was a bit premature.

     if ( mode == VM_Total ) {
-        return this->alpha * val1 + (1.-this->alpha) * val0;
-    } else if ( mode == VM_Velocity ) {
-        return (val1 - val0) / tStep->giveTimeIncrement();
-    } else if ( mode == VM_Incremental ) {
-        return val1 - val0;
+        return dof->giveUnknownsDictionaryValue(tStep, mode);

where the change should have been

-    if ( mode == VM_Total ) {
+    if ( mode == VM_TotalIntrinsic ) {

All the TM elements still ask for VM_Total, so right now TransientTransport acts as if  alpha = 1.
As I mentioned in the other thread,

Time integration applies to the unknown
x = a * x^(n+1) + (1-a) * x^n

and the residual should be computed

R(x) == R(a * x^(n+1) + (1-a) * x^n))

If, and only if, R is linear, you can split that into
a * R(x^(n+1)) + (1-a) * R(x^n)
but that doesn't hold for nonlinear problems. Even if the problem is linear, there are no downsides to just treating them all as potentially nonlinear (the code is actually simpler and equally fast).

This is why I think the approach in
NLTransientTransportProblem :: assembleAlgorithmicPartOfRhs
is wrong.

--------------------

A summary how i think it should act.
1. The elements should ask for VM_TotalIntrinsic when computing the internal forces.
2. We should apply the IC and BC's to the DofDistributedPrimaryField at the start, then let the time integration code compute these themselves.
3. VM_Velocity should not be stored in DofDistributedPrimaryField. It must be a result from the time-stepping algorithm.


Expanding on 2.
We have code like this in several places in the code:

        if ( dof->hasBc(tStep) ) { // boundary condition
            val = dof->giveBcValue(VM_Total, tStep);
        } else {

but I think this must change. It doesn't respect the time integration scheme. It also breaks FE^2 simulations, since I need to change BCs, and computing the increment from the previous step simply isn't there anymore.

Re: Sets in TransientTransport solver

Sorry for being so verbose here. I think a small piece of code exemplifies what I mean with the BCs and ICs. I think the MasterDof::giveUnknown method should look like this:

double MasterDof :: giveUnknown(ValueModeType mode, TimeStep *tStep)
{
    if ( dofManager->giveParallelMode() == DofManager_null ) {
        return 0.0;
    }
    return ( dofManager->giveDomain()->giveEngngModel()->
            giveUnknownComponent(mode, tStep, dofManager->giveDomain(), this) );
}

i think giveUnknownComponent should look like this

double TransientTransportProblem :: giveUnknownComponent(ValueModeType mode, TimeStep *tStep, Domain *d, Dof *dof)
{
    return this->field->giveUnknownValue(dof, mode, tStep);
}

and I think giveUnknownValue should look like this

double
DofDistributedPrimaryField :: giveUnknownValue(Dof *dof, ValueModeType mode, TimeStep *tStep)
{
    if ( mode == VM_Total ) {
        return dof->giveUnknownsDictionaryValue(tStep, mode);
    } else {
        double val1 = dof->giveUnknownsDictionaryValue(tStep, VM_Total);
        double val0 = dof->giveUnknownsDictionaryValue(tStep->givePreviousStep(), VM_Total);
        if ( mode == VM_TotalIntrinsic ) {
            return this->alpha * val1 - (1 - this->alpha) * val0;
        } else if ( mode == VM_Velocity ) {
            return (val1 - val0) / tStep->giveTimeIncrement();
        } else if ( mode == VM_Incremental ) {
            return val1 - val0;
        } else {
            OOFEM_ERROR("Unknown value mode requested");
            return 0;
        }
    }
    return 0;
}

The identical code-path no matter if the DOF has BCs, ICs, initial guesses or not. (these are simply set in the beginning of each time step, and are not updated during newton iterations)

13

Re: Sets in TransientTransport solver

Just to inform, I have merged the transient_tm branch into master one.
I fully agree with what has been written here by Mikael on how BC and IC should be treated, this would make the code much clearer and consistent.

Re: Sets in TransientTransport solver

Hi Borek.
Good to hear we are on the same page here.

I've tried to get closer to this implementation while trying to preserve backwards compatibility, but it has been very difficult.

At this point, I think we should just make a branch, implement the code how we think it should work, see what breaks and proceed to fix those. Doing it in one big commit is just not doable.

But, even as a branch, we don't want to let this drag on, as it will probably touch every single engineering model, and probably half the tests. If development proceeds on master, merging will be a huge pain.

Any thoughts on this plan Borek?
I have a good idea on what needs to be changed, and I'm willing to work overtime to try and get this done quickly and painlessly as possible.

15

Re: Sets in TransientTransport solver

Mikael, thank you.
The new branch is ok, there will be a lot of changes on many places that have to be sorted out first before going to master branch.
I just want to note, that I made some of the changes and merged them to the master:
- introduced VM_TotalIntrinsic
- when internal forces are evaluated, the VM_TotalIntrinsic is used (at least in tm module consistently)
- the tm tests have been changed accordingly.

what was not addressed yet is the support for VM_TotalIntrinsic in other modules (FM, SM).
I agree on avoiding to store velocity in the field, especially, when it could be computed from displacements and their history.
Perhaps we can also consider to convert some engineering models to use fields and perhaps even dof distributed field in this step.

I can take care of some work, but as the end of the year is quite busy period for me, I can start working on that in mid January.
Borek

Re: Sets in TransientTransport solver

How do we treat VM_Velocity (and VM_Acceleration)?
If we have a Newmark-type time-integration, we would have stored velocities at each time-step, and a type of "VM_VelocityIntrinsic" ?
If so, then the midpoint rule would always be VM_VelocityIntrinsic (since it is computed as an average over a time-step, not at the target time)

A simple way out (without having to introduce 2 different VM_Velocity, VM_Acceleration) is to let VM_Velocity represent the "intrinsic" value in the elements (the primary field is free to call the target values whatever it wants).

Similar complications apply to VM_Incremental. If we want to know how much a force has changed from one step to the next, then the increment must be the difference in VM_TotalIntrinsic.

17

Re: Sets in TransientTransport solver

For me, the option to have both Total and TotalIntrinsic for velocity sounds much general. At the end, we want to print all quantities at the end of time step and so there is clear need for Intrinsic as well as final values. Also, I think, that VM_Total should be called VM_TotalFinal.
We can also consider having just Total, Velocity ids and another type with ids  corresponding to Intrinsic and Final.

Re: Sets in TransientTransport solver

Well, the problem arise with generalized mid point rules; there simply isn't a VM_VelocityFinal to export.

Right now, we would export
1. GP values evaluated at intrinsic time, computed for a gradient based on VM_TotalIntrinsic
2. VM_VelocityIntrinsic
3. VM_TotalFinal  (this is actually the odd one out)

I totally understand why one wants TotalFinal, but I also think exporting VM_TotalIntrinsic would be consistent with all the other output. After all, it is the values used in the actual equilibrium solution (one simply has to think of all the outputs as being "over" a time step, and not the exact value at any given time)

This is why I have been rather torn on this issue. Perhaps leaving it to the user to pick is the only pratical choice.

--------------

I agree with changing VM_Total to something clearer, VM_TotalFinal or VM_TotalTarget comes to mind (borrowing from TimeStep, though this either can probably be misunderstood)

Though, perhaps we could hold back on this just a short while? I'm trying to merge back with oofem.org and a massive rename is bound to cause tons of conflicts.

Right now parallell lb0* tests are segfaulting in my repo, and I haven't found the problem yet.