Thank you Wolfgang.
As I changing built-in solver to PETScWrappers, I failed to converse a
"PETSc::Wrappers:Vector" variable to a "Vector<std::complex<double>>" one,
as shown in line 464 of my code. I tried
"Vector<std::complex<double>> a = b; (where b is a PETSc Vector)"
or
"Vector<std::complex<double>> a(b)"
but both failed. In step-17 however, both is good. Is it because the
complex type or other mistake inmy code?
Best
Wayne
在2022年9月26日星期一 UTC+8 23:32:40<Wolfgang Bangerth> 写道：
> On 9/26/22 09:22, 'yy.wayne' via deal.II User Group wrote:
> > I'm trying to apply a Krylov solver for complex-valued problem. Basiclly
> > I have modified step-16 into a complex-valued problem (do not seperate
> > real and imaginary parts) and solve with direct solver successfully.
> > However it failed on GMRES solver. Does deal.II built-in GMRES solver
> > only supprts double type?
>
> Yes. You probably want to take a look at step-58 as well, along with the
> discussion in the results section there.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
--
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].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/ffc4cafc-73b6-4b61-97e5-cc82b2be91c3n%40googlegroups.com.
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <iostream>
#include <fstream>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/base/mpi.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/base/timer.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/numerics/error_estimator.h>
// complex valued version of Step29
// adaptive mesh
namespace Step29
{
using namespace dealii;
template <int dim>
class DirichletBoundaryValues : public Function<dim, std::complex<double>>
{
public:
DirichletBoundaryValues()
: Function<dim, std::complex<double>>(1)
{}
virtual std::complex<double> value(const Point<dim> &p, const unsigned int component = 0) const override;
};
template <int dim>
std::complex<double> DirichletBoundaryValues<dim>::value(const Point<dim> &p, const unsigned int component) const
{
(void)component;
Assert(component==0, ExcIndexRange(component,0,1));
return std::complex<double>(1,0);
}
class ParameterReader : public Subscriptor
{
public:
ParameterReader(ParameterHandler &);
void read_parameters(const std::string &);
private:
void declare_parameters();
ParameterHandler &prm;
};
ParameterReader::ParameterReader(ParameterHandler ¶mhandler)
: prm(paramhandler)
{}
void ParameterReader::declare_parameters()
{
prm.enter_subsection("Mesh & geometry parameters");
{
prm.declare_entry("Number of refinements",
"6",
Patterns::Integer(0),
"Number of global mesh refinement steps "
"applied to initial coarse grid");
prm.declare_entry("Focal distance",
"0.3",
Patterns::Double(0),
"Distance of the focal point of the lens "
"to the x-axis");
}
prm.leave_subsection();
prm.enter_subsection("Physical constants");
{
prm.declare_entry("c", "1.5e5", Patterns::Double(0), "Wave speed");
prm.declare_entry("omega", "5.0e7", Patterns::Double(0), "Frequency");
}
prm.leave_subsection();
prm.enter_subsection("Output parameters");
{
prm.declare_entry("Output filename",
"solution",
Patterns::Anything(),
"Name of the output file (without extension)");
DataOutInterface<1>::declare_parameters(prm);
}
prm.leave_subsection();
}
void ParameterReader::read_parameters(const std::string ¶meter_file)
{
declare_parameters();
prm.parse_input(parameter_file);
}
template <int dim>
class ComputeIntensity : public DataPostprocessorScalar<dim>
{
public:
ComputeIntensity();
virtual void evaluate_vector_field(
const DataPostprocessorInputs::Vector<dim> &inputs,
std::vector<Vector<double>> &computed_quantities) const override;
};
template <int dim>
ComputeIntensity<dim>::ComputeIntensity()
: DataPostprocessorScalar<dim>("Intensity", update_values)
{}
template <int dim>
void ComputeIntensity<dim>::evaluate_vector_field(
const DataPostprocessorInputs::Vector<dim> &inputs,
std::vector<Vector<double>> & computed_quantities) const
{
AssertDimension(computed_quantities.size(), inputs.solution_values.size());
for (unsigned int i = 0; i < computed_quantities.size(); ++i)
{
AssertDimension(computed_quantities[i].size(), 1);
AssertDimension(inputs.solution_values[i].size(), 2);
const std::complex<double> u(inputs.solution_values[i](0),
inputs.solution_values[i](1));
computed_quantities[i](0) = std::abs(u);
}
}
template <int dim>
class UltrasoundProblem
{
public:
UltrasoundProblem(ParameterHandler &);
void run();
private:
// void make_grid();
void refine_grid();
void setup_system();
void assemble_system();
void solve();
void output_results(const unsigned int cycle) const;
ParameterHandler &prm;
Triangulation<dim> triangulation;
DoFHandler<dim> dof_handler;
FE_Q<dim> fe;
AffineConstraints<std::complex<double>> constraints;
// SparsityPattern sparsity_pattern;
// SparseMatrix<std::complex<double>> system_matrix;
// Vector<std::complex<double>> solution, system_rhs;
MPI_Comm mpi_communicator;
const unsigned int n_mpi_processes;
const unsigned int this_mpi_process;
PETScWrappers::MPI::SparseMatrix system_matrix;
PETScWrappers::MPI::Vector solution, system_rhs;
};
template <int dim>
UltrasoundProblem<dim>::UltrasoundProblem(ParameterHandler ¶m)
: prm(param)
, dof_handler(triangulation)
, fe(1)
, mpi_communicator(MPI_COMM_WORLD)
, n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
, this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
{}
#if 0
template <int dim>
void UltrasoundProblem<dim>::make_grid()
{
std::cout << "Generating grid... ";
Timer timer;
prm.enter_subsection("Mesh & geometry parameters");
const double focal_distance = prm.get_double("Focal distance");
const unsigned int n_refinements = prm.get_integer("Number of refinements");
prm.leave_subsection();
const Point<dim> transducer =
(dim == 2) ? Point<dim>(0.5, 0.0) : Point<dim>(0.5, 0.5, 0.0);
const Point<dim> focal_point = (dim == 2) ?
Point<dim>(0.5, focal_distance) :
Point<dim>(0.5, 0.5, focal_distance);
GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1);
for (auto &cell : triangulation.cell_iterators())
for (const auto &face : cell->face_iterators())
if (face->at_boundary() &&
((face->center() - transducer).norm_square() < 0.01))
{
face->set_boundary_id(1);
face->set_manifold_id(1);
}
triangulation.set_manifold(1, SphericalManifold<dim>(focal_point));
triangulation.refine_global(n_refinements);
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
#endif
template <int dim>
void UltrasoundProblem<dim>::refine_grid()
{
Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
KellyErrorEstimator<dim>::estimate(dof_handler,
QGauss<dim - 1>(fe.degree + 1),
{},
solution,
estimated_error_per_cell);
GridRefinement::refine_and_coarsen_fixed_number(triangulation,
estimated_error_per_cell,
0.3,
0.03);
triangulation.execute_coarsening_and_refinement();
}
template <int dim>
void UltrasoundProblem<dim>::setup_system()
{
std::cout << "Setting up system... ";
Timer timer;
GridTools::partition_triangulation(n_mpi_processes, triangulation);
dof_handler.distribute_dofs(fe);
DoFRenumbering::subdomain_wise(dof_handler);
constraints.clear();
DoFTools::make_hanging_node_constraints(dof_handler, constraints);
VectorTools::interpolate_boundary_values(dof_handler,
1,
DirichletBoundaryValues<dim>(),
constraints);
constraints.close();
DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, /*keep_constrained_dofs = */ false);
const std::vector<IndexSet> locally_owned_dofs_per_proc =
DoFTools::locally_owned_dofs_per_subdomain(dof_handler);
const IndexSet locally_owned_dofs =
locally_owned_dofs_per_proc[this_mpi_process];
system_matrix.reinit(locally_owned_dofs,
locally_owned_dofs,
dsp,
mpi_communicator);
solution.reinit(locally_owned_dofs, mpi_communicator);
system_rhs.reinit(locally_owned_dofs, mpi_communicator);
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
<< std::endl;
}
template <int dim>
void UltrasoundProblem<dim>::assemble_system()
{
std::cout << "Assembling system matrix... ";
Timer timer;
prm.enter_subsection("Physical constants");
const double omega = prm.get_double("omega"), c = prm.get_double("c");
prm.leave_subsection();
QGauss<dim> quadrature_formula(fe.degree + 1);
QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
const unsigned int n_q_points = quadrature_formula.size(),
n_face_q_points = face_quadrature_formula.size(),
dofs_per_cell = fe.n_dofs_per_cell();
std::cout << "fe.degree:" << fe.degree << std::endl;
std::cout << "n_dofs_per_cell:" << dofs_per_cell << std::endl;
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values | update_gradients |
update_JxW_values);
FEFaceValues<dim> fe_face_values(fe,
face_quadrature_formula,
update_values | update_JxW_values);
FullMatrix<std::complex<double>> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<std::complex<double>> cell_rhs(dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
for (const auto &cell : dof_handler.active_cell_iterators())
{
if (cell->subdomain_id() == this_mpi_process)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
for (unsigned int q_point = 0; q_point < n_q_points;
++q_point)
cell_matrix(i, j) +=
(((fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point)) *
(-omega * omega) +
(fe_values.shape_grad(i, q_point) *
fe_values.shape_grad(j, q_point)) *
c * c) *
fe_values.JxW(q_point));
}
cell_rhs(i) = 0;
}
for (const auto face_no : cell->face_indices())
if (cell->face(face_no)->at_boundary() &&
(cell->face(face_no)->boundary_id() == 0))
{
fe_face_values.reinit(cell, face_no);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
for (unsigned int q_point = 0; q_point < n_face_q_points;
++q_point)
cell_matrix(i, j) +=
fe_face_values.shape_value(i, q_point) *
fe_face_values.shape_value(j, q_point) * std::complex<double>(0,1) * c * omega *
fe_face_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(
cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
}
}
system_matrix.compress(VectorOperation::add);
system_rhs.compress(VectorOperation::add);
// std::map<types::global_dof_index, std::complex<double>> boundary_values;
// VectorTools::interpolate_boundary_values(dof_handler,
// 1,
// DirichletBoundaryValues<dim>(),
// boundary_values);
// VectorTools::interpolate_boundary_values(dof_handler,
// 1,
// Functions::ZeroFunction<dim, std::complex<double>>(),
// boundary_values);
// MatrixTools::apply_boundary_values(boundary_values,
// system_matrix,
// solution,
// system_rhs);
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
}
template <int dim>
void UltrasoundProblem<dim>::solve()
{
std::cout << "Solving linear system... ";
Timer timer;
#if 0
SparseDirectUMFPACK A_direct;
// A_direct.initialize(system_matrix);
// A_direct.vmult(solution, system_rhs);
A_direct.initialize(system_matrix);
A_direct.solve(system_matrix, system_rhs);
solution = system_rhs;
constraints.distribute(solution);
#else
//type of solution: PETScWrappers::MPI::Vector
SolverControl solver_control(1000,1e-4);
PETScWrappers::SolverGMRES gmres(solver_control);
PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
gmres.solve(system_matrix, solution, system_rhs, preconditioner);
Vector<std::complex<double>> localized_solution = solution;
constraints.distribute(localized_solution);
solution = localized_solution;
// constraints.distribute(solution);
#endif
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
std::cout << " Number of GMRES iterations: " << solver_control.last_step()
<< std::endl;
}
template <int dim>
void UltrasoundProblem<dim>::output_results(const unsigned int cycle) const
{
std::cout << "Generating output... ";
Timer timer;
ComputeIntensity<dim> intensities; //const DataPostprocessors::ComputeIntensity<dim> intensities;
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
prm.enter_subsection("Output parameters");
const std::string output_filename = prm.get("Output filename");
data_out.parse_parameters(prm);
prm.leave_subsection();
const std::string filename = output_filename + "_" + std::to_string(cycle) + data_out.default_suffix();
std::ofstream output(filename);
// std::vector<std::string> solution_names;
// solution_names.emplace_back("Re_u");
// solution_names.emplace_back("Im_u");
data_out.add_data_vector(solution, "u");
std::cout << "write solution u complete" << std::endl;
data_out.add_data_vector(solution, intensities);
data_out.build_patches();
data_out.write(output);
timer.stop();
std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
}
template <int dim>
void UltrasoundProblem<dim>::run()
{
std::cout << "Generating grid... ";
prm.enter_subsection("Mesh & geometry parameters");
const double focal_distance = prm.get_double("Focal distance");
const unsigned int n_refinements = prm.get_integer("Number of refinements");
prm.leave_subsection();
const Point<dim> transducer =
(dim == 2) ? Point<dim>(0.5, 0.0) : Point<dim>(0.5, 0.5, 0.0);
const Point<dim> focal_point = (dim == 2) ?
Point<dim>(0.5, focal_distance) :
Point<dim>(0.5, 0.5, focal_distance);
GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1);
for (auto &cell : triangulation.cell_iterators())
for (const auto &face : cell->face_iterators())
if (face->at_boundary() &&
((face->center() - transducer).norm_square() < 0.01))
{
face->set_boundary_id(1);
face->set_manifold_id(1);
}
triangulation.set_manifold(1, SphericalManifold<dim>(focal_point));
// triangulation.refine_global(1);
for (unsigned int cycle = 0; cycle < n_refinements + 1; ++ cycle)
{
if (cycle==0)
triangulation.refine_global(3);
else
refine_grid();
std::cout << " cycle: " << cycle << std::endl;
std::cout << " Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
setup_system();
assemble_system();
solve();
output_results(cycle);
}
}
} // namespace Step29
int main()
{
try
{
using namespace dealii;
using namespace Step29;
ParameterHandler prm;
ParameterReader param(prm);
param.read_parameters("step-29.prm");
UltrasoundProblem<2> ultrasound_problem(prm);
ultrasound_problem.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
