Assembly of system matrix for schrödinger equation with complex numbers

[Maxi Miller]

I have the following schrödinger equation:

which can be split up into the real and the imaginary part (as it is done in example 29):

This again can be solved using the RK-method by splitting the current and the future part up:

Rewriting it into the weak form results in

Thus I can create my right and my left side for the calculations in dealII as

This is done in code as (using FESystem for real and complex values):

        const FEValuesExtractors::Vector real_values(real_part), imag_values(imag_part);
        QGauss<dim>   quadrature_formula(fe.degree+2);
        FEValues<dim> fe_values (fe, quadrature_formula,
                                 update_values    | update_gradients |
                                 update_quadrature_points  | update_JxW_values);

        const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
        const unsigned int   n_q_points      = quadrature_formula.size();

        FullMatrix<double>   local_matrix_left (dofs_per_cell, dofs_per_cell), local_matrix_right(dofs_per_cell, dofs_per_cell);
        Vector<double>       local_rhs (dofs_per_cell);

        std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
        std::vector<double> rhs_values (n_q_points);
        Vector<double>        rhs_values_vector(n_q_points);

        std::vector<Tensor<1, dim>> fe_values_real_values(dofs_per_cell), fe_values_imag_values(dofs_per_cell);
        std::vector<Tensor<2, dim>> fe_values_real_gradients(dofs_per_cell), fe_values_imag_gradients(dofs_per_cell);


        for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end(); ++cell)
        {
            fe_values.reinit (cell);
            local_matrix_left = 0, local_matrix_right = 0;
            local_rhs = 0;
            fe_values.get_function_values(old_solution, rhs_values);


            for (unsigned int q=0; q<n_q_points; ++q)
            {
                for(size_t i = 0; i < dofs_per_cell; ++i)
                {
                    fe_values_real_values[i] = fe_values[real_values].value(i, q);
                    fe_values_imag_values[i] = fe_values[imag_values].value(i, q);
                    fe_values_real_gradients[i] = fe_values[real_values].gradient(i, q);
                    fe_values_imag_gradients[i] = fe_values[imag_values].gradient(i, q);
                }
                for (unsigned int i=0; i<dofs_per_cell; ++i)
                {
                    for (unsigned int j=0; j<dofs_per_cell; ++j)
                    {
                        local_matrix_left(i, j) += (fe_values_real_values[i]*fe_values_real_values[j]/time_step - 1/(4*wavenumber)*scalar_product(fe_values_imag_gradients[i], fe_values_imag_gradients[j]) + fe_values_imag_values[i]*fe_values_imag_values[j]/time_step + 1/(4*wavenumber)*scalar_product(fe_values_real_gradients[i], fe_values_real_gradients[j]))
                                * fe_values.JxW(q);
                        local_matrix_right(i, j) += (fe_values_real_values[i]*fe_values_real_values[j]/time_step + 1/(4*wavenumber)*scalar_product(fe_values_imag_gradients[i], fe_values_imag_gradients[j]) + fe_values_imag_values[i]*fe_values_imag_values[j]/time_step - 1/(4*wavenumber)*scalar_product(fe_values_real_gradients[i], fe_values_real_gradients[j]))
                                * fe_values.JxW(q);
                    }
                }
                rhs_values_vector[q] = rhs_values[q];
            }
            local_matrix_right.vmult(local_rhs, rhs_values_vector);
            cell->get_dof_indices (local_dof_indices);
            constraints.distribute_local_to_global(local_matrix_left, local_rhs, local_dof_indices, system_matrix, system_rhs);
        }

But now when solving that system using

        SolverControl           solver_control (10000, 1e-12*system_rhs.l2_norm());
        SolverCG<>              cg (solver_control);

        PreconditionSSOR<> preconditioner;
        preconditioner.initialize(system_matrix, 1.2);

        cg.solve (system_matrix, solution, system_rhs,
                  preconditioner);
        //solver_control.check(0, 0);

        std::cout << "   Needed steps for equation: " << solver_control.last_step()
                  << " CG iterations."
                  << std::endl;

and Initial/boundary conditions set to 0 my result looks like

As far as I am concerned, that looks wrong, but at which step did I make a mistake here?

[Wolfgang Bangerth]

On 08/21/2017 06:51 AM, 'Maxi Miller' via deal.II User Group wrote:
>
> As far as I am concerned, that looks wrong, but at which step did I make
> a mistake here?

Hard to tell -- you'll have to go through the problem step by step and
check every part.

But I do see two problems already:
* Your bilinear form (lhs) only couples phi and phi, and psi and psi. In
other words, there is no coupling between real and imaginary parts on
the left side, but that contradicts the strong form you previously showed.

* The matrix is not SPD. So you can't use SolverCG.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Maxi Miller]

Thank you for the answer! According to your statement

* Your bilinear form (lhs) only couples phi and phi, and psi and psi. In
other words, there is no coupling between real and imaginary parts on
the left side, but that contradicts the strong form you previously showed.

my left hand side is wrong, which in turn means that something in the creation of the weak form went wrong (else they would have been coupled). I think I found the error for that, and will fix it.

Concerning the solver: That has to be checked, too, I think GMRES will be more appropriate for that.
Thanks!