From: jurak <ju...@ma...> - 2008-12-08 18:21:02
|
I am using DifferentiableSystem class for some time now, together with EulerSolver class. Recently EulerSolver::element_residual() function was changed and now it calls _system.mass_residual() only once, instead of two times, as in 0.6.2 version. This has a side-effect that EulerSolver can now solve only equations linear in time derivative, that is, of the form dx/dt = F(x), while before one could also solve d/dt M(x) = F(x), with nonlinear M. The reason is that now the function forms residual R = dt F(x_theta) + M(x_old - x) while before it would form: R = dt F(x_theta) + M(x_old) - M(x). Is there any reason for this change except optimization? I would like to be able to treat nonlinearity in time derivative term. Regards, Mladen Jurak <http://web.math.hr/%7Ejurak> |
From: Mladen J. <ju...@ma...> - 2009-01-22 17:48:20
|
Hi, 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. 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"). 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. Thanks, Mladen |
From: Roy S. <roy...@ic...> - 2009-01-22 18:03:07
|
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 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 second-order accuracy in Crank-Nicholson, 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 |
From: Mladen J. <ju...@ma...> - 2009-01-22 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 two-phase flow equations by vertex-centered 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 convection-diffusion 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 second-order accuracy in Crank-Nicholson, 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 > > |
From: Roy S. <roy...@ic...> - 2009-01-22 19:00:20
|
On Thu, 22 Jan 2009, Mladen Jurak wrote: > 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. Would it suffice to have an elem_rate_solution available, which contained (u_n-u_p)/deltat on each variable? I might be adding that for one of my own applications and to make it simpler to write user code with nonlinear mass matrices again. --- Roy |
From: Mladen J. <ju...@ma...> - 2009-01-22 19:07:34
|
Roy Stogner wrote: > > Would it suffice to have an elem_rate_solution available, which > contained (u_n-u_p)/deltat on each variable? I might be adding that > for one of my own applications and to make it simpler to write user > code with nonlinear mass matrices again. > --- > Roy > > It would be nice to have "(u_n-u_p)/deltat" or "u_p". I would prefer some "elem_old_solution" in DifferentiableSystem class. Mladen |
From: Roy S. <roy...@ic...> - 2008-12-08 18:39:38
|
On Mon, 8 Dec 2008, jurak wrote: > I am using DifferentiableSystem class for some time now, together with > EulerSolver class. Recently EulerSolver::element_residual() function > was changed and now it calls _system.mass_residual() only once, instead > of two > times, as in 0.6.2 version. > > This has a side-effect that EulerSolver can now solve only equations > linear in time > derivative, that is, of the form > dx/dt = F(x), > while before one could also solve > d/dt M(x) = F(x), > with nonlinear M. The reason is that now the function forms residual > R = dt F(x_theta) + M(x_old - x) > while before it would form: > R = dt F(x_theta) + M(x_old) - M(x). > > Is there any reason for this change except optimization? Yes; it was necessary to get correct results on moving-mesh problems where the size of an element might change between x_old and x. > I would like to be able to treat nonlinearity in time derivative > term. You might have to turn on and use the "elem_fixed_solution" which we put in for use in stabilized formulations. That'll be equal to the theta_solution in an EulerSolver run. --- Roy |
From: jurak <ju...@ma...> - 2008-12-08 19:30:25
|
Roy Stogner wrote: > > > On Mon, 8 Dec 2008, jurak wrote: > >> I am using DifferentiableSystem class for some time now, together with >> EulerSolver class. Recently EulerSolver::element_residual() function >> was changed and now it calls _system.mass_residual() only once, instead >> of two >> times, as in 0.6.2 version. >> >> This has a side-effect that EulerSolver can now solve only equations >> linear in time >> derivative, that is, of the form >> dx/dt = F(x), >> while before one could also solve >> d/dt M(x) = F(x), >> with nonlinear M. The reason is that now the function forms residual >> R = dt F(x_theta) + M(x_old - x) >> while before it would form: >> R = dt F(x_theta) + M(x_old) - M(x). >> >> Is there any reason for this change except optimization? > > Yes; it was necessary to get correct results on moving-mesh problems > where the size of an element might change between x_old and x. > >> I would like to be able to treat nonlinearity in time derivative >> term. > > You might have to turn on and use the "elem_fixed_solution" which we > put in for use in stabilized formulations. That'll be equal to the > theta_solution in an EulerSolver run. > --- > Roy > I think that I understand. By using "elem_fixed_solution" I get "x_theta = theta*x+(1-theta)*x_old", and from "elem_solution" I get " x_old - x". Then one can recover "x_old" and "x" and "mass_residual()" function calculates M(x_old) - M(x). Right? Mladen - <http://web.math.hr/%7Ejurak> |
From: Roy S. <roy...@ic...> - 2008-12-08 20:03:44
|
On Mon, 8 Dec 2008, jurak wrote: > Roy Stogner wrote: > >> On Mon, 8 Dec 2008, jurak wrote: >>> while before one could also solve >>> d/dt M(x) = F(x), >>> with nonlinear M. The reason is that now the function forms residual >>> R = dt F(x_theta) + M(x_old - x) >>> while before it would form: >>> R = dt F(x_theta) + M(x_old) - M(x). >> >> Yes; it was necessary to get correct results on moving-mesh problems >> where the size of an element might change between x_old and x. >> >>> I would like to be able to treat nonlinearity in time derivative >>> term. >> >> You might have to turn on and use the "elem_fixed_solution" which we >> put in for use in stabilized formulations. That'll be equal to the >> theta_solution in an EulerSolver run. > > I think that I understand. > > By using "elem_fixed_solution" I get "x_theta = theta*x+(1-theta)*x_old", Right. > and from "elem_solution" I get " x_old - x". Right. > Then one can recover "x_old" and "x" and "mass_residual()" > function calculates M(x_old) - M(x). I'd thought of just calculating (dM/dx)(x_theta) * (x_old-x), which is more in the spirit of the EulerSolver; calculating M(x_old)-M(x) is more like the integration that occurs in the Euler2Solver. But as long as you're using a fixed mesh, doing either (or mixing and matching) should work, and if the nonlinearities in your problem are bad enough then M(x_old)-M(x) is likely to be more robust. On moving-mesh problems, I think the robust & accurate way to theta-integrate a nonlinear mass matrix would be ((1-theta)*(dM/dx)(x_old)+(theta)*(dM/dx)(x))*(x_old-x). The Euler2Solver + FEMSystem API doesn't currently make doing that very easy... This is one reason I've been procrastinating making the FEMContext API changes that FEMSystem needs; I'd like to figure out how to really get it right this time. --- Roy |