← Back to branch summary

~kroq-gar78/ubuntu/precise/meep/fix-990137

« back to all changes in this revision

Viewing changes to src/update_eh.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Thorsten Alteholz
  • Date: 2009-10-29 18:00:00 UTC
  • mfrom: (2.2.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20091029180000-wekpivz2ct4sggty
Tags: 1.1.1-3
http://bugs.debian.org/551822
* debian/rules: change configure script to produce meep-mpi paths
                everywhere (closes: #551822)
* debian/control: DM-Upload-Allowed added

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/update_eh.cpp
 
1
/* Copyright (C) 2005-2009 Massachusetts Institute of Technology
 
2
%
 
3
%  This program is free software; you can redistribute it and/or modify
 
4
%  it under the terms of the GNU General Public License as published by
 
5
%  the Free Software Foundation; either version 2, or (at your option)
 
6
%  any later version.
 
7
%
 
8
%  This program is distributed in the hope that it will be useful,
 
9
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
%  GNU General Public License for more details.
 
12
%
 
13
%  You should have received a copy of the GNU General Public License
 
14
%  along with this program; if not, write to the Free Software Foundation,
 
15
%  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 
16
*/
 
17
 
 
18
#include <string.h>
 
19
 
 
20
#include "meep.hpp"
 
21
#include "meep_internals.hpp"
 
22
 
 
23
namespace meep {
 
24
  
 
25
void fields::update_eh(field_type ft, bool skip_w_components) {
 
26
  if (ft != E_stuff && ft != H_stuff) abort("update_eh only works with E/H");
 
27
  for (int i=0;i<num_chunks;i++)
 
28
    if (chunks[i]->is_mine())
 
29
      if (chunks[i]->update_eh(ft, skip_w_components))
 
30
        chunk_connections_valid = false; // E/H allocated - reconnect chunks 
 
31
 
 
32
  /* synchronize to avoid deadlocks if one process decides it needs
 
33
     to allocate E or H ... */
 
34
  chunk_connections_valid = and_to_all(chunk_connections_valid);
 
35
}
 
36
 
 
37
bool fields_chunk::update_eh(field_type ft, bool skip_w_components) {
 
38
  field_type ft2 = ft == E_stuff ? D_stuff : B_stuff; // for sources etc.
 
39
  bool allocated_eh = false;
 
40
 
 
41
  bool have_int_sources = false;
 
42
  if (!doing_solve_cw) {
 
43
    for (src_vol *sv = sources[ft2]; sv; sv = sv->next)
 
44
      if (sv->t->is_integrated) {
 
45
        have_int_sources = true;
 
46
        break;
 
47
      }
 
48
  }
 
49
 
 
50
  FOR_FT_COMPONENTS(ft2, dc) DOCMP {
 
51
    if (f[dc][cmp] && (pols[ft] || have_int_sources)) {
 
52
      if (!f_minus_p[dc][cmp]) f_minus_p[dc][cmp] = new realnum[gv.ntot()];
 
53
    }
 
54
    else if (f_minus_p[dc][cmp]) { // remove unneeded f_minus_p
 
55
      delete[] f_minus_p[dc][cmp];
 
56
      f_minus_p[dc][cmp] = 0;
 
57
    }
 
58
  }
 
59
  bool have_f_minus_p = false;
 
60
  FOR_FT_COMPONENTS(ft2, dc) if (f_minus_p[dc][0]) {
 
61
    have_f_minus_p = true;
 
62
    break;
 
63
  }
 
64
 
 
65
  const int ntot = s->gv.ntot();
 
66
 
 
67
  if (have_f_minus_p && doing_solve_cw)
 
68
    abort("dispersive materials are not yet implemented for solve_cw");
 
69
 
 
70
  //////////////////////////////////////////////////////////////////////////
 
71
  // First, initialize f_minus_p to D - P, if necessary
 
72
 
 
73
  if (have_f_minus_p) {
 
74
    if (pols[ft]) {
 
75
      FOR_FT_COMPONENTS(ft, ec) if (f[ec][0]) {
 
76
        for (polarization *np=pols[ft],*op=olpols[ft]; np; 
 
77
             np=np->next, op=op->next) {
 
78
          if (np->energy[ec] && op->energy[ec]) {
 
79
            if (is_real) for (int i = 0; i < ntot; ++i) {
 
80
              np->energy[ec][i] = op->energy[ec][i] +
 
81
                (0.5)*(np->P[ec][0][i] - op->P[ec][0][i])
 
82
                * f[ec][0][i];
 
83
            }
 
84
            else for (int i = 0; i < ntot; ++i) {
 
85
              np->energy[ec][i] = op->energy[ec][i] +
 
86
                (0.5)*(np->P[ec][0][i] - op->P[ec][0][i])
 
87
                * f[ec][0][i] +
 
88
                (0.5)*(np->P[ec][1][i] - op->P[ec][1][i])
 
89
                * f[ec][1][i];
 
90
            }
 
91
          }
 
92
        }
 
93
        component dc = direction_component(first_field_component(ft2),
 
94
                                           component_direction(ec));
 
95
        DOCMP {
 
96
          realnum * fmp = f_minus_p[dc][cmp];
 
97
          memcpy(fmp, f[dc][cmp], sizeof(realnum) * ntot);
 
98
          for (polarization *p = pols[ft]; p; p = p->next) {
 
99
            const realnum * P = p->P[ec][cmp];
 
100
            for (int i=0;i<ntot;i++) fmp[i] -= P[i];
 
101
          }
 
102
        }
 
103
      }
 
104
    }
 
105
    else {
 
106
      FOR_FT_COMPONENTS(ft2, dc) if (f[dc][0]) DOCMP
 
107
        memcpy(f_minus_p[dc][cmp], f[dc][cmp], ntot * sizeof(realnum));
 
108
    }
 
109
  }
 
110
 
 
111
  //////////////////////////////////////////////////////////////////////////
 
112
  // Next, subtract time-integrated sources (i.e. polarizations, not currents)
 
113
 
 
114
  if (have_f_minus_p && !doing_solve_cw) {
 
115
    for (src_vol *sv = sources[ft2]; sv; sv = sv->next) {  
 
116
      if (sv->t->is_integrated && f[sv->c][0] && ft == type(sv->c)) {
 
117
        component c = field_type_component(ft2, sv->c);
 
118
        for (int j = 0; j < sv->npts; ++j) { 
 
119
          const complex<double> A = sv->dipole(j);
 
120
          DOCMP {
 
121
            f_minus_p[c][cmp][sv->index[j]] -= 
 
122
              (cmp) ? imag(A) :  real(A);
 
123
          }
 
124
        }
 
125
      }
 
126
    }
 
127
  }
 
128
 
 
129
  //////////////////////////////////////////////////////////////////////////
 
130
  // Finally, compute E = chi1inv * D
 
131
  
 
132
  realnum *dmp[NUM_FIELD_COMPONENTS][2];
 
133
  if (have_f_minus_p) {
 
134
    FOR_FT_COMPONENTS(ft2,dc) DOCMP2 dmp[dc][cmp] = f_minus_p[dc][cmp];
 
135
  } else {
 
136
    FOR_FT_COMPONENTS(ft2,dc) DOCMP2 dmp[dc][cmp] = f[dc][cmp];
 
137
  }
 
138
 
 
139
  DOCMP FOR_FT_COMPONENTS(ft,ec) if (f[ec][cmp]) {
 
140
    if (type(ec) != ft) abort("bug in FOR_FT_COMPONENTS");
 
141
    component dc = field_type_component(ft2, ec);
 
142
    const direction d_ec = component_direction(ec);
 
143
    const int s_ec = gv.stride(d_ec) * (ft == H_stuff ? -1 : +1);
 
144
    const direction d_1 = cycle_direction(gv.dim, d_ec, 1);
 
145
    const component dc_1 = direction_component(dc,d_1);
 
146
    const int s_1 = gv.stride(d_1) * (ft == H_stuff ? -1 : +1);
 
147
    const direction d_2 = cycle_direction(gv.dim, d_ec, 2);
 
148
    const component dc_2 = direction_component(dc,d_2);
 
149
    const int s_2 = gv.stride(d_2) * (ft == H_stuff ? -1 : +1);
 
150
 
 
151
    direction dsigw0 = d_ec;
 
152
    direction dsigw = s->sigsize[dsigw0] > 1 ? dsigw0 : NO_DIRECTION;
 
153
 
 
154
    // lazily allocate any E/H fields that are needed (H==B initially)
 
155
    if (f[ec][cmp] == f[dc][cmp]
 
156
        && (s->chi1inv[ec][d_ec] || have_f_minus_p || dsigw != NO_DIRECTION)) {
 
157
      f[ec][cmp] = new realnum[gv.ntot()];
 
158
      memcpy(f[ec][cmp], f[dc][cmp], gv.ntot() * sizeof(realnum));
 
159
      allocated_eh = true;
 
160
    }
 
161
    
 
162
    // lazily allocate W auxiliary field
 
163
    if (!f_w[ec][cmp] && dsigw != NO_DIRECTION) {
 
164
      f_w[ec][cmp] = new realnum[gv.ntot()];
 
165
      memcpy(f_w[ec][cmp], f[ec][cmp], gv.ntot() * sizeof(realnum));
 
166
    }
 
167
 
 
168
    // for solve_cw, when W exists we get W and E from special variables
 
169
    if (f_w[ec][cmp] && skip_w_components) continue;
 
170
 
 
171
    if (f[ec][cmp] != f[dc][cmp])
 
172
      STEP_UPDATE_EDHB(f[ec][cmp], ec, gv, 
 
173
                       dmp[dc][cmp], dmp[dc_1][cmp], dmp[dc_2][cmp],
 
174
                       s->chi1inv[ec][d_ec], dmp[dc_1][cmp]?s->chi1inv[ec][d_1]:NULL, dmp[dc_2][cmp]?s->chi1inv[ec][d_2]:NULL,
 
175
                       s_ec, s_1, s_2, s->chi2[ec], s->chi3[ec],
 
176
                       f_w[ec][cmp], dsigw, s->sig[dsigw]);
 
177
  }
 
178
 
 
179
  /* Do annoying special cases for r=0 in cylindrical coords.  Note
 
180
     that this only really matters for field output; the Ez and Ep
 
181
     components at r=0 don't usually affect the fields elsewhere
 
182
     because of the form of Maxwell's equations in cylindrical coords. */
 
183
  // (FIXME: handle Kerr case?  Do we care about auxiliary PML fields here?)
 
184
  if (gv.dim == Dcyl && gv.origin_r() == 0.0)
 
185
    DOCMP FOR_FT_COMPONENTS(ft,ec) if (f[ec][cmp] && (ec == Ep || ec == Ez
 
186
                                                      || ec == Hr)) {
 
187
      component dc = field_type_component(ft2, ec);
 
188
      if (f[ec][cmp] == f[dc][cmp]) continue;
 
189
      const int yee_idx = gv.yee_index(ec);
 
190
      const int d_ec = component_direction(ec);
 
191
      const int sR = gv.stride(R), nZ = gv.num_direction(Z);
 
192
      realnum *E = f[ec][cmp];
 
193
      const realnum *D = have_f_minus_p ? f_minus_p[dc][cmp] : f[dc][cmp];
 
194
      const realnum *chi1inv = s->chi1inv[ec][d_ec];
 
195
      if (chi1inv)
 
196
        for (int iZ=0; iZ<nZ; iZ++) {
 
197
          const int i = yee_idx + iZ - sR;
 
198
          E[i] = chi1inv[i] * D[i];
 
199
        }
 
200
      else
 
201
        for (int iZ=0; iZ<nZ; iZ++) {
 
202
          const int i = yee_idx + iZ - sR;
 
203
          E[i] = D[i];
 
204
        }
 
205
    }
 
206
  
 
207
  return allocated_eh;
 
208
}
 
209
 
 
210
} // namespace meep