This is my workaround. I dialed the clock back and forth around
equation_systems.reinit().
I am not sure if there is a more effective way (in performance) to do so.
system.time = 0;
Create a DirichletBoundary object and add it to equation_system;
equation_systems.init ();
for each time step {
system.time += dt;
...
system.reinit_constraints();
*system.old_local_solution = *system.current_local_solution;
...
for (unsigned int r_step=0; r_step<max_r_steps; r_step++) {
system.solve();
if (r_step+1 != max_r_steps) {
mesh_refinement.uniformly_refine(1);
system.time -= dt;
equation_systems.reinit ();
system.time += dt;
system.reinit_constraints();
}
}
}
--Junchao Zhang
On Wed, Jan 27, 2016 at 4:03 PM, Junchao Zhang <junchao.zh...@gmail.com>
wrote:
> Thanks for the info. Let me try the workaround.
> By the way, in the last few hours, I debugged the code and found these
> comments in TransientSystem<Base>::reinit (). I traced through the class
> hierarchy from TransientSystem, LinearImplicitSystem, ImplicitSystem,
> ExplicitSystem to System, but found nothing related to old solution
> projection.
>
> template <class Base>
> void TransientSystem<Base>::reinit ()
> {
> // initialize parent data
> Base::reinit();
>
> // Project the old & older vectors to the new mesh
> // The System::reinit handles this now
> // this->project_vector (*old_local_solution);
> // this->project_vector (*older_local_solution);
> }
>
> --Junchao Zhang
>
> On Wed, Jan 27, 2016 at 3:43 PM, Roy Stogner <royst...@ices.utexas.edu>
> wrote:
>
>>
>> On Wed, 27 Jan 2016, Junchao Zhang wrote:
>>
>> I slightly modified adaptivity_ex2.C ( see the outline below). I find
>>> after the mesh is refined, the old solution is not correctly projected to
>>> the new mesh.
>>> I print out old solution in each refinement step. At the beginning, time
>>> = 0.025, r_step = 0, old_local_solution->print() prints the solution at
>>> time = 0, which is correct. However, in next refinement (i.e, r_step =
>>> 1), old_local_solution->print()'s output has boundary values at
>>> time=0.025,
>>> which instead should be values at time=0, I think. The wrong
>>> old_local_solution will affect matrix assembly, where we call
>>> system.old_solution()
>>> to compute u_old and grad_u_old. I tried to debug the code, but it is
>>> hard
>>> to get how equation_systems.reinit () works in libmesh.
>>> Is there something missing in my code? Thank you!
>>>
>>
>> Nothing missing in your code, just a problem with the way we handle
>> Dirichlet boundary conditions, I'm afraid. You might be the first
>> person to combine time-varying Dirichlet conditions with transient AMR
>> and pay close enough attention to notice an O(h*deltat) flaw in the
>> results.
>>
>> Dirichlet boundary conditions get converted into heterogeneous
>> constraint equations, and we only store one set of constraint
>> equations at a time, but we assume that these equations are applicable
>> to all vectors in need of constraint. Transient AMR breaks that
>> assumption.
>>
>> We do have an option to forgo projecting a particular vector, but I
>> don't think we have an option to project a vector without also
>> applying constraint equations afterward... and we can't just add that
>> option, because often applying constraint equations is necessary, so
>> that e.g. an adaptive coarsening doesn't leave a non-conforming value
>> on a hanging node. Even applying only hanging node constraints
>> wouldn't be an option, because we expand our constraint equations, and
>> that expansion process brings in Dirichlet values in the cases of
>> hanging nodes with neighboring nodes on a Dirichlet boundary.
>>
>> Gah.
>>
>> We currently do allow for each adjoint solution vector to have its own
>> independent set of heterogeneous constraint values. Perhaps we could
>> extend that to any solution vector, old_solution included? We'd also
>> then need to re-evaluate the Dirichlet function after mesh refinement.
>> None of that would be easily fixable from user code without tricky
>> library changes.
>>
>> That's nothing I'd be able to get to myself in the immediate future,
>> although obviously "subtle inaccuracy when mixing time-dependent
>> Dirichlet + AMR" ought to be a high priority to fix.
>>
>> As a workaround from user code: you could project while using a
>> "time lagged" boundary function, then rerun reinit_constraints() with
>> the correct function, then enforce_constraints_exactly() on the
>> solution, which will still leave the old_solution in the old state.
>> ---
>> Roy
>>
>
>
------------------------------------------------------------------------------
Site24x7 APM Insight: Get Deep Visibility into Application Performance
APM + Mobile APM + RUM: Monitor 3 App instances at just $35/Month
Monitor end-to-end web transactions and take corrective actions now
Troubleshoot faster and improve end-user experience. Signup Now!
http://pubads.g.doubleclick.net/gampad/clk?id=267308311&iu=/4140
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
