From: Mladen Jurak <jurak@ma...>  20090122 18:40:17

Roy Stogner wrote: > Mladen Jurak wrote: > >> I wander about the use of "elem_fixed_solution" in >> "EulerSolver::element_residual" method. Fixed solution is >> set to "theta_solution" and therefore during the call to >> "_system.element_time_derivative" both "elem_fixed_solution" >> and "elem_solution" have the same data. > > Right. > >> It happens to me >> that in "element_time_derivative" I need the solution from the >> preceding time level ("old_elem_solution" in >> "EulerSolver::element_residual"). > > May I ask why? > My use of the library is not standard. I solve compressible twophase flow equations by vertexcentered finite volume method, by using some additional classes (made by David Trujillo from University of Pau). The system in fractional flow formulation is composed of pressure equation which gives the total flow, and saturation equation (of convectiondiffusion type) in which the total flow is used. The problem is on the boundary where I need to calculate total flow in order to complete saturation equation. This boundary flow can be calculated only by using local mass conservation and in that way I introduce time derivative into spatial discretization. My need for old solution comes from this time derivative, but, perhaps, I could shift these boundary terms containing (p^new p^old)/dt to mass residual part of the residual. Sorry for this confusing explanation. >> My question is couldn't "elem_fixed_solution" be set to >> "old_elem_solution" in "EulerSolver::element_residual". Then one >> would have more information while assembling the system. >> I don't see any obstacle, but in the other hand I am not sure that >> I understand well the concept of fixed solution. > > Suppose you've got a system like: > > d(M(u))/dt = F(u) > > A theta method solve, from previous timestep u_p to next timestep u_n > via theta interpolant u_t, might look like: > > (M(u_n)  M(u_p))/deltat = F(u_t) > > And when discretizing in space by taking only the weighted products > against test functions: > > (M(u_n,v)  M(u_p,v))/deltat = F(u_t,v) > > So for most people, all F() needs to evaluate is the derivative at > u_t, not u_p or u_n. > > The reason for adding elem_fixed_solution was, when we wanted to do > some stabilized formulations, we had v = v_c + w(u_f)  the test > function was now a function v(u) of the solution. But that u had to > be fixed  evaluating M(u_n,v(u_n))  M(u_p,v(u_p)) is inconsistent. > To avoid losing secondorder accuracy in CrankNicholson, the obvious > choice for EulerSolver was v(u_f) where u_f = u_t. > > But while that's the right thing to do for EulerSolver, it sounds like > it's incompatible with what you want  hence my "why?" question. The > right thing to do may be to modify a cloned EulerSolver (like I did > with Euler2Solver when I wanted trapezoidal time integration), but > it's impossible to be sure without knowing what you plan to use the > old solution for. >  > Roy > > 