Skip to content

Commit c965dd6

Browse files
roystgnrroystgnr
authored
Merge pull request #3603 from loganharbour/l2_norm_diff
Add in-place L2 norm difference method
2 parents c2728e6 + 5a0d7a2 commit c965dd6

3 files changed

Lines changed: 85 additions & 12 deletions

File tree

‎include/numerics/numeric_vector.h‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -308,6 +308,12 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
308308
*/
309309
virtual Real subset_linfty_norm (const std::set<numeric_index_type> & indices) const;
310310

311+
/**
312+
* \returns The \f$ \ell_2 \f$-norm of \f$ \vec{u} - \vec{v} \f$, where
313+
* \f$ \vec{u} \f$ is \p this.
314+
*/
315+
Real l2_norm_diff (const NumericVector<T> & other_vec) const;
316+
311317
/**
312318
* \returns The size of the vector.
313319
*/
@@ -742,6 +748,18 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
742748
*/
743749
virtual std::size_t max_allowed_id() const = 0;
744750

751+
/**
752+
* \returns whether or not this vector is able to be read for
753+
* global operations
754+
*/
755+
bool readable() const;
756+
757+
/**
758+
* \returns whether or not this vector and the vector \p v are
759+
* able to be read for global operations with one-another
760+
*/
761+
bool compatible(const NumericVector<T> & v) const;
762+
745763
protected:
746764

747765
/**

‎src/numerics/numeric_vector.C‎

Lines changed: 53 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,8 @@ template <typename T>
8484
void NumericVector<T>::insert (const T * v,
8585
const std::vector<numeric_index_type> & dof_indices)
8686
{
87+
libmesh_assert (v);
88+
8789
for (auto i : index_range(dof_indices))
8890
this->set (dof_indices[i], v[i]);
8991
}
@@ -95,6 +97,7 @@ void NumericVector<T>::insert (const NumericVector<T> & V,
9597
const std::vector<numeric_index_type> & dof_indices)
9698
{
9799
libmesh_assert_equal_to (V.size(), dof_indices.size());
100+
libmesh_assert (V.readable());
98101

99102
for (auto i : index_range(dof_indices))
100103
this->set (dof_indices[i], V(i));
@@ -106,10 +109,7 @@ template <typename T>
106109
int NumericVector<T>::compare (const NumericVector<T> & other_vector,
107110
const Real threshold) const
108111
{
109-
libmesh_assert (this->initialized());
110-
libmesh_assert (other_vector.initialized());
111-
libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
112-
libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
112+
libmesh_assert(this->compatible(other_vector));
113113

114114
int first_different_i = std::numeric_limits<int>::max();
115115
numeric_index_type i = first_local_index();
@@ -137,10 +137,7 @@ template <typename T>
137137
int NumericVector<T>::local_relative_compare (const NumericVector<T> & other_vector,
138138
const Real threshold) const
139139
{
140-
libmesh_assert (this->initialized());
141-
libmesh_assert (other_vector.initialized());
142-
libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
143-
libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
140+
libmesh_assert(this->compatible(other_vector));
144141

145142
int first_different_i = std::numeric_limits<int>::max();
146143
numeric_index_type i = first_local_index();
@@ -170,10 +167,7 @@ template <typename T>
170167
int NumericVector<T>::global_relative_compare (const NumericVector<T> & other_vector,
171168
const Real threshold) const
172169
{
173-
libmesh_assert (this->initialized());
174-
libmesh_assert (other_vector.initialized());
175-
libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
176-
libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
170+
libmesh_assert(this->compatible(other_vector));
177171

178172
int first_different_i = std::numeric_limits<int>::max();
179173
numeric_index_type i = first_local_index();
@@ -311,6 +305,8 @@ return rvalue;
311305
template <class T>
312306
Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indices) const
313307
{
308+
libmesh_assert (this->readable());
309+
314310
const NumericVector<T> & v = *this;
315311

316312
Real norm = 0;
@@ -326,6 +322,8 @@ Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indi
326322
template <class T>
327323
Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indices) const
328324
{
325+
libmesh_assert (this->readable());
326+
329327
const NumericVector<T> & v = *this;
330328

331329
Real norm = 0;
@@ -341,6 +339,8 @@ Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indi
341339
template <class T>
342340
Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> & indices) const
343341
{
342+
libmesh_assert (this->readable());
343+
344344
const NumericVector<T> & v = *this;
345345

346346
Real norm = 0;
@@ -359,10 +359,28 @@ Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> &
359359

360360

361361

362+
template <class T>
363+
Real NumericVector<T>::l2_norm_diff (const NumericVector<T> & v) const
364+
{
365+
libmesh_assert(this->compatible(v));
366+
367+
Real norm = 0;
368+
for (const auto i : make_range(this->first_local_index(), this->last_local_index()))
369+
norm += TensorTools::norm_sq((*this)(i) - v(i));
370+
371+
this->comm().sum(norm);
372+
373+
return std::sqrt(norm);
374+
}
375+
376+
377+
362378
template <typename T>
363379
void NumericVector<T>::add_vector (const T * v,
364380
const std::vector<numeric_index_type> & dof_indices)
365381
{
382+
libmesh_assert(v);
383+
366384
for (auto i : index_range(dof_indices))
367385
this->add (dof_indices[i], v[i]);
368386
}
@@ -373,6 +391,8 @@ template <typename T>
373391
void NumericVector<T>::add_vector (const NumericVector<T> & v,
374392
const std::vector<numeric_index_type> & dof_indices)
375393
{
394+
libmesh_assert(v.readable());
395+
376396
const std::size_t n = dof_indices.size();
377397
libmesh_assert_equal_to(v.size(), n);
378398
for (numeric_index_type i=0; i != n; i++)
@@ -385,11 +405,32 @@ template <typename T>
385405
void NumericVector<T>::add_vector (const NumericVector<T> & v,
386406
const ShellMatrix<T> & a)
387407
{
408+
libmesh_assert(this->compatible(v));
409+
388410
a.vector_mult_add(*this,v);
389411
}
390412

391413

392414

415+
template <typename T>
416+
bool NumericVector<T>::readable () const
417+
{
418+
return this->initialized() && this->closed();
419+
}
420+
421+
422+
template <typename T>
423+
bool NumericVector<T>::compatible (const NumericVector<T> & v) const
424+
{
425+
return this->readable() && v.readable() &&
426+
this->size() == v.size() &&
427+
this->local_size() == v.local_size() &&
428+
this->first_local_index() == v.first_local_index() &&
429+
this->last_local_index() == v.last_local_index();
430+
}
431+
432+
433+
393434
//------------------------------------------------------------------
394435
// Explicit instantiations
395436
template class LIBMESH_EXPORT NumericVector<Number>;

‎tests/numerics/numeric_vector_test.h‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -174,6 +174,20 @@ class NumericVectorTest : public CppUnit::TestCase {
174174
libMesh::TOLERANCE*libMesh::TOLERANCE);
175175
LIBMESH_ASSERT_FP_EQUAL(v.linfty_norm(), 2*(global_size-libMesh::Real(1)),
176176
libMesh::TOLERANCE*libMesh::TOLERANCE);
177+
178+
auto u_ptr = std::make_unique<Derived>(*my_comm, global_size, local_size);
179+
Base & u = *u_ptr;
180+
for (libMesh::dof_id_type n=first; n != last; n++)
181+
u.set (n, -static_cast<libMesh::Number>(n * n));
182+
u.close();
183+
184+
Derived diff_derived(*my_comm, global_size, local_size);
185+
Base & diff = diff_derived;
186+
diff = u;
187+
diff -= v;
188+
diff.close();
189+
LIBMESH_ASSERT_FP_EQUAL(diff.l2_norm(), u.l2_norm_diff(v),
190+
libMesh::TOLERANCE*libMesh::TOLERANCE);
177191
}
178192

179193
template <class Base, class Derived>

0 commit comments

Comments
 (0)