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
