disable default decimation for nonlinear materials

NanoComp · Feb 25, 2023 · 0953062 · 0953062
1 parent ccb3c37
commit 0953062
Show file tree

Hide file tree

Showing 4 changed files with 26 additions and 2 deletions.
diff --git a/src/dft.cpp b/src/dft.cpp
@@ -204,7 +204,7 @@ dft_chunk *fields::add_dft(component c, const volume &where, const double *freq,
     double freq_max = 0;
     for (size_t i = 0; i < Nfreq; ++i)
       freq_max = std::max(freq_max, std::abs(freq[i]));
-    if ((freq_max > 0) && (src_freq_max > 0))
+    if ((freq_max > 0) && (src_freq_max > 0) && !has_nonlinearities())
       decimation_factor = std::max(1, int(std::floor(1 / (dt * (freq_max + src_freq_max)))));
     else
       decimation_factor = 1;

diff --git a/src/fields.cpp b/src/fields.cpp
@@ -663,6 +663,13 @@ void fields_chunk::use_real_fields() {
     }
 }
 
+bool fields::has_nonlinearities() const {
+  bool nonlinear = false;
+  for (int i = 0; i < num_chunks; i++)
+    if (chunks[i]->is_mine()) nonlinear = nonlinear || chunks[i]->s->has_nonlinearities();
+  return or_to_all(nonlinear);
+}
+
 int fields::phase_in_material(const structure *snew, double time) {
   if (snew->num_chunks != num_chunks)
     meep::abort("Can only phase in similar sets of chunks: %d vs %d\n", snew->num_chunks,

diff --git a/src/meep.hpp b/src/meep.hpp
@@ -105,7 +105,9 @@ class susceptibility {
   int get_id() const { return id; }
   bool operator==(const susceptibility &s) const { return id == s.id; };
 
-  // Returns the 1st order nonlinear susceptibility (generic)
+  virtual bool has_nonlinearities() const { return next ? next->has_nonlinearities() : false; }
+
+  // Returns the 1st order (linear) susceptibility (generic)
   virtual std::complex<realnum> chi1(realnum freq, realnum sigma = 1);
 
   // update all of the internal polarization state given the W field
@@ -345,6 +347,8 @@ class multilevel_susceptibility : public susceptibility {
   virtual susceptibility *clone() const { return new multilevel_susceptibility(*this); }
   virtual ~multilevel_susceptibility();
 
+  virtual bool has_nonlinearities() const { return true; }
+
   virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
                         realnum *W_prev[NUM_FIELD_COMPONENTS][2], realnum dt, const grid_volume &gv,
                         void *P_internal_data) const;
@@ -616,6 +620,8 @@ class structure_chunk {
                pml_profile_func pml_profile, void *pml_profile_data, double pml_profile_integral,
                double pml_profile_integral_u);
 
+  bool has_nonlinearities() const;
+
   void add_susceptibility(material_function &sigma, field_type ft, const susceptibility &sus);
 
   void mix_with(const structure_chunk *, double);
@@ -1753,6 +1759,7 @@ class fields {
   void reset();
   void log(const char *prefix = "");
   void change_m(double new_m);
+  bool has_nonlinearities() const;
 
   // time.cpp
   std::vector<double> time_spent_on(time_sink sink);

diff --git a/src/structure.cpp b/src/structure.cpp
@@ -545,6 +545,16 @@ bool structure_chunk::has_chi1inv(component c, direction d) const {
   return is_mine() && chi1inv[c][d] && !trivial_chi1inv[c][d];
 }
 
+bool structure_chunk::has_nonlinearities() const {
+  bool nonlinear = false;
+  if (!is_mine()) return false;
+  FOR_COMPONENTS(c) { nonlinear = nonlinear || chi2[c] || chi3[c]; }
+  FOR_FIELD_TYPES(ft) {
+    if (chiP[ft]) nonlinear = nonlinear || chiP[ft]->has_nonlinearities();
+  }
+  return nonlinear;
+}
+
 void structure::mix_with(const structure *oth, double f) {
   if (num_chunks != oth->num_chunks)
     meep::abort("You can't phase materials with different chunk topologies...\n");