poject_solution of multiple dimension

[neosimsim]

Dear libmesh-Team

I struggle for a few days to solve the wave equation using reduction of order. So I have to solve a equation in dimension 2 for several time steps. I think I got things working pretty far but I'm stuck projecting the initial condition.

I created a SolutionFunction class

class SolutionFunciton : public FunctionBase<Number>
{
    SolutionFunciton(const unsigned int u_var) : _u_var(u_var)
    {
    }

    Number operator() (const Point &, const Real = 0)
    {
        libmesh_not_implemented();
    }
    void operator() (const Point & p, const Real time, DenseVector<Number> & output)
    {
        output.zero();
        const Real x=p(0), y=p(1);

        std::cout << "Output size: " << output.size() << std::endl; // this gives me 1
        output(_u_var) = cos(M_PI/2*x)*cos(M_PI/2*y)*cos(M_PI/sqrt(2)*time);
        output(_u_var+1) = 0;
    }
        //...
}

and call

SolutionFunciton *ic = new SolutionFunciton(system.variable_number("u")); // u is LAGRANGE_VEC
system.project_solution(ic);

in my init_cd function. Unfortunately the solution vector looks like

(x,y,z)=(0.000401774, -0.0190865, 0) 

on each node and the log of the output.size() gives we 1 every time. So I thing my system.project_solution doesn't now about u being in dimension 2.

Is there a way to project solution in such a way?

Thank you very much Alex

PS: If you need more Information like gits please feel free to ask.

[jwpeterson]

On Tue, Jul 12, 2016 at 8:13 AM, Alexander Ben Nasrallah < notifications@github.com> wrote:

Dear libmesh-Team

I struggle for a few days to solve the wave equation using reduction of order. So I have to solve a equation in dimension 2 for several time steps. I think I got things working pretty far but I'm stuck projecting the initial condition.

I created a SolutionFunction class

class SolutionFunciton : public FunctionBase { SolutionFunciton(const unsigned int u_var) : _u_var(u_var) { }

Number operator() (const Point &, const Real = 0)
{
    libmesh_not_implemented();
}
void operator() (const Point & p, const Real time, DenseVector<Number> & output)
{
    output.zero();
    const Real x=p(0), y=p(1);

    std::cout << "Output size: " << output.size() << std::endl; // this gives me 1
    output(_u_var) = cos(M_PI/2*x)*cos(M_PI/2*y)*cos(M_PI/sqrt(2)*time);
    output(_u_var+1) = 0;
}
    //...

}

and call

SolutionFunciton *ic = new SolutionFunciton(system.variable_number("u")); // u is LAGRANGE_VEC system.project_solution(ic);

in my init_cd function. Unfortunately the solution vector looks like

(x,y,z)=(0.000401774, -0.0190865, 0)

on each node and the log of the output.size() gives we 1 every time. So I thing my system.project_solution doesn't now about u being in dimension 2.

Is there a way to project solution in such a way?

I haven't used the LAGRANGE_VEC that much, but their output type is "RealGradient". From fe.h:

template<> struct FEOutputType { typedef RealGradient type; };

So I think at the very least you will need to inherit from FunctionBase. The the size of the "output" vector will still be 1, since there is only 1 variable, but you will fill up the components of the RealGradient with the necessary values.

Paul Bauman might be able to comment on this as well.

John

[pbauman]

So I think at the very least you will need to inherit from FunctionBase. The the size of the "output" vector will still be 1, since there is only 1 variable, but you will fill up the components of the RealGradient with the necessary values.

This is correct. Note, however, that not all the projection functionality is there yet for vector-valued shape functions, particularly for AMR. We do have a couple of LAGRANGE_VEC examples: here and FEMSystem-based here that may be helpful.

[neosimsim]

Thank you for your quick response.

Using RealGradient makes sense. Maybe I got the examples I read so war wrong, which were using u_var+n. I thought the n-th value of the vector is addressed by u_var+n.

I tried RealGradient now, but it seems you can't use System.project_solution() with it since it expects FunctionBase<Number>.

Is there another simple way to project initial conditions onto System.current_local_solution? I have to admit I'm still confused by dof_mapper. Until then I try to work myself through more examples.

Thanks Alex

[neosimsim]

I tried to use RealGradient and failed, since is used in several location deep inside libmesh. So I tried to project the solution myself solving a simple linear equation $Ax=b$ where $A_{i,j} = \phix_i(xj)$, $A{i+1,j} = \phiy_i(x_j)$, $b_i = solution_x(xi)$ and $b{i+1} = solution_y(x_i)$. So M is 2n_n matrix and b is of length 2_n. Now I have the problem, that PETSC is complaining

[0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Must be square matrix, rows 1024 columns 512 [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.6.2, Oct, 02, 2015 [0]PETSC ERROR: ./project_initial_condition on a arch-linux-c-opt named Turtok by neosimsim Thu Jul 21 13:53:46 2016 [0]PETSC ERROR: Configure options --prefix=/opt/petsc/linux-c-opt --PETSC_ARCH=arch-linux-c-opt --with-shared-libraries=1 --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --with-hdf5=1 --with-suitesparse=1 [0]PETSC ERROR: #1 MatGetOrdering() line 253 in /tmp/yaourt-tmp-neosimsim/aur-petsc/src/petsc-3.6.2/src/mat/order/sorder.c [0]PETSC ERROR: #2 PCSetUp_ILU() line 194 in /tmp/yaourt-tmp-neosimsim/aur-petsc/src/petsc-3.6.2/src/ksp/pc/impls/factor/ilu/ilu.c [0]PETSC ERROR: #3 PCSetUp() line 983 in /tmp/yaourt-tmp-neosimsim/aur-petsc/src/petsc-3.6.2/src/ksp/pc/interface/precon.c [0]PETSC ERROR: #4 KSPSetUp() line 332 in /tmp/yaourt-tmp-neosimsim/aur-petsc/src/petsc-3.6.2/src/ksp/ksp/interface/itfunc.c [0]PETSC ERROR: #5 KSPSolve() line 546 in /tmp/yaourt-tmp-neosimsim/aur-petsc/src/petsc-3.6.2/src/ksp/ksp/interface/itfunc.c

Is there any approach to solve non-square problems? Otherwise I would be really stuck and desperate.

Here is the link to my gist.

[pbauman]

Is there another simple way to project initial conditions onto System.current_local_solution? I have to admit I'm still confused by dof_mapper. Until then I try to work myself through more examples.

Sorry no one responded to you here. What problem are you trying to solve? Since, unfortunately, not all the needed functionality is there for LAGRANGE_VEC, you can just treat each component as it's own separate LAGRANGE variable and you should have the full functionality of the library available to you.

[neosimsim]

For now I try to solve the linear wave equation using reduction of order. So I have the position in the first and the velocity in the second component.

You mean it should be possible to add two variable "u" and "v" to the system, one for the position and one for the velocity and treat them the way?

[pbauman]

You mean it should be possible to add two variable "u" and "v" to the system, one for the position and one for the velocity and treat them the way?

That's correct. It might help to have a look at this example, where we solve the incompressible Navier-Stokes equations this way - using a different scalar variable for each component of velocity.

[neosimsim]

It seems like this approach got it working.

Thank you very much!!!

[pbauman]

Great! Please let us know if you have any other questions or problems.