dolfin team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting mspieg <mspieg@xxxxxxxxxxxxxxxxx>:

> 
> On Jun 21, 2008, at 9:50 AM, Kristian Oelgaard wrote:
> 
> > Quoting mspieg <mspieg@xxxxxxxxxxxxxxxxx>:
> >
> >> Hi all,
> >>      does anyone have an example  of how to access the coordinates
> >> and set the Dofs (quadrature points) for a user defined Function
> >> defined with a Quadrature element?  From what I can glean from the
> >> documentation, this should be possible, but I can't quite figure out
> >> how to do it.
> >
> > To set the values of you user defined function, just overload the  
> > eval()
> > function like you would for a function defined on a 'normal' element.
> > When the function is being integrated in the form it will be  
> > evaluated at
> > quadrature points, meaning that whatever coordinates are entering  
> > eval() (the x
> > argument) will be the location of the quadrature point (the 'dof').
> >
> 
> Thanks Kristian,
> 
> I figured out how to do this and it all works...the only issue is  
> that this function is a source term in a time-dependent non-linear  
> solve that is rather expensive to compute.  Right now it is being  
> evaluated every time the residual vector is being assembled, but  
> really only needs to be calculated once per time step.

So what you have is (in sloppy DOLFIN syntax):

ExpensiveFunction : public Function
{
  // Define some complex function
}

  ExpensiveFunction F_ex(mesh);

In the nonlinear problem you have:

  RHSLinearForm L(F_ex);

In main.cpp:

  loop time
  {
    // F_ex is evaluated for all Newton iterations
    newton_solver.solve(nonlinear_problem, x);
  }
 
> I'm curious to know if there is a straightforward way to evaluate it  
> once (at the quadrature points), outside of the form assembly, and  
> then just pass the tabulated function to the form within the non- 
> linear solve.  (I'm getting excellent quadratic convergence in the  
> newton solve, but this would speed things up by about a factor of 3-4).

Happy to hear that it converges.

I see too ways to do what you want (someone else might see better ways?):

Create a discrete function and for each time step compute the entries of the
underlying vector:

  Vector vec;
  Function F_ex(mesh, vec, form);

  loop time
  {
    // evaluate function according to the definition
    evaluate_complex_function(vec);

    newton_solver.solve(nonlinear_problem, x);
  }

Or you can let FFC/DOLFIN handle things for you by defining a simple form and
let the assembler take care of things for you

In FFC:
  L = v*f*dx # Functions defined on QuadratureElement

and in main.cpp:

  Vector vec;
  Function F_simple(mesh, vec, form);
  ExpensiveFunction F_ex;
  SimpleFormLinearForm L_simple(F_ex);

  loop time
  {
    // The assembler will evaluate the expensive function once for each time
    // step and store the values in the simple function
    assemble(vec, L_simple, mesh);

    newton_solver.solve(nonlinear_problem, x);
  }

  Note: Your RHS form in nonlinear problem should use the F_simple function
  which is evaluated once per time step.

  I think this should work.

Kristian

 
> All help greatly appreciated
> cheers
> marc
> >>
> 
> 
> >> ----------------------------------------------------
> >> Marc Spiegelman
> >> Lamont-Doherty Earth Observatory
> >> Dept. of Applied Physics/Applied Math
> >> Columbia University
> >> http://www.ldeo.columbia.edu/~mspieg
> >> tel: 845 704 2323 (SkypeIn)
> >> ----------------------------------------------------
> >>
> >>
> >>
> >
> 
> ----------------------------------------------------
> Marc Spiegelman
> Lamont-Doherty Earth Observatory
> Dept. of Applied Physics/Applied Math
> Columbia University
> http://www.ldeo.columbia.edu/~mspieg
> tel: 845 704 2323 (SkypeIn)
> ----------------------------------------------------
> 
> 
>