Skip to content

Commit

Permalink
Unify assembly patterns
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller
gassmoeller committed Nov 22, 2024
1 parent 7f93e8a commit 6a41f4e
Show file tree
Hide file tree
Showing 9 changed files with 163 additions and 107 deletions.
21 changes: 11 additions & 10 deletions source/mesh_deformation/diffusion.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,10 +281,10 @@ namespace aspect
cell_matrix = 0;

// Loop over the quadrature points of the current face
for (unsigned int point=0; point<n_fs_face_q_points; ++point)
for (unsigned int q=0; q<n_fs_face_q_points; ++q)
{
// Get the gravity vector to compute the outward direction of displacement
Tensor<1,dim> direction = -(this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(point)));
Tensor<1,dim> direction = -(this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(q)));
// Normalize direction vector
if (direction.norm() > 0.0)
direction *= 1./direction.norm();
Expand All @@ -299,15 +299,16 @@ namespace aspect

// Compute the total displacement in the gravity direction,
// i.e. the initial topography + any additional mesh displacement.
const double displacement = direction * (displacement_values[point] + initial_topography_values[point]);
const double displacement = direction * (displacement_values[q] + initial_topography_values[q]);

// To project onto the tangent space of the surface,
// we define the projection P:= I- n x n,
// with I the unit tensor and n the unit normal to the surface.
// The surface gradient then is P times the usual gradient of the shape functions.
const Tensor<2, dim, double> projection = unit_symmetric_tensor<dim>() -
outer_product(fs_fe_face_values.normal_vector(point), fs_fe_face_values.normal_vector(point));
outer_product(fs_fe_face_values.normal_vector(q), fs_fe_face_values.normal_vector(q));

const double JxW = fs_fe_face_values.JxW(q);

// The shape values for the i-loop
std::vector<double> phi(dofs_per_cell);
Expand All @@ -327,12 +328,12 @@ namespace aspect
if (mesh_deformation_dof_handler.get_fe().system_to_component_index(i).first == dim-1)
{
// Precompute shape values and projected shape value gradients
phi[i] = fs_fe_face_values.shape_value (i, point);
projected_grad_phi[i] = projection * fs_fe_face_values.shape_grad(i, point);
phi[i] = fs_fe_face_values.shape_value (i, q);
projected_grad_phi[i] = projection * fs_fe_face_values.shape_grad(i, q);

// Assemble the RHS
// RHS = M*H_old
cell_vector(i) += phi[i] * displacement * fs_fe_face_values.JxW(point);
cell_vector(i) += phi[i] * displacement * JxW;

for (unsigned int j=0; j<dofs_per_cell; ++j)
{
Expand All @@ -344,12 +345,12 @@ namespace aspect
// Matrix := (M+dt*K) = (M+dt*B^T*kappa*B)
cell_matrix(i,j) +=
(
phi[i] * fs_fe_face_values.shape_value (j, point) +
phi[i] * fs_fe_face_values.shape_value (j, q) +
this->get_timestep() * diffusivity *
projected_grad_phi[i] *
(projection * fs_fe_face_values.shape_grad(j, point))
(projection * fs_fe_face_values.shape_grad(j, q))
)
* fs_fe_face_values.JxW(point);
* JxW;
}
}
}
Expand Down
24 changes: 15 additions & 9 deletions source/mesh_deformation/free_surface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,31 +156,37 @@ namespace aspect

cell_vector = 0;
cell_matrix = 0;
for (unsigned int point=0; point<n_face_q_points; ++point)
for (unsigned int q=0; q<n_face_q_points; ++q)
{
// Select the direction onto which to project the velocity solution
Tensor<1,dim> direction;
if ( advection_direction == SurfaceAdvection::normal ) // project onto normal vector
direction = fs_fe_face_values.normal_vector(point);
direction = fs_fe_face_values.normal_vector(q);
else if ( advection_direction == SurfaceAdvection::vertical ) // project onto local gravity
direction = this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(point));
direction = this->get_gravity_model().gravity_vector(fs_fe_face_values.quadrature_point(q));
else
AssertThrow(false, ExcInternalError());

direction *= ( direction.norm() > 0.0 ? 1./direction.norm() : 0.0 );

const double JxW = fs_fe_face_values.JxW(q);

small_vector<Tensor<1,dim>> phi_u(dofs_per_cell);
for (unsigned int i=0; i<dofs_per_cell; ++i)
phi_u[i] = fs_fe_face_values[extract_vel].value(i,q);

for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix(i,j) += (fs_fe_face_values[extract_vel].value(j,point) *
fs_fe_face_values[extract_vel].value(i,point) ) *
fs_fe_face_values.JxW(point);
cell_matrix(i,j) += (phi_u[j] *
phi_u[i]) *
JxW;
}

cell_vector(i) += (fs_fe_face_values[extract_vel].value(i,point) * direction)
* (velocity_values[point] * direction)
* fs_fe_face_values.JxW(point);
cell_vector(i) += (phi_u[i] * direction)
* (velocity_values[q] * direction)
* JxW;
}
}

Expand Down
43 changes: 28 additions & 15 deletions source/mesh_deformation/interface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,39 +109,49 @@ namespace aspect

this->get_material_model().evaluate(scratch.face_material_model_inputs, scratch.face_material_model_outputs);

for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
for (unsigned int q = 0; q < n_face_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.face_finite_element_values[introspection.extractors.velocities].value(i, q_point);
scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection.extractors.velocities].value(i, q);
++i_stokes;
}
++i;
}

const Tensor<1,dim>
gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q_point));
gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q));
const double g_norm = gravity.norm();

// construct the relevant vectors
const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point);
const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q);
const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm);

const double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] *
small_vector<double> phi_u_times_g_hat(stokes_dofs_per_cell);
small_vector<double> phi_u_times_n_hat(stokes_dofs_per_cell);
for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
{
phi_u_times_g_hat[i] = scratch.phi_u[i] * g_hat;
phi_u_times_n_hat[i] = scratch.phi_u[i] * n_hat;
}

const double pressure_perturbation = scratch.face_material_model_outputs.densities[q] *
this->get_timestep() *
free_surface_theta *
g_norm;

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

// see Kaus et al 2010 for details of the stabilization term
for (unsigned int i=0; i< stokes_dofs_per_cell; ++i)
for (unsigned int j=0; j< stokes_dofs_per_cell; ++j)
{
// The fictive stabilization stress is (phi_u[i].g)*(phi_u[j].n)
const double stress_value = -pressure_perturbation*
(scratch.phi_u[i]*g_hat) * (scratch.phi_u[j]*n_hat)
*scratch.face_finite_element_values.JxW(q_point);
phi_u_times_g_hat[i] * phi_u_times_n_hat[j]
* JxW;

data.local_matrix(i,j) += stress_value;
}
Expand Down Expand Up @@ -867,14 +877,17 @@ namespace aspect

cell_vector = 0;
cell_matrix = 0;
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j) += scalar_product( fe_values[extract_vel].gradient(i,point),
fe_values[extract_vel].gradient(j,point) ) *
fe_values.JxW(point);
}
for (unsigned int q=0; q<n_q_points; ++q)
{
const double JxW = fe_values.JxW(q);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
for (unsigned int j=0; j<dofs_per_cell; ++j)
cell_matrix(i,j) += scalar_product( fe_values[extract_vel].gradient(i,q),
fe_values[extract_vel].gradient(j,q) ) *
JxW;
}
}

mesh_velocity_constraints.distribute_local_to_global (cell_matrix, cell_vector,
cell_dof_indices, mesh_matrix, rhs, false);
Expand Down
28 changes: 18 additions & 10 deletions source/postprocess/dynamic_topography.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,7 @@ namespace aspect
const double density = out_volume.densities[q];
const bool is_compressible = this->get_material_model().is_compressible();
const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(in_volume.position[q]);
const double JxW = fe_volume_values.JxW(q);

// Set up shape function values
for (unsigned int k=0; k<dofs_per_cell; ++k)
Expand All @@ -173,20 +174,26 @@ namespace aspect
{
// Viscous stress part
local_vector(i) += 2.0 * eta * ( epsilon_phi_u[i] * in_volume.strain_rate[q]
- (is_compressible ? 1./3. * div_phi_u[i] * div_solution[q] : 0.0) ) * fe_volume_values.JxW(q);
- (is_compressible ? 1./3. * div_phi_u[i] * div_solution[q] : 0.0) ) * JxW;
// Pressure and compressibility parts
local_vector(i) -= div_phi_u[i] * in_volume.pressure[q] * fe_volume_values.JxW(q);
local_vector(i) -= div_phi_u[i] * in_volume.pressure[q] * JxW;
// Force part
local_vector(i) -= density * gravity * phi_u[i] * fe_volume_values.JxW(q);
local_vector(i) -= density * gravity * phi_u[i] * JxW;
}
}
// Assemble the mass matrix for cell face. Since we are using GLL
// quadrature, the mass matrix will be diagonal, and we can just assemble it into a vector.
for (unsigned int q=0; q < n_face_q_points; ++q)
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_mass_matrix(i) += fe_face_values[this->introspection().extractors.velocities].value(i,q) *
fe_face_values[this->introspection().extractors.velocities].value(i,q) *
fe_face_values.JxW(q);
{
const double JxW = fe_face_values.JxW(q);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const Tensor<1,dim> phi_u = fe_face_values[this->introspection().extractors.velocities].value(i,q);
local_mass_matrix(i) += phi_u *
phi_u *
JxW;
}
}

cell->distribute_local_to_global(local_vector, rhs_vector);
cell->distribute_local_to_global(local_mass_matrix, mass_matrix);
Expand Down Expand Up @@ -330,21 +337,22 @@ namespace aspect
const Point<dim> point = fe_output_values.quadrature_point(q);
const Tensor<1,dim> normal = fe_output_values.normal_vector(q);
const double gravity_norm = this->get_gravity_model().gravity_vector(point).norm();
const double JxW = fe_output_values.JxW(q);

if (at_upper_surface)
{
const double delta_rho = out_output.densities[q] - density_above;
dynamic_topography += (-stress_output_values[q]*normal - surface_pressure)
/ delta_rho / gravity_norm * fe_output_values.JxW(q);
/ delta_rho / gravity_norm * JxW;
}
else
{
const double delta_rho = out_output.densities[q] - density_below;

dynamic_topography += (-stress_output_values[q]*normal - bottom_pressure)
/ delta_rho / gravity_norm * fe_output_values.JxW(q);
/ delta_rho / gravity_norm * JxW;
}
face_area += fe_output_values.JxW(q);
face_area += JxW;
}
// Get the average dynamic topography for the cell
dynamic_topography = dynamic_topography * (backward_advection ? -1. : 1.) / face_area;
Expand Down
28 changes: 16 additions & 12 deletions source/simulator/assemblers/stokes.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -908,16 +908,19 @@ namespace aspect
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;

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);
data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * scratch.finite_element_values.JxW(q);
++i_stokes;
}
++i;
}
{
const double JxW = scratch.finite_element_values.JxW(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);
data.local_pressure_shape_function_integrals(i_stokes) += scratch.phi_p[i_stokes] * JxW;
++i_stokes;
}
++i;
}
}
}


Expand Down Expand Up @@ -951,13 +954,14 @@ namespace aspect
scratch.face_finite_element_values.quadrature_point(q),
scratch.face_finite_element_values.normal_vector(q));

const double JxW = scratch.face_finite_element_values.JxW(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))
{
data.local_rhs(i_stokes) += scratch.face_finite_element_values[introspection.extractors.velocities].value(i,q) *
traction *
scratch.face_finite_element_values.JxW(q);
traction * JxW;
++i_stokes;
}
++i;
Expand Down
Loading

0 comments on commit 6a41f4e

Please sign in to comment.