Skip to content

Clear up duplicate Newton Stokes assemblers #6860

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
YiminJin wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Clear up duplicate Newton Stokes assemblers #6860

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 0 additions & 68 deletions include/aspect/newton.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -307,74 +307,6 @@ namespace aspect
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
};

/**
* This class assembles the right-hand-side term of the Newton Stokes system
* that is caused by the compressibility in the mass conservation equation.
* This function approximates this term as
* $- \nabla \mathbf{u} = \frac{1}{\rho} \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u}$
*/
template <int dim>
class NewtonStokesReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
public:
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
};

/**
* This class assembles the compressibility term of the Newton Stokes system
* that is caused by the compressibility in the mass conservation equation.
* It includes this term implicitly in the matrix,
* which is therefore not longer symmetric.
* This function approximates this term as
* $ - \nabla \mathbf{u} - \frac{1}{\rho} \frac{\partial rho}{\partial z} \frac{\mathbf{g}}{||\mathbf{g}||} \cdot \mathbf{u} = 0$
*/
template <int dim>
class NewtonStokesImplicitReferenceDensityCompressibilityTerm : public Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
public:
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
};

/**
* This class assembles the right-hand-side term of the Newton Stokes system
* that is caused by the compressibility in the mass conservation equation.
* This function approximates this term as
* $ - \nabla \mathbf{u} = \frac{1}{\rho} \frac{\partial rho}{\partial p} \rho \mathbf{g} \cdot \mathbf{u}$
*/
template <int dim>
class NewtonStokesIsentropicCompressionTerm : public Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
public:
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const override;
};

/**
* This class assembles the right-hand-side term of the Stokes equation
* that is caused by the variable density in the mass conservation equation.
* This class approximates this term as
* $ - \nabla \cdot \mathbf{u} = \frac{1}{\rho} \frac{\partial \rho}{\partial t} + \frac{1}{\rho} \nabla \rho \cdot \mathbf{u}$
* where the right-hand side velocity is explicitly taken from the last timestep,
* and the density is taken from a compositional field of the type 'density'.
*/
template <int dim>
class NewtonStokesProjectedDensityFieldTerm : public Assemblers::Interface<dim>,
public SimulatorAccess<dim>
{
public:
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
internal::Assembly::CopyData::CopyDataBase<dim> &data) const override;
};
}
}

Expand Down
259 changes: 1 addition & 258 deletions source/simulator/assemblers/newton_stokes.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -992,259 +992,6 @@ namespace aspect
}
}
}



template <int dim>
void
NewtonStokesReferenceDensityCompressibilityTerm<dim>::
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
{
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);

// assemble RHS of:
// - div u = 1/rho * drho/dz g/||g||* u
Assert(this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::reference_density_profile,
ExcInternalError());

const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
const double pressure_scaling = this->get_pressure_scaling();

for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
++i_stokes;
}
++i;
}

const Tensor<1,dim>
gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
const double drho_dz_u = scratch.reference_densities_depth_derivative[q]
* (gravity * scratch.velocity_values[q]) / gravity.norm();
const double one_over_rho = 1.0/scratch.reference_densities[q];
const double JxW = scratch.finite_element_values.JxW(q);

for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
data.local_rhs(i) += (pressure_scaling *
one_over_rho * drho_dz_u * scratch.phi_p[i])
* JxW;
}
}



template <int dim>
void
NewtonStokesImplicitReferenceDensityCompressibilityTerm<dim>::
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
{
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);

// assemble compressibility term of:
// - div u - 1/rho * drho/dz g/||g||* u = 0
Assert(this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile,
ExcInternalError());

if (!scratch.rebuild_stokes_matrix)
return;

const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
const double pressure_scaling = this->get_pressure_scaling();

for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.phi_u[i_stokes] = scratch.finite_element_values[introspection.extractors.velocities].value (i,q);
scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i,q);
++i_stokes;
}
++i;
}

const Tensor<1,dim>
gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));
const Tensor<1,dim> drho_dz = scratch.reference_densities_depth_derivative[q]
* gravity / gravity.norm();
const double one_over_rho = 1.0/scratch.reference_densities[q];
const double JxW = scratch.finite_element_values.JxW(q);

for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
data.local_matrix(i,j) += (pressure_scaling *
one_over_rho * drho_dz * scratch.phi_u[j] * scratch.phi_p[i])
* JxW;
}
}



template <int dim>
void
NewtonStokesIsentropicCompressionTerm<dim>::
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
{
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);

// assemble RHS of:
// - div u = 1/rho * drho/dp rho * g * u
Assert(this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::isentropic_compression,
ExcInternalError());

const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
const double pressure_scaling = this->get_pressure_scaling();

for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
++i_stokes;
}
++i;
}

const Tensor<1,dim>
gravity = this->get_gravity_model().gravity_vector (scratch.finite_element_values.quadrature_point(q));

const double compressibility
= scratch.material_model_outputs.compressibilities[q];

const double density = scratch.material_model_outputs.densities[q];
const double JxW = scratch.finite_element_values.JxW(q);

for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
data.local_rhs(i) += (
// add the term that results from the compressibility. compared
// to the manual, this term seems to have the wrong sign, but this
// is because we negate the entire equation to make sure we get
// -div(u) as the adjoint operator of grad(p)
(pressure_scaling *
compressibility * density *
(scratch.velocity_values[q] * gravity) *
scratch.phi_p[i])
)
* JxW;
}
}



template <int dim>
void
NewtonStokesProjectedDensityFieldTerm<dim>::
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
{
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>&> (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>&> (data_base);

// assemble RHS of:
// $ - \nabla \cdot \mathbf{u} = \frac{1}{\rho} \frac{\partial \rho}{\partial t} + \frac{1}{\rho} \nabla \rho \cdot \mathbf{u}$

// Compared to the manual, this term seems to have the wrong sign, but
// this is because we negate the entire equation to make sure we get
// -div(u) as the adjoint operator of grad(p)

Assert(this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::projected_density_field,
ExcInternalError());

const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
const double pressure_scaling = this->get_pressure_scaling();
const unsigned int density_idx = this->introspection().find_composition_type(CompositionalFieldDescription::density);

const double time_step = this->get_timestep();
const double old_time_step = this->get_old_timestep();

std::vector<Tensor<1,dim>> density_gradients(n_q_points);
std::vector<double> density(n_q_points);
std::vector<double> density_old(n_q_points);
std::vector<double> density_old_old(n_q_points);

scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_gradients (this->get_current_linearization_point(),
density_gradients);
scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_values (this->get_current_linearization_point(),
density);
scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_values (this->get_old_solution(),
density_old);
scratch.finite_element_values[introspection.extractors.compositional_fields[density_idx]].get_function_values (this->get_old_old_solution(),
density_old_old);

for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int i=0, i_stokes=0; i_stokes<stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch.phi_p[i_stokes] = scratch.finite_element_values[introspection.extractors.pressure].value (i, q);
++i_stokes;
}
++i;
}

double drho_dt;

if (this->get_timestep_number() > 1)
drho_dt = (1.0/time_step) *
(density[q] *
(2*time_step + old_time_step) / (time_step + old_time_step)
-
density_old[q] *
(1 + time_step/old_time_step)
+
density_old_old[q] *
(time_step * time_step) / (old_time_step * (time_step + old_time_step)));
else if (this->get_timestep_number() == 1)
drho_dt =
(density[q] - density_old[q]) / time_step;
else
drho_dt = 0.0;

const double JxW = scratch.finite_element_values.JxW(q);

for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
data.local_rhs(i) += (
(pressure_scaling *
(1.0 / density[q]) *
(density_gradients[q] *
scratch.velocity_values[q]) *
scratch.phi_p[i])
+ pressure_scaling * (1.0 / density[q]) * drho_dt * scratch.phi_p[i]
)
* JxW;
}
}
}
} // namespace aspect

Expand All @@ -1258,11 +1005,7 @@ namespace aspect
template class NewtonStokesPreconditioner<dim>; \
template class NewtonStokesCompressiblePreconditioner<dim>; \
template class NewtonStokesIncompressibleTerms<dim>; \
template class NewtonStokesCompressibleStrainRateViscosityTerm<dim>; \
template class NewtonStokesReferenceDensityCompressibilityTerm<dim>; \
template class NewtonStokesImplicitReferenceDensityCompressibilityTerm<dim>; \
template class NewtonStokesIsentropicCompressionTerm<dim>; \
template class NewtonStokesProjectedDensityFieldTerm<dim>;
template class NewtonStokesCompressibleStrainRateViscosityTerm<dim>;

ASPECT_INSTANTIATE(INSTANTIATE)

Expand Down
8 changes: 4 additions & 4 deletions source/simulator/newton.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,13 +67,13 @@ namespace aspect
Parameters<dim>::Formulation::MassConservation::implicit_reference_density_profile)
{
assemblers.stokes_system.push_back(
std::make_unique<aspect::Assemblers::NewtonStokesImplicitReferenceDensityCompressibilityTerm<dim>>());
std::make_unique<aspect::Assemblers::StokesImplicitReferenceDensityCompressibilityTerm<dim>>());
}
else if (this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::reference_density_profile)
{
assemblers.stokes_system.push_back(
std::make_unique<aspect::Assemblers::NewtonStokesReferenceDensityCompressibilityTerm<dim>>());
std::make_unique<aspect::Assemblers::StokesReferenceDensityCompressibilityTerm<dim>>());
}
else if (this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::incompressible)
Expand All @@ -84,14 +84,14 @@ namespace aspect
Parameters<dim>::Formulation::MassConservation::isentropic_compression)
{
assemblers.stokes_system.push_back(
std::make_unique<aspect::Assemblers::NewtonStokesIsentropicCompressionTerm<dim>>());
std::make_unique<aspect::Assemblers::StokesIsentropicCompressionTerm<dim>>());
}
else if (this->get_parameters().formulation_mass_conservation ==
Parameters<dim>::Formulation::MassConservation::projected_density_field)
{
CitationInfo::add("pda");
assemblers.stokes_system.push_back(
std::make_unique<aspect::Assemblers::NewtonStokesProjectedDensityFieldTerm<dim>>());
std::make_unique<aspect::Assemblers::StokesProjectedDensityFieldTerm<dim>>());
}
else
AssertThrow(false,
Expand Down