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.