1 /*============================================================================
2 * Handle Maxwell module with CDO schemes
3 *============================================================================*/
4
5 /*
6 This file is part of Code_Saturne, a general-purpose CFD tool.
7
8 Copyright (C) 1998-2021 EDF S.A.
9
10 This program is free software; you can redistribute it and/or modify it under
11 the terms of the GNU General Public License as published by the Free Software
12 Foundation; either version 2 of the License, or (at your option) any later
13 version.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
18 details.
19
20 You should have received a copy of the GNU General Public License along with
21 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22 Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24
25 /*----------------------------------------------------------------------------*/
26
27 #include "cs_defs.h"
28
29 /*----------------------------------------------------------------------------
30 * Standard C library headers
31 *----------------------------------------------------------------------------*/
32
33 #include <assert.h>
34
35 /*----------------------------------------------------------------------------
36 * Local headers
37 *----------------------------------------------------------------------------*/
38
39 #include <bft_error.h>
40 #include <bft_mem.h>
41
42 #include "cs_hodge.h"
43 #include "cs_post.h"
44
45 /*----------------------------------------------------------------------------
46 * Header for the current file
47 *----------------------------------------------------------------------------*/
48
49 #include "cs_maxwell.h"
50
51 /*----------------------------------------------------------------------------*/
52
53 BEGIN_C_DECLS
54
55 /*=============================================================================
56 * Additional doxygen documentation
57 *============================================================================*/
58
59 /*!
60 \file cs_maxwell.h
61
62 \brief Structure and functions handling the Maxwell module dedicated to
63 the resolution of electro-magnetic equations
64
65 */
66
67 /*=============================================================================
68 * Local Macro definitions
69 *============================================================================*/
70
71 #define CS_MAXWELL_DBG 0
72
73 /*============================================================================
74 * Type definitions
75 *============================================================================*/
76
77 /* Set of parameters related to the Maxwell module */
78
79 struct _maxwell_t {
80
81 cs_flag_t model; /* Modelling for the Maxwell module */
82 cs_flag_t options; /* Flag dedicated to general options to handle
83 * the Maxwell module */
84 cs_flag_t post_flag; /* Flag dedicated to the post-processing
85 * of the Maxwell module */
86
87 /* Properties associated to this module */
88 cs_real_t e_perm_ref; /* Reference value of the electric
89 * permeability */
90 cs_property_t *e_permeability; /* Electric permeability oftenly denoted
91 * by epsilon */
92
93 cs_real_t m_perm_ref; /* Reference value of the magnetic
94 * permittivity */
95 cs_property_t *m_permittivity; /* Magnetic permittivity oftenly denoted
96 * by mu */
97
98 cs_property_t *conductivity; /* Conductivity in Ohm's law oftenly denoted
99 * by sigma */
100
101 /* Fields associated to this module */
102 cs_field_t *scal_pot; /* Scalar potential at vertices called
103 "electric_potential" */
104
105 cs_field_t *vect_pot; /* Vector potential */
106
107 cs_field_t *e_field; /* E: Electric field */
108 cs_real_t *e_field_array; /* E: Electric field along edges */
109
110 cs_field_t *d_field; /* D: Electric induction field
111 electric flux density (As/m^2) */
112 cs_real_t *d_field_array; /* D: Electric induction field across dual
113 faces */
114
115 cs_field_t *h_field; /* H = Magnetic field (faces) */
116 cs_real_t *h_field_array; /* H along dual edges */
117
118 cs_field_t *b_field; /* B = Magnetic induction field */
119 cs_real_t *b_field_array; /* B across faces */
120
121 cs_field_t *j_field; /* J = density flux field */
122 cs_real_t *j_field_array; /* J across dual faces */
123
124 cs_field_t *joule_effect; /* Joule effect power (W.m^-3) */
125 };
126
127 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
128
129 /*============================================================================
130 * Static global variables
131 *============================================================================*/
132
133 static const char _err_empty_maxwell[] =
134 " Stop execution.\n"
135 " The structure related to the Maxwell module is empty.\n"
136 " Please check your settings.\n";
137
138 static cs_maxwell_t *cs_maxwell_structure = NULL;
139
140 /* Vacuum magnetic permeability constant (H/m) */
141 static const cs_real_t cs_maxwell_vacuum_m_permeability = 1.25663706143591e-6;
142
143 /* Vacuum permittivity constant (F/m) */
144 static const cs_real_t cs_maxwell_vacuum_e_permittivity = 8.85418782e-12;
145
146 /*============================================================================
147 * Private static inline function prototypes
148 *============================================================================*/
149
150 /*============================================================================
151 * Private function prototypes
152 *============================================================================*/
153
154 /*----------------------------------------------------------------------------*/
155 /*!
156 * \brief Compute cellwise constant vector-valued approximation of fields
157 *
158 * \param[in] quant pointer to a cs_cdo_quantities_t structure
159 * \param[in] c2e pointer to a cs_adjacency_t structure
160 * \param[in] ep_values array of edge values
161 * \param[in] df_values array of dual face values
162 * \param[in, out] c_ep_values array vector-valued cell arrays
163 * \param[in, out] c_fd_values array vector-valued cell arrays
164 */
165 /*----------------------------------------------------------------------------*/
166
167 static void
_build_edge_based_vector_fields(const cs_cdo_quantities_t * quant,const cs_adjacency_t * c2e,const cs_real_t * ep_values,const cs_real_t * fd_values,cs_real_t * c_ep_values,cs_real_t * c_fd_values)168 _build_edge_based_vector_fields(const cs_cdo_quantities_t *quant,
169 const cs_adjacency_t *c2e,
170 const cs_real_t *ep_values,
171 const cs_real_t *fd_values,
172 cs_real_t *c_ep_values,
173 cs_real_t *c_fd_values)
174 {
175 assert(ep_values != NULL && fd_values != NULL);
176 assert(c_ep_values != NULL && c_fd_values != NULL);
177
178 memset(c_ep_values, 0, 3*quant->n_cells*sizeof(cs_real_t));
179 memset(c_fd_values, 0, 3*quant->n_cells*sizeof(cs_real_t));
180
181 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
182
183 cs_real_t *c_ep = c_ep_values + 3*c_id;
184 cs_real_t *c_fd = c_fd_values + 3*c_id;
185
186 for (cs_lnum_t j = c2e->idx[c_id]; j <c2e->idx[c_id+1]; j++) {
187
188 const cs_lnum_t e_id = c2e->ids[j];
189 const cs_real_t e_val = ep_values[e_id];
190 const cs_real_t *e_vect = quant->edge_vector + 3*e_id;
191 const cs_real_t *dface = quant->dface_normal + 3*j;
192 for (int k = 0; k < 3; k++) {
193 c_fd[k] += fd_values[j] * e_vect[k];
194 c_ep[k] += e_val * dface[k];
195 }
196
197 }
198
199 const double invvol = 1/quant->cell_vol[c_id];
200 for (int k = 0; k < 3; k++) {
201 c_fd[k] *= invvol;
202 c_ep[k] *= invvol;
203 }
204
205 } /* Loop on cells */
206
207 }
208
209 /*----------------------------------------------------------------------------*/
210 /*!
211 * \brief Compute cellwise constant vector-valued approximation of fields
212 * for fields associated to (1) primal faces and seen as a flux
213 * (fp_values) and to (2) dual edges and seen as a circulation
214 * (ed_values)
215 *
216 * \param[in] quant pointer to a cs_cdo_quantities_t structure
217 * \param[in] c2f pointer to a cs_adjacency_t structure
218 * \param[in] fp_values array of (primal) face values
219 * \param[in] ed_values array of dual edge values
220 * \param[in, out] c_fp_values array vector-valued cell arrays
221 * \param[in, out] c_ed_values array vector-valued cell arrays
222 */
223 /*----------------------------------------------------------------------------*/
224
225 static void
_build_face_based_vector_fields(const cs_cdo_quantities_t * quant,const cs_adjacency_t * c2f,const cs_real_t * fp_values,const cs_real_t * ed_values,cs_real_t * c_fp_values,cs_real_t * c_ed_values)226 _build_face_based_vector_fields(const cs_cdo_quantities_t *quant,
227 const cs_adjacency_t *c2f,
228 const cs_real_t *fp_values,
229 const cs_real_t *ed_values,
230 cs_real_t *c_fp_values,
231 cs_real_t *c_ed_values)
232 {
233 assert(fp_values != NULL && ed_values != NULL);
234 assert(c_fp_values != NULL && c_ed_values != NULL);
235
236 memset(c_fp_values, 0, 3*quant->n_cells*sizeof(cs_real_t));
237 memset(c_ed_values, 0, 3*quant->n_cells*sizeof(cs_real_t));
238
239 for (cs_lnum_t c_id = 0; c_id < quant->n_cells; c_id++) {
240
241 cs_real_t *c_fp = c_fp_values + 3*c_id;
242 cs_real_t *c_ed = c_ed_values + 3*c_id;
243
244 for (cs_lnum_t j = c2f->idx[c_id]; j <c2f->idx[c_id+1]; j++) {
245
246 const cs_lnum_t f_id = c2f->ids[j];
247 const cs_nvec3_t pfq = cs_quant_set_face_nvec(f_id, quant);
248 const cs_real_t ed_coef = ed_values[j] * pfq.meas;
249 const cs_real_t *ed_vect = quant->dedge_vector + 3*j;
250 const cs_real_t f_val = fp_values[f_id];
251
252 for (int k = 0; k < 3; k++) {
253 c_ed[k] += ed_coef * pfq.unitv[k];
254 c_fp[k] += f_val * ed_vect[k];
255 }
256
257 }
258
259 const double invvol = 1/quant->cell_vol[c_id];
260 for (int k = 0; k < 3; k++) {
261 c_ed[k] *= invvol;
262 c_fp[k] *= invvol;
263 }
264
265 } /* Loop on cells */
266
267 }
268
269 /*----------------------------------------------------------------------------*/
270 /*!
271 * \brief Create a structure dedicated to manage the Maxwell module
272 *
273 * \return a pointer to a new allocated cs_maxwell_t structure
274 */
275 /*----------------------------------------------------------------------------*/
276
277 static cs_maxwell_t *
_maxwell_create(void)278 _maxwell_create(void)
279 {
280 cs_maxwell_t *mxl = NULL;
281
282 BFT_MALLOC(mxl, 1, cs_maxwell_t);
283
284 /* Default initialization */
285 mxl->model = 0;
286 mxl->options = 0;
287 mxl->post_flag = 0;
288
289 /* Properties */
290 mxl->e_perm_ref = cs_maxwell_vacuum_e_permittivity;
291 mxl->e_permeability = NULL;
292
293 mxl->m_perm_ref = cs_maxwell_vacuum_m_permeability;
294 mxl->m_permittivity = NULL;
295
296 mxl->conductivity = NULL;
297
298 /* Fields and related arrays */
299 mxl->scal_pot = NULL;
300
301 mxl->vect_pot = NULL;
302
303 mxl->e_field = NULL;
304 mxl->e_field_array = NULL; /* different location that e_field */
305
306 mxl->d_field = NULL;
307 mxl->d_field_array = NULL; /* different location that d_field */
308
309 mxl->h_field = NULL;
310 mxl->h_field_array = NULL; /* different location that h_field */
311
312 mxl->b_field = NULL;
313 mxl->b_field_array = NULL; /* different location that b_field */
314
315 mxl->j_field = NULL;
316 mxl->j_field_array = NULL; /* different location that j_field */
317
318 /* Additional quantities (source terms for instance) */
319 mxl->joule_effect = NULL;
320
321 return mxl;
322 }
323
324 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
325
326 /*============================================================================
327 * Public function prototypes
328 *============================================================================*/
329
330 /*----------------------------------------------------------------------------*/
331 /*!
332 * \brief Test if the computation of Maxwell equations is activated
333 */
334 /*----------------------------------------------------------------------------*/
335
336 bool
cs_maxwell_is_activated(void)337 cs_maxwell_is_activated(void)
338 {
339 if (cs_maxwell_structure == NULL)
340 return false;
341 else
342 return true;
343 }
344
345 /*----------------------------------------------------------------------------*/
346 /*!
347 * \brief Activate the future computation of the Maxwell equations
348 *
349 * \param[in] model type of modelling
350 * \param[in] options flag to handle optional parameters
351 *
352 * \return a pointer to a new allocated Maxwell structure
353 */
354 /*----------------------------------------------------------------------------*/
355
356 cs_maxwell_t *
cs_maxwell_activate(cs_flag_t model,cs_flag_t options)357 cs_maxwell_activate(cs_flag_t model,
358 cs_flag_t options)
359 {
360 if (model < 1)
361 bft_error(__FILE__, __LINE__, 0,
362 " %s: Invalid modelling. Model = %d\n", __func__, model);
363
364 /* Allocate an empty structure */
365 cs_maxwell_t *mxl = _maxwell_create();
366
367 /* Set members of the structure according to the given settings */
368 mxl->model = model;
369 mxl->options = options;
370
371 if (model & CS_MAXWELL_MODEL_ELECTROSTATIC) {
372
373 cs_equation_t *e_static = cs_equation_add(CS_MAXWELL_ESTATIC_EQNAME,
374 "electric_potential",
375 CS_EQUATION_TYPE_MAXWELL,
376 1,
377 CS_PARAM_BC_HMG_NEUMANN);
378 cs_equation_param_t *eqp = cs_equation_get_param(e_static);
379
380 mxl->e_permeability = cs_property_add("electric_permeability",
381 CS_PROPERTY_ISO);
382
383 /* By default, set the reference permeability */
384 cs_property_def_iso_by_value(mxl->e_permeability, NULL, mxl->e_perm_ref);
385
386 cs_equation_add_diffusion(eqp, mxl->e_permeability);
387
388 /* Should be symmetric */
389 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_vb");
390 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_ALGO, "bubble");
391 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_COEF, "frac23");
392 cs_equation_param_set(eqp, CS_EQKEY_SOLVER_FAMILY, "cs");
393 cs_equation_param_set(eqp, CS_EQKEY_PRECOND, "amg");
394 cs_equation_param_set(eqp, CS_EQKEY_ITSOL, "cg");
395 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_EPS, "1e-6");
396 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_RESNORM_TYPE, "filtered");
397
398 }
399
400 if (model & CS_MAXWELL_MODEL_MAGNETOSTATIC) {
401
402 cs_equation_t *m_static = cs_equation_add(CS_MAXWELL_MSTATIC_EQNAME,
403 "magnetic_potential",
404 CS_EQUATION_TYPE_MAXWELL,
405 3,
406 CS_PARAM_BC_HMG_DIRICHLET);
407
408 cs_equation_param_t *eqp = cs_equation_get_param(m_static);
409
410 mxl->m_permittivity = cs_property_add("magnetic_permittivity",
411 CS_PROPERTY_ISO);
412
413 /* By default, set the reference permeability */
414 cs_property_def_iso_by_value(mxl->m_permittivity, NULL, mxl->m_perm_ref);
415
416 cs_equation_add_curlcurl(eqp, mxl->m_permittivity,
417 1); /* Inverse of the property is requested */
418
419 /* Should be symmetric */
420 cs_equation_param_set(eqp, CS_EQKEY_SPACE_SCHEME, "cdo_eb");
421 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_ALGO, "cost");
422 cs_equation_param_set(eqp, CS_EQKEY_HODGE_DIFF_COEF, "dga");
423 cs_equation_param_set(eqp, CS_EQKEY_SOLVER_FAMILY, "cs");
424 cs_equation_param_set(eqp, CS_EQKEY_PRECOND, "amg");
425 cs_equation_param_set(eqp, CS_EQKEY_ITSOL, "cg");
426 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_EPS, "1e-6");
427 cs_equation_param_set(eqp, CS_EQKEY_ITSOL_RESNORM_TYPE, "filtered");
428
429 }
430
431 cs_maxwell_structure = mxl;
432
433 return mxl;
434 }
435
436 /*----------------------------------------------------------------------------*/
437 /*!
438 * \brief Free the main structure related to the Maxwell module
439 *
440 * \return a NULL pointer
441 */
442 /*----------------------------------------------------------------------------*/
443
444 cs_maxwell_t *
cs_maxwell_destroy_all(void)445 cs_maxwell_destroy_all(void)
446 {
447 if (cs_maxwell_structure == NULL)
448 return NULL;
449
450 cs_maxwell_t *mxl = cs_maxwell_structure;
451
452 /* The lifecycle of properties and fields is not managed by the current
453 structure.
454 Free only arrays which are owned by this structure */
455
456 BFT_FREE(mxl->e_field_array);
457 BFT_FREE(mxl->d_field_array);
458 BFT_FREE(mxl->h_field_array);
459 BFT_FREE(mxl->b_field_array);
460 BFT_FREE(mxl->j_field_array);
461
462 BFT_FREE(mxl);
463
464 return NULL;
465 }
466
467 /*----------------------------------------------------------------------------*/
468 /*!
469 * \brief Setup equations/properties related to the Maxwell module
470 */
471 /*----------------------------------------------------------------------------*/
472
473 void
cs_maxwell_init_setup(void)474 cs_maxwell_init_setup(void)
475 {
476 cs_maxwell_t *mxl = cs_maxwell_structure;
477
478 /* Sanity checks */
479 if (mxl == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_maxwell));
480
481 const int field_mask = CS_FIELD_INTENSIVE | CS_FIELD_CDO;
482 const int log_key = cs_field_key_id("log");
483 const int post_key = cs_field_key_id("post_vis");
484
485 if (mxl->model & CS_MAXWELL_MODEL_ELECTROSTATIC) {
486
487 mxl->e_field = cs_field_create(CS_MAXWELL_EFIELD_NAME,
488 field_mask,
489 CS_MESH_LOCATION_CELLS,
490 3,
491 true);
492
493 cs_field_set_key_int(mxl->e_field, log_key, 1);
494 cs_field_set_key_int(mxl->e_field, post_key, 1);
495
496 mxl->d_field = cs_field_create(CS_MAXWELL_DFIELD_NAME,
497 field_mask,
498 CS_MESH_LOCATION_CELLS,
499 3,
500 true);
501
502 cs_field_set_key_int(mxl->d_field, log_key, 1);
503 cs_field_set_key_int(mxl->d_field, post_key, 1);
504
505 }
506
507 if (mxl->model & CS_MAXWELL_MODEL_MAGNETOSTATIC) {
508
509 mxl->b_field = cs_field_create(CS_MAXWELL_BFIELD_NAME,
510 field_mask,
511 CS_MESH_LOCATION_CELLS,
512 3,
513 true);
514
515 cs_field_set_key_int(mxl->b_field, log_key, 1);
516 cs_field_set_key_int(mxl->b_field, post_key, 1);
517
518 mxl->h_field = cs_field_create(CS_MAXWELL_MFIELD_NAME,
519 field_mask,
520 CS_MESH_LOCATION_CELLS,
521 3,
522 true);
523
524 cs_field_set_key_int(mxl->h_field, log_key, 1);
525 cs_field_set_key_int(mxl->h_field, post_key, 1);
526
527 }
528
529 /* Optional settings */
530 if (mxl->options & CS_MAXWELL_JOULE_EFFECT) {
531
532 mxl->joule_effect = cs_field_create(CS_MAXWELL_JEFFECT_NAME,
533 field_mask,
534 CS_MESH_LOCATION_CELLS,
535 1,
536 true);
537 cs_field_set_key_int(mxl->joule_effect, log_key, 1);
538 cs_field_set_key_int(mxl->joule_effect, post_key, 1);
539
540 }
541
542 /* Add default post-processing related to the Maxwell module */
543 cs_post_add_time_mesh_dep_output(cs_maxwell_extra_post, mxl);
544 }
545
546 /*----------------------------------------------------------------------------*/
547 /*!
548 * \brief Finalize the setup stage for equations related to the Maxwell module
549 *
550 * \param[in] connect pointer to a cs_cdo_connect_t structure
551 * \param[in] quant pointer to a cs_cdo_quantities_t structure
552 */
553 /*----------------------------------------------------------------------------*/
554
555 void
cs_maxwell_finalize_setup(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)556 cs_maxwell_finalize_setup(const cs_cdo_connect_t *connect,
557 const cs_cdo_quantities_t *quant)
558 {
559 cs_maxwell_t *mxl = cs_maxwell_structure;
560
561 /* Sanity checks */
562 if (mxl == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_maxwell));
563
564 if (mxl->model & CS_MAXWELL_MODEL_ELECTROSTATIC) {
565
566 cs_equation_t *es_eq = cs_equation_by_name(CS_MAXWELL_ESTATIC_EQNAME);
567 assert(cs_equation_get_space_scheme(es_eq) == CS_SPACE_SCHEME_CDOVB);
568
569 mxl->scal_pot = cs_equation_get_field(es_eq);
570
571 /* Electric field array along edges */
572 BFT_MALLOC(mxl->e_field_array, quant->n_edges, cs_real_t);
573 memset(mxl->e_field_array, 0, quant->n_edges*sizeof(cs_real_t));
574
575 /* Electric induction (flux density) across dual faces */
576 const cs_adjacency_t *c2e = connect->c2e;
577 const cs_lnum_t array_size = c2e->idx[quant->n_cells];
578 BFT_MALLOC(mxl->d_field_array, array_size, cs_real_t);
579 memset(mxl->d_field_array, 0, array_size*sizeof(cs_real_t));
580
581 }
582
583 if (mxl->model & CS_MAXWELL_MODEL_MAGNETOSTATIC) {
584
585 cs_equation_t *ms_eq = cs_equation_by_name(CS_MAXWELL_MSTATIC_EQNAME);
586 assert(cs_equation_get_space_scheme(ms_eq) == CS_SPACE_SCHEME_CDOEB);
587
588 mxl->vect_pot = cs_equation_get_field(ms_eq);
589
590 /* Magnetic field array along dual edges */
591 const cs_adjacency_t *c2f = connect->c2f;
592 const cs_lnum_t array_size = c2f->idx[quant->n_cells];
593 BFT_MALLOC(mxl->h_field_array, array_size, cs_real_t);
594 memset(mxl->h_field_array, 0, array_size*sizeof(cs_real_t));
595
596 /* Magnetic induction (flux density) across primal faces */
597 BFT_MALLOC(mxl->b_field_array, quant->n_faces, cs_real_t);
598 memset(mxl->b_field_array, 0, quant->n_faces*sizeof(cs_real_t));
599
600 }
601
602 }
603
604 /*----------------------------------------------------------------------------*/
605 /*!
606 * \brief Log a summary of the Maxwell module
607 */
608 /*----------------------------------------------------------------------------*/
609
610 void
cs_maxwell_log_setup(void)611 cs_maxwell_log_setup(void)
612 {
613 cs_maxwell_t *mxl = cs_maxwell_structure;
614
615 if (mxl == NULL)
616 return;
617
618 cs_log_printf(CS_LOG_SETUP, "\nSummary of the Maxwell module\n");
619 cs_log_printf(CS_LOG_SETUP, "%s\n", cs_sep_h1);
620
621 cs_log_printf(CS_LOG_SETUP, " * Maxwell | Model:");
622 if (mxl->model & CS_MAXWELL_MODEL_ELECTROSTATIC)
623 cs_log_printf(CS_LOG_SETUP, " Electro-static");
624 if (mxl->model & CS_MAXWELL_MODEL_MAGNETOSTATIC)
625 cs_log_printf(CS_LOG_SETUP, "+ Magneto-static");
626 cs_log_printf(CS_LOG_SETUP, "\n");
627
628 if (mxl->options & CS_MAXWELL_JOULE_EFFECT)
629 cs_log_printf(CS_LOG_SETUP, " * Maxwell | Joule effect\n");
630
631 }
632
633 /*----------------------------------------------------------------------------*/
634 /*!
635 * \brief Solve if needed the steady-state equations related to the Maxwell
636 * module
637 *
638 * \param[in] mesh pointer to a cs_mesh_t structure
639 * \param[in] time_step pointer to a cs_time_step_t structure
640 * \param[in] connect pointer to a cs_cdo_connect_t structure
641 * \param[in] quant pointer to a cs_cdo_quantities_t structure
642 */
643 /*----------------------------------------------------------------------------*/
644
645 void
cs_maxwell_compute_steady_state(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)646 cs_maxwell_compute_steady_state(const cs_mesh_t *mesh,
647 const cs_time_step_t *time_step,
648 const cs_cdo_connect_t *connect,
649 const cs_cdo_quantities_t *quant)
650 {
651 cs_maxwell_t *mxl = cs_maxwell_structure;
652
653 /* Sanity checks */
654 if (mxl == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_maxwell));
655
656 if (mxl->model & CS_MAXWELL_MODEL_ELECTROSTATIC) {
657
658 cs_equation_t *es_eq = cs_equation_by_name(CS_MAXWELL_ESTATIC_EQNAME);
659
660 /* Sanity checks */
661 assert(es_eq != NULL);
662 assert(cs_equation_uses_new_mechanism(es_eq));
663
664 cs_equation_solve_steady_state(mesh, es_eq);
665
666 }
667
668 if (mxl->model & CS_MAXWELL_MODEL_MAGNETOSTATIC) {
669
670 cs_equation_t *ms_eq = cs_equation_by_name(CS_MAXWELL_MSTATIC_EQNAME);
671
672 /* Sanity checks */
673 assert(ms_eq != NULL);
674 assert(cs_equation_uses_new_mechanism(ms_eq));
675
676 cs_equation_solve_steady_state(mesh, ms_eq);
677
678 }
679
680 /* Update fields and properties which are related to solved variables */
681 cs_maxwell_update(mesh, connect, quant, time_step,
682 true); /* operate current to previous ? */
683 }
684
685 /*----------------------------------------------------------------------------*/
686 /*!
687 * \brief Solve equations related to the Maxwell module
688 *
689 * \param[in] mesh pointer to a cs_mesh_t structure
690 * \param[in] time_step pointer to a cs_time_step_t structure
691 * \param[in] connect pointer to a cs_cdo_connect_t structure
692 * \param[in] quant pointer to a cs_cdo_quantities_t structure
693 */
694 /*----------------------------------------------------------------------------*/
695
696 void
cs_maxwell_compute(const cs_mesh_t * mesh,const cs_time_step_t * time_step,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)697 cs_maxwell_compute(const cs_mesh_t *mesh,
698 const cs_time_step_t *time_step,
699 const cs_cdo_connect_t *connect,
700 const cs_cdo_quantities_t *quant)
701 {
702 cs_maxwell_t *mxl = cs_maxwell_structure;
703
704 /* Sanity checks */
705 if (mxl == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_maxwell));
706
707 /* Add equations to be solved at each time step */
708
709 /* Update fields and properties which are related to solved variables */
710 cs_maxwell_update(mesh, connect, quant, time_step,
711 true); /* operate current to previous ? */
712 }
713
714 /*----------------------------------------------------------------------------*/
715 /*!
716 * \brief Update/initialize the Maxwell module according to the settings
717 *
718 * \param[in] mesh pointer to a cs_mesh_t structure
719 * \param[in] connect pointer to a cs_cdo_connect_t structure
720 * \param[in] quant pointer to a cs_cdo_quantities_t structure
721 * \param[in] ts pointer to a cs_time_step_t structure
722 * \param[in] cur2prev true or false
723 */
724 /*----------------------------------------------------------------------------*/
725
726 void
cs_maxwell_update(const cs_mesh_t * mesh,const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant,const cs_time_step_t * ts,bool cur2prev)727 cs_maxwell_update(const cs_mesh_t *mesh,
728 const cs_cdo_connect_t *connect,
729 const cs_cdo_quantities_t *quant,
730 const cs_time_step_t *ts,
731 bool cur2prev)
732 {
733 CS_UNUSED(mesh);
734
735 cs_maxwell_t *mxl = cs_maxwell_structure;
736
737 /* Sanity checks */
738 if (mxl == NULL) bft_error(__FILE__, __LINE__, 0, _(_err_empty_maxwell));
739
740 if (mxl->model & CS_MAXWELL_MODEL_ELECTROSTATIC) {
741
742 cs_equation_t *es_eq = cs_equation_by_name(CS_MAXWELL_ESTATIC_EQNAME);
743
744 /* Retrieve the scalar electric potential */
745 const cs_real_t *pot = cs_equation_get_vertex_values(es_eq, false);
746
747 /* Compute the electric field: E = -grad(scal_pot) */
748 const cs_adjacency_t *e2v = connect->e2v;
749 for (cs_lnum_t e = 0; e < quant->n_edges; e++) {
750
751 const cs_lnum_t *v_ids = e2v->ids + 2*e;
752
753 /* E = -grad(scal_pot) */
754 mxl->e_field_array[e] = e2v->sgn[2*e]*(pot[v_ids[0]] - pot[v_ids[1]]);
755
756 } /* Loop on edges */
757
758 cs_equation_compute_diff_flux_cellwise(es_eq,
759 cs_flag_dual_face_byc,
760 ts->t_cur,
761 mxl->d_field_array);
762
763 /* Update related vector-valued fields at cell centers */
764 if (cur2prev) {
765 cs_field_current_to_previous(mxl->e_field);
766 cs_field_current_to_previous(mxl->d_field);
767 }
768
769 _build_edge_based_vector_fields(quant,
770 connect->c2e,
771 mxl->e_field_array, /* in */
772 mxl->d_field_array, /* in */
773 mxl->e_field->val,
774 mxl->d_field->val);
775
776 } /* Electrostatic updates */
777
778 if (mxl->model & CS_MAXWELL_MODEL_MAGNETOSTATIC) {
779
780 cs_equation_t *ms_eq = cs_equation_by_name(CS_MAXWELL_MSTATIC_EQNAME);
781 cs_equation_param_t *ms_eqp = cs_equation_get_param(ms_eq);
782
783 /* Retrieve the scalar electric potential */
784 const cs_real_t *pot = cs_equation_get_edge_values(ms_eq, false);
785
786 /* Compute the magnetic induction field: B = curl(vect_pot) */
787 cs_cdo_connect_discrete_curl(connect, pot, &(mxl->b_field_array));
788
789 /* Compute the magnetic field using Hodge operator */
790 cs_hodge_circulation_from_flux(connect, quant, ts->t_cur,
791 ms_eqp->curlcurl_hodgep,
792 ms_eqp->curlcurl_property,
793 mxl->b_field_array,
794 mxl->h_field_array);
795
796
797 /* Update related vector-valued fields at cell centers */
798 if (cur2prev) {
799 cs_field_current_to_previous(mxl->b_field);
800 cs_field_current_to_previous(mxl->h_field);
801 }
802
803 _build_face_based_vector_fields(quant,
804 connect->c2f,
805 mxl->b_field_array, /* in */
806 mxl->h_field_array, /* in */
807 mxl->b_field->val,
808 mxl->h_field->val);
809
810 } /* Magnetostatic updates */
811
812 }
813
814 /*----------------------------------------------------------------------------*/
815 /*!
816 * \brief Predefined extra-operations for the Maxwell module
817 *
818 * \param[in] connect pointer to a cs_cdo_connect_t structure
819 * \param[in] quant pointer to a cs_cdo_quantities_t structure
820 */
821 /*----------------------------------------------------------------------------*/
822
823 void
cs_maxwell_extra_op(const cs_cdo_connect_t * connect,const cs_cdo_quantities_t * quant)824 cs_maxwell_extra_op(const cs_cdo_connect_t *connect,
825 const cs_cdo_quantities_t *quant)
826 {
827 CS_UNUSED(connect);
828 CS_UNUSED(quant);
829
830 cs_maxwell_t *mxl = cs_maxwell_structure;
831
832 if (mxl == NULL)
833 return;
834
835 /* TODO */
836 }
837
838 /*----------------------------------------------------------------------------*/
839 /*!
840 * \brief Predefined post-processing output for the Maxwell module.
841 * Prototype of this function is fixed since it is a function pointer
842 * defined in cs_post.h (\ref cs_post_time_mesh_dep_output_t)
843 *
844 * \param[in, out] input pointer to an optional structure (here a
845 * cs_gwf_t structure)
846 * \param[in] mesh_id id of the output mesh for the current call
847 * \param[in] cat_id category id of the output mesh for this call
848 * \param[in] ent_flag indicate global presence of cells (ent_flag[0]),
849 * interior faces (ent_flag[1]), boundary faces
850 * (ent_flag[2]), particles (ent_flag[3]) or probes
851 * (ent_flag[4])
852 * \param[in] n_cells local number of cells of post_mesh
853 * \param[in] n_i_faces local number of interior faces of post_mesh
854 * \param[in] n_b_faces local number of boundary faces of post_mesh
855 * \param[in] cell_ids list of cells (0 to n-1)
856 * \param[in] i_face_ids list of interior faces (0 to n-1)
857 * \param[in] b_face_ids list of boundary faces (0 to n-1)
858 * \param[in] time_step pointer to a cs_time_step_t struct.
859 */
860 /*----------------------------------------------------------------------------*/
861
862 void
cs_maxwell_extra_post(void * input,int mesh_id,int cat_id,int ent_flag[5],cs_lnum_t n_cells,cs_lnum_t n_i_faces,cs_lnum_t n_b_faces,const cs_lnum_t cell_ids[],const cs_lnum_t i_face_ids[],const cs_lnum_t b_face_ids[],const cs_time_step_t * time_step)863 cs_maxwell_extra_post(void *input,
864 int mesh_id,
865 int cat_id,
866 int ent_flag[5],
867 cs_lnum_t n_cells,
868 cs_lnum_t n_i_faces,
869 cs_lnum_t n_b_faces,
870 const cs_lnum_t cell_ids[],
871 const cs_lnum_t i_face_ids[],
872 const cs_lnum_t b_face_ids[],
873 const cs_time_step_t *time_step)
874 {
875 CS_UNUSED(mesh_id);
876 CS_UNUSED(cat_id);
877 CS_UNUSED(ent_flag);
878 CS_UNUSED(n_cells);
879 CS_UNUSED(n_i_faces);
880 CS_UNUSED(n_b_faces);
881 CS_UNUSED(cell_ids);
882 CS_UNUSED(i_face_ids);
883 CS_UNUSED(b_face_ids);
884 CS_UNUSED(time_step);
885
886 if (input == NULL)
887 return;
888
889 /* TODO */
890 }
891
892 /*----------------------------------------------------------------------------*/
893
894 END_C_DECLS
895