• Registered: 28-07-2008
  • Posts: 71

Topic: Status of PrimaryField :: applyInitialCondition

I'm working through the new engineering model for the fluid-film bearing code I'm developing. I have been using nonstationarytransport and transienttransport as templates/examples. These use fields to store the solution, so I have written my code to do the same.

The problem I am implementing needs to have a nonzero initial guess, which I have been treating like an initial condition using sets. I am running into problems with PrimaryField :: applyInitialCondition. I get an error at the line:

      int indxm1 = this->resolveIndx(tStep, -1)

I think this is because time index = tStep-1 is undefined at the first step in my implementation (???)

Thus, I have two questions:

1) Is PrimaryField :: applyInitialCondition known to be working correctly? I assume it has been debugged, and the problem is with my code. However, the relevant non-zero initial condition test case input files I found look like they would be using a DOFDistributedPrmaryField, not a PrimaryField. In which case, there could be unresolved bugs in the PrimaryField code.

2) Is there a better choice for an engineering model to use as an example for a non-linear problem that requires nonzero initial conditions, and which will be solved using the NR solver?

Thanks,
Erik

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: Status of PrimaryField :: applyInitialCondition

Short answer:
1. I know several serious bugs in primaryfield. I would like to see primaryfield removed. Many of the bugs can't be fixed because existing solvers rely on strange behvaior. See below for a bit of a rant.
2. No.

I recommend using (a very recent) TransientTransportProblem as the template, but even in that code I have to put in a lot of temporary work-arounds otherwise I would break backwards compatibility with some input files.
http://www.oofem.org/forum/viewtopic.php?id=1598

Long answer;
The big problem is that a lot of code relies on inconsistent means of applying initial & boundary conditions.
These should be stored in the field, and when the elements request the values to compute strains/fluxes these should fetch all whatever is stored in the field.

This does not happen. See MasterDof :: giveUnknown

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

*************** THIS IS INCONSISTENT TIME DISCRETIZATION ************
    // first try if IC apply
    if ( tStep->giveNumber() == dofManager->giveDomain()->giveEngngModel()->giveNumberOfTimeStepWhenIcApply() ) { // step when Ic apply
        if ( this->hasIcOn(mode) ) {
            return this->giveIc()->give(mode);
        } else if ( this->hasBc(tStep) ) {
            return this->giveBcValue(mode, tStep);
        } else {
            return 0.;
        }
    }

    if ( !dofManager->giveDomain()->giveEngngModel()->requiresUnknownsDictionaryUpdate() && this->hasBc(tStep) ) {
        return this->giveBcValue(mode, tStep);
    }

    return ( dofManager->giveDomain()->giveEngngModel()->
            giveUnknownComponent(mode, tStep, dofManager->giveDomain(), this) );
}

and what I am suggesting is that it would simply change to

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) );
}

We then let the engineering model (or the field) handle the time discretization for all values. I simply store the initial condition and b.c. in the field (only reading the IC/BC once), and it is all consistently treated by the time discretization.

However, I cannot simply do these fixes, because it breaks tests that relies on inconsistent time discretization / inconsistent BC IC.
I'm trying to get nonstationarytransport and nltrasienttransport deprecated so that I can start to fix these things, but I need to check with the other developers before, but so far, noone has replied (see thread  http://www.oofem.org/forum/viewtopic.php?id=1598 )



There is also no reason for PrimaryField to exist. Changing b.c.s, adaptivity, xfem, contact b.c.s, load balancing etc., and all of these nice things will require something like DofDistributedPrimaryField anyway (or need to duplicate the functionality), so DofDistributedPrimaryField has a very real potential to actually be *faster* than PrimaryField as well (because we read many times, only store the value once per newton-iteration), though most likely cost of either field is beyond negligible.
Storing the "vectorized" solution arrays is only suitable for the brief "u=K^(-1) * f" step. Codes that has options for either PrimaryField or DofDistributedField just makes the incorrect assumption that it's costly to copy the solution vector once per iteration (it is of course a very very cheap operation).

  • Registered: 28-07-2008
  • Posts: 71

Re: Status of PrimaryField :: applyInitialCondition

Thank-you for the detailed discussion - that helps. I'll take a look at the newest TransientTransportProblem, and be aware that IC may require some extra attention to get everything to work correctly.

If I understand it, I think the solution you are proposing results in more modular code. I like the idea of classes that manage their own data as much as possible. Code always gets harder to test and understand each time a class relies on code in a different class to do something that  (conceptually) should to be done internally.

For this code development, I have been using the Google Gtest unit testing framework to test new classes somewhat in isolation as I develop them. So far it is working really well. It also seems to work well within a CMake file, which I think is required to make it practical. Once I get the engineering model test working, I'll write it up for the group to consider. A key thing to make unit testing work well, is classes that are as independent as possible.

-Erik

  • From: Gothenburg, Sweden
  • Registered: 08-09-2009
  • Posts: 1,419

Re: Status of PrimaryField :: applyInitialCondition

Actually, while I'm always in favor of more modular code, ny main reason for these transient problems was to deal with the inconsistent treatment of dof values.
This meant that when you have a simple trasient problem, you might have
v = (u_0 + u_1)/dt
but for the DOFs that have BC we instead do something like
v = u_BC * f'((t_0 + t_1)/2)
This means that the heaviside function, the derivative is always zero, so you can change u, without constant zero v.
Clearly, even though we might have the possibility to compute v "analytically" it's a bad idea. Combining this with initial conditions, turning on/off Dirichlet b.c.s, and changing Dirichlet b.c.s (for example, FE^2 simulations, the macroscale modifies the subscale b.c. repeatedly) it is clear that we cannot do it this way.