<https://lh3.googleusercontent.com/-2uVGk2INmT8/WZrXBNHw4rI/AAAAAAAACRs/AK7JGees11o2KzsaUJTme-JAr8xRltYPwCLcBGAs/s1600/spectacles.png>
I have the following schrödinger equation: <https://lh3.googleusercontent.com/-qziNEIWQZ9k/WZrUPz95pOI/AAAAAAAACRM/UTjx29KCu1sSiLspuTsCqbxZWZ9_Ch_zgCEwYBhgL/s1600/schroedinger.png> which can be split up into the real and the imaginary part (as it is done in example 29): <https://lh3.googleusercontent.com/-CmPxA39KaJ4/WZrU9G-qvvI/AAAAAAAACRQ/0HeQ2QzLA_sqxNVpxZl_3sVabGYRW21FwCLcBGAs/s1600/schroedinger_split.png> This again can be solved using the RK-method by splitting the current and the future part up: <https://lh3.googleusercontent.com/-C_rfdxdr3Nc/WZrVBrIYffI/AAAAAAAACRU/G-GaIU80UvkFh3z2wf_O_nMzYMlD2Iu4ACLcBGAs/s1600/schroedinger_rk4.png> Rewriting it into the weak form results in <https://lh3.googleusercontent.com/-qQF3D1oEbFk/WZrVGxC29vI/AAAAAAAACRY/HzqDap2HxDcg__SmSVCP35YrkaGcIJ4AACLcBGAs/s1600/schroedinger_rk42.png> Thus I can create my right and my left side for the calculations in dealII as <https://lh3.googleusercontent.com/-NigH49prD-Y/WZrVUex6_AI/AAAAAAAACRc/nJtp5uw7Q-U-QwPd6WszTtYd6Cgw96BwQCLcBGAs/s1600/schroedinger_final.png> 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 <https://lh3.googleusercontent.com/-2uVGk2INmT8/WZrXBNHw4rI/AAAAAAAACRs/AK7JGees11o2KzsaUJTme-JAr8xRltYPwCLcBGAs/s1600/spectacles.png> As far as I am concerned, that looks wrong, but at which step did I make a mistake here? -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. For more options, visit https://groups.google.com/d/optout.
