BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
powerlaws.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California, through
4 * Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
53
54#ifndef BELFEM_POWERLAWS_HPP
55#define BELFEM_POWERLAWS_HPP
56
57#include "cl_Material.hpp"
58
59namespace belfem
60{
61//------------------------------------------------------------------------------
62// Power-law resistivity model
63//------------------------------------------------------------------------------
64
65 inline real
66 Material::jc_eval( const real T, const real normB, const real angleNxB ) const
67 {
68 return ( mJcFunction != nullptr )
69 ? ( mJcFunction->depends_on( material::JcParameter::T )
70 ? mJcFunction->eval( normB, angleNxB, T )
71 : mJcFunction->eval( normB, angleNxB ) )
72 : this->constant_property( MaterialProperty::jc ) ;
73 }
74
75 inline real
76 Material::n_eval_raw( const real T, const real normB, const real angleNxB ) const
77 {
78 return ( mNFunction != nullptr )
79 ? ( mNFunction->depends_on( material::JcParameter::T )
80 ? mNFunction->eval( normB, angleNxB, T )
81 : mNFunction->eval( normB, angleNxB ) )
82 : this->constant_property( MaterialProperty::n ) ;
83 }
84
85 inline real
86 Material::n_eval( const real T, const real normB, const real angleNxB ) const
87 {
88 // ohmic floor ( 2026-08-27 ): measured tables soften through n = 1
89 // near T_crit; below that the raw power law is sub-ohmic and its
90 // J -> 0 limit flips. At n = 1 the tape is a plain resistor
91 // ec/jc. dn_eval_dB / dn_eval_dT return 0 while the floor binds,
92 // so the tangents differentiate the SAME clamped law
93 return std::max( this->n_eval_raw( T, normB, angleNxB ), 1.0 ) ;
94 }
95
96 // derivative routing mirrors jc_eval / n_eval: with no function attached
97 // the value is a constant and its derivative is EXACTLY zero — which also
98 // makes constant-jc decks bit-identical, since the Newton consumer
99 // early-outs on zero ( audited 2026-08-13 )
100
101 inline real
102 Material::djc_eval_dB( const real T, const real normB, const real angleNxB ) const
103 {
104 return ( mJcFunction != nullptr )
105 ? mJcFunction->deval_dB( normB, angleNxB, T )
106 : 0.0 ;
107 }
108
109 inline real
110 Material::dn_eval_dB( const real T, const real normB, const real angleNxB ) const
111 {
112 if ( mNFunction == nullptr ) return 0.0 ;
113
114 // n_eval clamps at 1; the clamped law is locally constant there
115 if ( ! ( this->n_eval_raw( T, normB, angleNxB ) > 1.0 ) ) return 0.0 ;
116
117 return mNFunction->deval_dB( normB, angleNxB, T ) ;
118 }
119
120 // T-leg twins ( audited 2026-08-13 ): same null-check
121 // routing. Exact zero for constants, ModifiedKim and 2-arg functions;
122 // for a 3-arg ( T-dependent ) UserDefined WITHOUT a deval_dT override
123 // the base-class 0.0 is a conservative fallback, not an exact derivative
124
125 inline real
126 Material::djc_eval_dT( const real T, const real normB, const real angleNxB ) const
127 {
128 return ( mJcFunction != nullptr )
129 ? mJcFunction->deval_dT( normB, angleNxB, T )
130 : 0.0 ;
131 }
132
133 inline real
134 Material::dn_eval_dT( const real T, const real normB, const real angleNxB ) const
135 {
136 if ( mNFunction == nullptr ) return 0.0 ;
137
138 // n_eval clamps at 1; the clamped law is locally constant there
139 if ( ! ( this->n_eval_raw( T, normB, angleNxB ) > 1.0 ) ) return 0.0 ;
140
141 return mNFunction->deval_dT( normB, angleNxB, T ) ;
142 }
143
183 inline real
184 Material::rho_powerlaw( const real normJ ) const
185 {
186 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
187 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
188
189 real ec = this->constant_property( MaterialProperty::ec ) ;
190 real jc = this->constant_property( MaterialProperty::jc ) ;
191 real n = this->constant_property( MaterialProperty::n ) ;
192
193 real rhon = this->rho( gTbulk ) ;
194 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
195
196 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
197 }
198
204 inline real
205 Material::rho_powerlaw( const real normJ, const real x, const real y, const real z, const real t ) const
206 {
207 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
208 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
209
210 real ec = this->constant_property( MaterialProperty::ec ) ;
211 real jc = this->constant_property( MaterialProperty::jc )*((this->mDefectFunction) (x,y,z,t)) ;
212 real n = this->constant_property( MaterialProperty::n ) ;
213
214 real rhon = this->rho( gTbulk ) ;
215 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
216
217 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
218 }
219
225 inline real
226 Material::rho_powerlaw( const real normJ, const real normB, const real angleNxB ) const
227 {
228 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
229 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
230
231
232 BELFEM_ASSERT(
233 mJcFunction->depends_on( material::JcParameter::normB ) &&
234 mJcFunction->depends_on( material::JcParameter::angleNxB ) &&
235 ! mJcFunction->depends_on( material::JcParameter::T ), "wrong powerlaw for material %s", mLabel.c_str() );
236
237
238 real ec = this->constant_property( MaterialProperty::ec ) ;
239 real jc = mJcFunction->eval( normB, angleNxB ) ;
240 real n = mNFunction->eval( normB, angleNxB ) ;
241
242 real rhon = this->rho( gTbulk ) ;
243 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
244
245 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
246 }
247
252 inline real
253 Material::rho_powerlaw( const real normJ, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
254 {
255 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
256 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
257
258
259 BELFEM_ASSERT(
260 mJcFunction->depends_on( material::JcParameter::normB ) &&
261 mJcFunction->depends_on( material::JcParameter::angleNxB ) &&
262 ! mJcFunction->depends_on( material::JcParameter::T ), "wrong powerlaw for material %s", mLabel.c_str() );
263
264
265 real ec = this->constant_property( MaterialProperty::ec ) ;
266 real jc = mJcFunction->eval( normB, angleNxB )*((this->mDefectFunction) (x,y,z,t)) ; ;
267 real n = mNFunction->eval( normB, angleNxB ) ;
268
269 real rhon = this->rho( gTbulk ) ;
270 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
271
272 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
273 }
274
282 inline real
283 Material::rho_powerlaw( const real normJ, const real T, const real normB, const real angleNxB ) const
284 {
285 // dependency-routed evaluation with constants-fallback ( O1 policy )
286 real ec = this->constant_property( MaterialProperty::ec ) ;
287 real jc = this->jc_eval( T, normB, angleNxB ) ;
288 real n = this->n_eval( T, normB, angleNxB ) ;
289
290 real rhon = this->rho( T ) ;
291 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
292
293 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
294 }
295
301 inline real
302 Material::rho_powerlaw( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
303 {
304 // dependency-routed evaluation with constants-fallback ( O1 policy )
305 real ec = this->constant_property( MaterialProperty::ec ) ;
306 real jc = this->jc_eval( T, normB, angleNxB ) * ((this->mDefectFunction) (x,y,z,t)) ;
307 real n = this->n_eval( T, normB, angleNxB ) ;
308
309 real rhon = this->rho( T ) ;
310 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
311
312 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
313 }
314
323 inline real
324 Material::rho_powerlaw( const real normJ, const real T ) const
325 {
326
327 real ec = this->constant_property( MaterialProperty::ec ) ;
328 real jc = this->jc_custom( T ) ;
329 real n = this->n_custom( T ) ;
330
331 real rhon = this->rho( T ) ;
332 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
333
334 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
335 }
336
341 inline real
342 Material::rho_powerlaw( const real normJ, const real T, const real x, const real y, const real z, const real t ) const
343 {
344 real ec = this->constant_property( MaterialProperty::ec ) ;
345 real jc = this->jc_custom( T ) * ((this->mDefectFunction) (x,y,z,t)) ;
346 real n = this->n_custom( T ) ;
347
348 real rhon = this->rho( T ) ;
349 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
350
351 return 1.0/((1.0/rhon) + (1.0/rhoPL)) ;
352 }
353
354//------------------------------------------------------------------------------
355// Three-regime piecewise resistivity model
356//------------------------------------------------------------------------------
357
411 inline real
412 Material::rho_piecewise( const real normJ ) const
413 {
414 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
415 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
416
417 real rhon = this->rho( gTbulk ) ;
418
419 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) ) // (assumed critical temperature)
420 {
421 return rhon ;
422 }
423
424 real ec = this->constant_property( MaterialProperty::ec ) ;
425 real jc = this->constant_property( MaterialProperty::jc ) ;
426 real n = this->constant_property( MaterialProperty::n ) ;
427
428 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
429 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
430
431 //Power law resistivity
432 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
433
434 //Power-law regime
435 real j1 = jc * std::pow(10.0,2.5/n) ;
436 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
437
438 if (normJ <= j1)
439 {
440 return rhoPL ;
441 }
442
443 //Normal regime
444 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
445 real rho3 = rhon ;
446
447 if ( normJ > j3 )
448 {
449 return rhon ;
450 }
451
452 //Flux-flow regime
453 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
454 real rho2 = rhon ;
455
456 //Bezier control points
457 real logj1 = std::log10(j1) ;
458 real logj2 = std::log10(j2) ;
459
460 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
461 real b = logj1 - logj2 ;
462 real c = std::log10(normJ) - logj1 ;
463
464 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
465 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
466 ( double ) std::abs( a ), mLabel.c_str() ) ;
467
468 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
469
470 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
471
472 return rhoFF ;
473 }
474
480 inline real
481 Material::rho_piecewise( const real normJ, const real x, const real y, const real z, const real t ) const
482 {
483 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
484 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
485
486 real rhon = this->rho( gTbulk ) ;
487
488 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
489 {
490 return rhon ;
491 }
492
493 real ec = this->constant_property( MaterialProperty::ec ) ;
494 real jc = this->constant_property( MaterialProperty::jc )*((this->mDefectFunction) (x,y,z,t)) ;
495 real n = this->constant_property( MaterialProperty::n ) ;
496
497 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
498 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
499
500 //Power law resistivity
501 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
502
503 //Power-law regime
504 real j1 = jc * std::pow(10.0,2.5/n) ;
505 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
506
507 if (normJ <= j1)
508 {
509 return rhoPL ;
510 }
511
512 //Normal regime
513 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
514 real rho3 = rhon ;
515
516 if ( normJ > j3 )
517 {
518 return rhon ;
519 }
520
521 //Flux-flow regime
522 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
523 real rho2 = rhon ;
524
525 //Bezier control points
526 real logj1 = std::log10(j1) ;
527 real logj2 = std::log10(j2) ;
528
529 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
530 real b = logj1 - logj2 ;
531 real c = std::log10(normJ) - logj1 ;
532
533 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
534 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
535 ( double ) std::abs( a ), mLabel.c_str() ) ;
536
537 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
538
539 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
540
541 return rhoFF ;
542 }
543
548 inline real
549 Material::rho_piecewise( const real normJ, const real normB, const real angleNxB ) const
550 {
551 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
552 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
553
554
555 BELFEM_ASSERT(
556 mJcFunction->depends_on( material::JcParameter::normB ) &&
557 mJcFunction->depends_on( material::JcParameter::angleNxB ) &&
558 ! mJcFunction->depends_on( material::JcParameter::T ), "wrong powerlaw for material %s", mLabel.c_str() );
559
560 real rhon = this->rho( gTbulk ) ;
561
562 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
563 {
564 return rhon ;
565 }
566
567
568 real ec = this->constant_property( MaterialProperty::ec ) ;
569 real jc = mJcFunction->eval( normB, angleNxB ) ;
570 real n = mNFunction->eval( normB, angleNxB ) ;
571
572 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
573 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
574
575 //Power law resistivity
576 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
577
578 //Power-law regime
579 real j1 = jc * std::pow(10.0,2.5/n) ;
580 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
581
582 if (normJ <= j1)
583 {
584 return rhoPL ;
585 }
586
587 //Normal regime
588 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
589 real rho3 = rhon ;
590
591 if ( normJ > j3 )
592 {
593 return rhon ;
594 }
595
596 //Flux-flow regime
597 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
598 real rho2 = rhon ;
599
600 //Bezier control points
601 real logj1 = std::log10(j1) ;
602 real logj2 = std::log10(j2) ;
603
604 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
605 real b = logj1 - logj2 ;
606 real c = std::log10(normJ) - logj1 ;
607
608 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
609 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
610 ( double ) std::abs( a ), mLabel.c_str() ) ;
611
612 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
613
614 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
615
616 return rhoFF ;
617 }
618
623 inline real
624 Material::rho_piecewise( const real normJ, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
625 {
626 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
627 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
628
629
630 BELFEM_ASSERT(
631 mJcFunction->depends_on( material::JcParameter::normB ) &&
632 mJcFunction->depends_on( material::JcParameter::angleNxB ) &&
633 ! mJcFunction->depends_on( material::JcParameter::T ), "wrong powerlaw for material %s", mLabel.c_str() );
634
635 real rhon = this->rho( gTbulk ) ;
636
637 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
638 {
639 return rhon ;
640 }
641
642
643 real ec = this->constant_property( MaterialProperty::ec ) ;
644 real jc = mJcFunction->eval( normB, angleNxB )*((this->mDefectFunction) (x,y,z,t)) ; ;
645 real n = mNFunction->eval( normB, angleNxB ) ;
646
647 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
648 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
649
650 //Power law resistivity
651 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
652
653 //Power-law regime
654 real j1 = jc * std::pow(10.0,2.5/n) ;
655 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
656
657 if (normJ <= j1)
658 {
659 return rhoPL ;
660 }
661
662 //Normal regime
663 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
664 real rho3 = rhon ;
665
666 if ( normJ > j3 )
667 {
668 return rhon ;
669 }
670
671 //Flux-flow regime
672 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
673 real rho2 = rhon ;
674
675 //Bezier control points
676 real logj1 = std::log10(j1) ;
677 real logj2 = std::log10(j2) ;
678
679 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
680 real b = logj1 - logj2 ;
681 real c = std::log10(normJ) - logj1 ;
682
683 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
684 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
685 ( double ) std::abs( a ), mLabel.c_str() ) ;
686
687 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
688
689 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
690
691 return rhoFF ;
692 }
693
700 inline real
701 Material::rho_piecewise( const real normJ, const real T, const real normB, const real angleNxB ) const
702 {
703 real rhon = this->rho( T ) ;
704
705 if ( T > this->constant_property( MaterialProperty::T_crit ) )
706 {
707 return rhon ;
708 }
709
710 // dependency-routed evaluation with constants-fallback ( O1 policy )
711 real ec = this->constant_property( MaterialProperty::ec ) ;
712 real jc = this->jc_eval( T, normB, angleNxB ) ;
713 real n = this->n_eval( T, normB, angleNxB ) ;
714
715 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
716 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
717
718 //Power law resistivity
719 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin);
720
721 //Power-law regime
722 real j1 = jc * std::pow(10.0,2.5/n) ;
723 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
724
725 if (normJ <= j1)
726 {
727 return rhoPL ;
728 }
729
730 //Normal regime
731 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
732 real rho3 = rhon ;
733
734 if ( normJ > j3 )
735 {
736 return rhon ;
737 }
738
739 //Flux-flow regime
740 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
741 real rho2 = rhon ;
742
743 //Bezier control points
744 real logj1 = std::log10(j1) ;
745 real logj2 = std::log10(j2) ;
746
747 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
748 real b = logj1 - logj2 ;
749 real c = std::log10(normJ) - logj1 ;
750
751 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
752 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
753 ( double ) std::abs( a ), mLabel.c_str() ) ;
754
755 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
756
757 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
758
759 return rhoFF ;
760 }
761
767 inline real
768 Material::rho_piecewise( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
769 {
770 real rhon = this->rho( T ) ;
771
772 if ( T > this->constant_property( MaterialProperty::T_crit ) ) // (assumed critical temperature)
773 {
774 return rhon ;
775 }
776
777 // dependency-routed evaluation with constants-fallback ( O1 policy )
778 real ec = this->constant_property( MaterialProperty::ec ) ;
779 real jc = this->jc_eval( T, normB, angleNxB ) * ((this->mDefectFunction) (x,y,z,t)) ;
780 real n = this->n_eval( T, normB, angleNxB ) ;
781
782 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
783 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
784
785 //Power law resistivity
786 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
787
788 //Power-law regime
789 real j1 = jc * std::pow(10.0,2.5/n) ;
790 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
791
792 if (normJ <= j1)
793 {
794 return rhoPL ;
795 }
796
797 //Normal regime
798 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
799 real rho3 = rhon ;
800
801 if ( normJ > j3 )
802 {
803 return rhon ;
804 }
805
806 //Flux-flow regime
807 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
808 real rho2 = rhon ;
809
810 //Bezier control points
811 real logj1 = std::log10(j1) ;
812 real logj2 = std::log10(j2) ;
813
814 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
815 real b = logj1 - logj2 ;
816 real c = std::log10(normJ) - logj1 ;
817
818 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
819 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
820 ( double ) std::abs( a ), mLabel.c_str() ) ;
821
822 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
823
824 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
825
826 return rhoFF ;
827 }
828
834 inline real
835 Material::rho_piecewise( const real normJ, const real T ) const
836 {
837 real rhon = this->rho( T ) ;
838
839 if ( T > this->constant_property( MaterialProperty::T_crit ) )
840 {
841 return rhon ;
842 }
843
844 real ec = this->constant_property( MaterialProperty::ec ) ;
845 real jc = this->jc_custom( T ) ;
846 real n = this->n_custom( T ) ;
847
848 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
849 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
850
851 //Power law resistivity
852 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
853
854 //Power-law regime
855 real j1 = jc * std::pow(10.0,2.5/n) ;
856 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
857
858 if (normJ <= j1)
859 {
860 return rhoPL ;
861 }
862
863 //Normal regime
864 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
865 real rho3 = rhon ;
866
867 if ( normJ > j3 )
868 {
869 return rhon ;
870 }
871
872 //Flux-flow regime
873 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
874 real rho2 = rhon ;
875
876 //Bezier control points
877 real logj1 = std::log10(j1) ;
878 real logj2 = std::log10(j2) ;
879
880 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
881 real b = logj1 - logj2 ;
882 real c = std::log10(normJ) - logj1 ;
883
884 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
885 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
886 ( double ) std::abs( a ), mLabel.c_str() ) ;
887
888 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
889
890 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
891
892 return rhoFF ;
893 }
894
899 inline real
900 Material::rho_piecewise( const real normJ, const real T, const real x, const real y, const real z, const real t ) const
901 {
902 real rhon = this->rho( T ) ;
903
904 if ( T > this->constant_property( MaterialProperty::T_crit ) ) // (assumed critical temperature)
905 {
906 return rhon ;
907 }
908
909 real ec = this->constant_property( MaterialProperty::ec ) ;
910 real jc = this->jc_custom( T ) * ((this->mDefectFunction) (x,y,z,t)) ;
911 real n = this->n_custom( T ) ;
912
913 BELFEM_ASSERT( jc > 0.0, "Piecewise power law requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
914 BELFEM_ASSERT( n > 1.0, "Piecewise power law requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
915
916 //Power law resistivity
917 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
918
919 //Power-law regime
920 real j1 = jc * std::pow(10.0,2.5/n) ;
921 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
922
923 if (normJ <= j1)
924 {
925 return rhoPL ;
926 }
927
928 //Normal regime
929 real j3 = j1 * std::pow((rhon/rho1),1.0/(mNff)) ;
930 real rho3 = rhon ;
931
932 if ( normJ > j3 )
933 {
934 return rhon ;
935 }
936
937 //Flux-flow regime
938 real j2 = j1 * std::pow((rhon/rho1),1.0/(n-1.0)) ;
939 real rho2 = rhon ;
940
941 //Bezier control points
942 real logj1 = std::log10(j1) ;
943 real logj2 = std::log10(j2) ;
944
945 real a = logj1 - 2.0*logj2 + std::log10(j3) ;
946 real b = logj1 - logj2 ;
947 real c = std::log10(normJ) - logj1 ;
948
949 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
950 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
951 ( double ) std::abs( a ), mLabel.c_str() ) ;
952
953 real tParam = (b+std::pow(b*b+a*c,0.5))/a ;
954
955 real rhoFF = std::pow(10.0,std::log10(rho1)*(1.0-tParam)*(1.0-tParam) + std::log10(rho2)*2.0*(1.0-tParam)*tParam + std::log10(rho3)*tParam*tParam );
956
957 return rhoFF ;
958 }
959
960//------------------------------------------------------------------------------
961// Function-backed Jc and n accessors
962//------------------------------------------------------------------------------
963
975 inline real
976 Material::n( const real normB, const real angleNxB, const real T ) const
977 {
978 BELFEM_ASSERT( mNFunction != nullptr,
979 "Material %s does not have an N function", mLabel.c_str() ) ;
980
981 return mNFunction->eval( normB, angleNxB, T ) ;
982 }
983
988 inline real
989 Material::n( const real normB, const real angleNxB ) const
990 {
991 BELFEM_ASSERT( mNFunction != nullptr,
992 "Material %s does not have an N function", mLabel.c_str() ) ;
993
994 return mNFunction->eval( normB, angleNxB ) ;
995 }
996
1011 inline real
1012 Material::jc( const real B, const real angleNxB, const real T ) const
1013 {
1014 // Constants-fallback: when jc was set as a plain constant via the
1015 // MaterialFactory constants path, mJcFunction is null but
1016 // constant_property(jc) is valid. See rho_powerlaw(...) for the
1017 // matching pattern in the assembly path.
1018 return ( mJcFunction != nullptr )
1019 ? mJcFunction->eval( B, angleNxB, T )
1020 : this->constant_property( MaterialProperty::jc ) ;
1021 }
1022
1027 inline real
1028 Material::jc( const real B, const real angleNxB ) const
1029 {
1030 // Constants-fallback (see jc(B,angle,T)): have(jc) is true even for a
1031 // plain-constant jc, where mJcFunction is null.
1032 return ( mJcFunction != nullptr )
1033 ? mJcFunction->eval( B, angleNxB )
1034 : this->constant_property( MaterialProperty::jc ) ;
1035 }
1036
1042 inline real
1043 Material::jc( const real B, const real angleNxB, const real T, const real x, const real y, const real z, const real t ) const
1044 {
1045 BELFEM_ASSERT( mJcFunction != nullptr,
1046 "Material %s does not have a Jc function", mLabel.c_str() ) ;
1047
1048 return mJcFunction->eval( B, angleNxB, T )*((this->mDefectFunction) (x,y,z,t)) ;
1049 }
1050
1056 inline real
1057 Material::jc( const real B, const real angleNxB, const real x, const real y, const real z, const real t ) const
1058 {
1059 BELFEM_ASSERT( mJcFunction != nullptr,
1060 "Material %s does not have a Jc function", mLabel.c_str() ) ;
1061
1062 return mJcFunction->eval( B, angleNxB )*((this->mDefectFunction) (x,y,z,t)) ;
1063 }
1064
1065
1066//------------------------------------------------------------------------------
1067// Power-law resistivity derivatives for Newton-Raphson
1068//------------------------------------------------------------------------------
1069
1092 inline real
1093 Material::drho_powerlaw_dJ( const real normJ ) const
1094 {
1095 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
1096 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1097
1098 real ec = this->constant_property( MaterialProperty::ec ) ;
1099 real jc = this->constant_property( MaterialProperty::jc ) ;
1100 real n = this->constant_property( MaterialProperty::n ) ;
1101
1102 if (normJ < BELFEM_EPSILON) return 0.0;
1103
1104 real rhon = this->rho( gTbulk ) ;
1105 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
1106
1107 // dρ_PL/dJ = (Ec/Jc^n)·(n-1)·J^(n-2)
1108 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1109
1110 // dρ/dJ = (ρ²/ρ_PL²)·dρ_PL/dJ
1111 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1112 }
1113
1118 inline real
1119 Material::drho_powerlaw_dJ( const real normJ, const real x, const real y, const real z, const real t ) const
1120 {
1121 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
1122 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1123
1124 real ec = this->constant_property( MaterialProperty::ec ) ;
1125 real jc = this->constant_property( MaterialProperty::jc ) * ((this->mDefectFunction) (x,y,z,t)) ;
1126 real n = this->constant_property( MaterialProperty::n ) ;
1127
1128 if (normJ < BELFEM_EPSILON) return 0.0;
1129
1130 real rhon = this->rho( gTbulk ) ;
1131 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
1132
1133 // dρ_PL/dJ = (Ec/Jc^n)·(n-1)·J^(n-2)
1134 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1135
1136 // dρ/dJ = (ρ²/ρ_PL²)·dρ_PL/dJ
1137 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1138 }
1139
1145 inline real
1146 Material::drho_powerlaw_dJ( const real normJ, const real T ) const
1147 {
1148 real ec = this->constant_property( MaterialProperty::ec ) ;
1149 real jc = this->jc_custom( T ) ;
1150 real n = this->n_custom( T ) ;
1151
1152 if (normJ < BELFEM_EPSILON) return 0.0;
1153
1154 real rhon = this->rho( T ) ;
1155 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
1156
1157 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1158
1159 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1160 }
1161
1166 inline real
1167 Material::drho_powerlaw_dJ( const real normJ, const real T, const real x, const real y, const real z, const real t ) const
1168 {
1169 real ec = this->constant_property( MaterialProperty::ec ) ;
1170 real jc = this->jc_custom( T ) * ((this->mDefectFunction) (x,y,z,t)) ;
1171 real n = this->n_custom( T ) ;
1172
1173 if (normJ < BELFEM_EPSILON) return 0.0;
1174
1175 real rhon = this->rho( T ) ;
1176 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1177
1178 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1179
1180 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1181 }
1182
1187 inline real
1188 Material::drho_powerlaw_dJ( const real normJ, const real normB, const real angleNxB ) const
1189 {
1190 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
1191 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1192 BELFEM_ASSERT( ! mJcFunction->depends_on( material::JcParameter::T ),
1193 "wrong powerlaw derivative for material %s", mLabel.c_str() ) ;
1194
1195 real ec = this->constant_property( MaterialProperty::ec ) ;
1196 real jc = mJcFunction->eval( normB, angleNxB ) ;
1197 real n = mNFunction->eval( normB, angleNxB ) ;
1198
1199 if (normJ < BELFEM_EPSILON) return 0.0;
1200
1201 real rhon = this->rho( gTbulk ) ;
1202 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1203
1204 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1205
1206 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1207 }
1208
1214 inline real
1215 Material::drho_powerlaw_dJ( const real normJ, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
1216 {
1217 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
1218 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1219
1220 real ec = this->constant_property( MaterialProperty::ec ) ;
1221 real jc = mJcFunction->eval( normB, angleNxB ) * ((this->mDefectFunction) (x,y,z,t)) ;
1222 real n = mNFunction->eval( normB, angleNxB ) ;
1223
1224 if (normJ < BELFEM_EPSILON) return 0.0;
1225
1226 real rhon = this->rho( gTbulk ) ;
1227 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1228
1229 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1230
1231 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1232 }
1233
1239 inline real
1240 Material::drho_powerlaw_dJ( const real normJ, const real T, const real normB, const real angleNxB ) const
1241 {
1242 // dependency-routed evaluation with constants-fallback ( O1 policy )
1243 real ec = this->constant_property( MaterialProperty::ec ) ;
1244 real jc = this->jc_eval( T, normB, angleNxB ) ;
1245 real n = this->n_eval( T, normB, angleNxB ) ;
1246
1247 if (normJ < BELFEM_EPSILON) return 0.0;
1248
1249 real rhon = this->rho( T ) ;
1250 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1251
1252 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1253
1254 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1255 }
1256
1262 inline real
1263 Material::drho_powerlaw_dJ( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
1264 {
1265 // dependency-routed evaluation with constants-fallback ( O1 policy )
1266 real ec = this->constant_property( MaterialProperty::ec ) ;
1267 real jc = this->jc_eval( T, normB, angleNxB ) * ((this->mDefectFunction) (x,y,z,t)) ;
1268 real n = this->n_eval( T, normB, angleNxB ) ;
1269
1270 if (normJ < BELFEM_EPSILON) return 0.0;
1271
1272 real rhon = this->rho( T ) ;
1273 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
1274
1275 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1276
1277 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1278 }
1279
1280//------------------------------------------------------------------------------
1281// Piecewise resistivity derivatives for Newton-Raphson
1282//------------------------------------------------------------------------------
1283
1322 inline real
1323 Material::drho_piecewise_dJ( const real normJ ) const
1324 {
1325 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
1326 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1327
1328 real rhon = this->rho( gTbulk ) ;
1329
1330 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
1331 {
1332 return 0.0 ;
1333 }
1334
1335 real ec = this->constant_property( MaterialProperty::ec ) ;
1336 real jc = this->constant_property( MaterialProperty::jc ) ;
1337 real n = this->constant_property( MaterialProperty::n ) ;
1338
1339 if (normJ < BELFEM_EPSILON) return 0.0;
1340
1341 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1342 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1343
1344 // Power-law regime thresholds
1345 real j1 = jc * std::pow(10.0, 2.5/n) ;
1346 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1347
1348 if (normJ <= j1)
1349 {
1350 // Power-law region: same derivative as power-law
1351 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1352 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1353 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1354 }
1355
1356 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1357
1358 if ( normJ > j3 )
1359 {
1360 // Normal regime: constant resistivity
1361 return 0.0 ;
1362 }
1363
1364 // Flux-flow regime: derivative of Bezier curve in log-space
1365 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1366 real rho2 = rhon ;
1367
1368 real logj1 = std::log10(j1) ;
1369 real logj2 = std::log10(j2) ;
1370 real logj3 = std::log10(j3) ;
1371
1372 real a = logj1 - 2.0*logj2 + logj3 ;
1373 real b = logj1 - logj2 ;
1374 real c = std::log10(normJ) - logj1 ;
1375
1376 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1377 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1378 ( double ) std::abs( a ), mLabel.c_str() ) ;
1379
1380 // t parameter from quadratic formula
1381 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1382
1383 // Derivative dt/dJ
1384 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1385 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1386 real dt_dJ = dt_dc * dc_dJ ;
1387
1388 // Bezier curve: log(ρ) = (1-t)² log(rho1) + 2(1-t)t log(rho2) + t² log(rho3)
1389 // d(log(ρ))/dt = -2(1-t) log(rho1) + 2(1-2t) log(rho2) + 2t log(rho3)
1390 real logrho1 = std::log10(rho1) ;
1391 real logrho2 = std::log10(rho2) ;
1392 real logrho3 = std::log10(rhon) ;
1393
1394 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1395
1396 // ρ = 10^(log(ρ))
1397 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1398
1399 // dρ/dJ = ρ · ln(10) · d(log(ρ))/dt · dt/dJ
1400 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1401 }
1402
1407 inline real
1408 Material::drho_piecewise_dJ( const real normJ, const real x, const real y, const real z, const real t ) const
1409 {
1410 BELFEM_ASSERT( this->is_constant( MaterialProperty::jc ), "Material %s does not have a constant jc parameter defined", mLabel.c_str() ) ;
1411 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1412
1413 real rhon = this->rho( gTbulk ) ;
1414
1415 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
1416 {
1417 return 0.0 ;
1418 }
1419
1420 real ec = this->constant_property( MaterialProperty::ec ) ;
1421 real jc = this->constant_property( MaterialProperty::jc ) * ((this->mDefectFunction) (x,y,z,t)) ;
1422 real n = this->constant_property( MaterialProperty::n ) ;
1423
1424 if (normJ < BELFEM_EPSILON) return 0.0;
1425
1426 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1427 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1428
1429 // Power-law regime thresholds
1430 real j1 = jc * std::pow(10.0, 2.5/n) ;
1431 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1432
1433 if (normJ <= j1)
1434 {
1435 // Power-law region: same derivative as power-law
1436 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1437 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1438 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1439 }
1440
1441 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1442
1443 if ( normJ > j3 )
1444 {
1445 // Normal regime: constant resistivity
1446 return 0.0 ;
1447 }
1448
1449 // Flux-flow regime: derivative of Bezier curve in log-space
1450 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1451 real rho2 = rhon ;
1452
1453 real logj1 = std::log10(j1) ;
1454 real logj2 = std::log10(j2) ;
1455 real logj3 = std::log10(j3) ;
1456
1457 real a = logj1 - 2.0*logj2 + logj3 ;
1458 real b = logj1 - logj2 ;
1459 real c = std::log10(normJ) - logj1 ;
1460
1461 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1462 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1463 ( double ) std::abs( a ), mLabel.c_str() ) ;
1464
1465 // t parameter from quadratic formula
1466 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1467
1468 // Derivative dt/dJ
1469 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1470 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1471 real dt_dJ = dt_dc * dc_dJ ;
1472
1473 // Bezier curve: log(ρ) = (1-t)² log(rho1) + 2(1-t)t log(rho2) + t² log(rho3)
1474 // d(log(ρ))/dt = -2(1-t) log(rho1) + 2(1-2t) log(rho2) + 2t log(rho3)
1475 real logrho1 = std::log10(rho1) ;
1476 real logrho2 = std::log10(rho2) ;
1477 real logrho3 = std::log10(rhon) ;
1478
1479 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1480
1481 // ρ = 10^(log(ρ))
1482 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1483
1484 // dρ/dJ = ρ · ln(10) · d(log(ρ))/dt · dt/dJ
1485 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1486 }
1487
1493 inline real
1494 Material::drho_piecewise_dJ( const real normJ, const real T ) const
1495 {
1496 real rhon = this->rho( T ) ;
1497
1498 if ( T > this->constant_property( MaterialProperty::T_crit ) )
1499 {
1500 return 0.0 ;
1501 }
1502
1503 real ec = this->constant_property( MaterialProperty::ec ) ;
1504 real jc = this->jc_custom( T ) ;
1505 real n = this->n_custom( T ) ;
1506
1507 if (normJ < BELFEM_EPSILON) return 0.0;
1508
1509 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1510 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1511
1512 real j1 = jc * std::pow(10.0, 2.5/n) ;
1513 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1514
1515 if (normJ <= j1)
1516 {
1517 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1518 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1519 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1520 }
1521
1522 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1523
1524 if ( normJ > j3 )
1525 {
1526 return 0.0 ;
1527 }
1528
1529 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1530
1531 real logj1 = std::log10(j1) ;
1532 real logj2 = std::log10(j2) ;
1533 real logj3 = std::log10(j3) ;
1534
1535 real a = logj1 - 2.0*logj2 + logj3 ;
1536 real b = logj1 - logj2 ;
1537 real c = std::log10(normJ) - logj1 ;
1538
1539 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1540 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1541 ( double ) std::abs( a ), mLabel.c_str() ) ;
1542
1543 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1544
1545 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1546 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1547 real dt_dJ = dt_dc * dc_dJ ;
1548
1549 real logrho1 = std::log10(rho1) ;
1550 real logrho2 = std::log10(rhon) ;
1551 real logrho3 = std::log10(rhon) ;
1552
1553 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1554
1555 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1556
1557 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1558 }
1559
1564 inline real
1565 Material::drho_piecewise_dJ( const real normJ, const real T, const real x, const real y, const real z, const real t ) const
1566 {
1567 real rhon = this->rho( T ) ;
1568
1569 if ( T > this->constant_property( MaterialProperty::T_crit ) )
1570 {
1571 return 0.0 ;
1572 }
1573
1574 real ec = this->constant_property( MaterialProperty::ec ) ;
1575 real jc = this->jc_custom( T ) * ((this->mDefectFunction) (x,y,z,t)) ;
1576 real n = this->n_custom( T ) ;
1577
1578 if (normJ < BELFEM_EPSILON) return 0.0;
1579
1580 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1581 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1582
1583 real j1 = jc * std::pow(10.0, 2.5/n) ;
1584 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1585
1586 if (normJ <= j1)
1587 {
1588 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1589 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1590 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1591 }
1592
1593 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1594
1595 if ( normJ > j3 )
1596 {
1597 return 0.0 ;
1598 }
1599
1600 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1601
1602 real logj1 = std::log10(j1) ;
1603 real logj2 = std::log10(j2) ;
1604 real logj3 = std::log10(j3) ;
1605
1606 real a = logj1 - 2.0*logj2 + logj3 ;
1607 real b = logj1 - logj2 ;
1608 real c = std::log10(normJ) - logj1 ;
1609
1610 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1611 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1612 ( double ) std::abs( a ), mLabel.c_str() ) ;
1613
1614 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1615
1616 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1617 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1618 real dt_dJ = dt_dc * dc_dJ ;
1619
1620 real logrho1 = std::log10(rho1) ;
1621 real logrho2 = std::log10(rhon) ;
1622 real logrho3 = std::log10(rhon) ;
1623
1624 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1625
1626 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1627
1628 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1629 }
1630
1635 inline real
1636 Material::drho_piecewise_dJ( const real normJ, const real normB, const real angleNxB ) const
1637 {
1638 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
1639 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1640 BELFEM_ASSERT( ! mJcFunction->depends_on( material::JcParameter::T ),
1641 "wrong powerlaw derivative for material %s", mLabel.c_str() ) ;
1642
1643 real rhon = this->rho( gTbulk ) ;
1644
1645 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
1646 {
1647 return 0.0 ;
1648 }
1649
1650 real ec = this->constant_property( MaterialProperty::ec ) ;
1651 real jc = mJcFunction->eval( normB, angleNxB ) ;
1652 real n = mNFunction->eval( normB, angleNxB ) ;
1653
1654 if (normJ < BELFEM_EPSILON) return 0.0;
1655
1656 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1657 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1658
1659 real j1 = jc * std::pow(10.0, 2.5/n) ;
1660 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1661
1662 if (normJ <= j1)
1663 {
1664 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1665 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1666 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1667 }
1668
1669 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1670
1671 if ( normJ > j3 )
1672 {
1673 return 0.0 ;
1674 }
1675
1676 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1677
1678 real logj1 = std::log10(j1) ;
1679 real logj2 = std::log10(j2) ;
1680 real logj3 = std::log10(j3) ;
1681
1682 real a = logj1 - 2.0*logj2 + logj3 ;
1683 real b = logj1 - logj2 ;
1684 real c = std::log10(normJ) - logj1 ;
1685
1686 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1687 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1688 ( double ) std::abs( a ), mLabel.c_str() ) ;
1689
1690 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1691
1692 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1693 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1694 real dt_dJ = dt_dc * dc_dJ ;
1695
1696 real logrho1 = std::log10(rho1) ;
1697 real logrho2 = std::log10(rhon) ;
1698 real logrho3 = std::log10(rhon) ;
1699
1700 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1701
1702 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1703
1704 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1705 }
1706
1712 inline real
1713 Material::drho_piecewise_dJ( const real normJ, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
1714 {
1715 BELFEM_ASSERT( mJcFunction != nullptr, "Material %s does not have Jc function", mLabel.c_str() ) ;
1716 BELFEM_ASSERT( this->is_constant( MaterialProperty::n ), "Material %s does not have a constant n parameter defined", mLabel.c_str() ) ;
1717
1718 real rhon = this->rho( gTbulk ) ;
1719
1720 if ( gTbulk > this->constant_property( MaterialProperty::T_crit ) )
1721 {
1722 return 0.0 ;
1723 }
1724
1725 real ec = this->constant_property( MaterialProperty::ec ) ;
1726 real jc = mJcFunction->eval( normB, angleNxB ) * ((this->mDefectFunction) (x,y,z,t)) ;
1727 real n = mNFunction->eval( normB, angleNxB ) ;
1728
1729 if (normJ < BELFEM_EPSILON) return 0.0;
1730
1731 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1732 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1733
1734 real j1 = jc * std::pow(10.0, 2.5/n) ;
1735 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1736
1737 if (normJ <= j1)
1738 {
1739 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin );
1740 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1741 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1742 }
1743
1744 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1745
1746 if ( normJ > j3 )
1747 {
1748 return 0.0 ;
1749 }
1750
1751 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1752
1753 real logj1 = std::log10(j1) ;
1754 real logj2 = std::log10(j2) ;
1755 real logj3 = std::log10(j3) ;
1756
1757 real a = logj1 - 2.0*logj2 + logj3 ;
1758 real b = logj1 - logj2 ;
1759 real c = std::log10(normJ) - logj1 ;
1760
1761 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1762 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1763 ( double ) std::abs( a ), mLabel.c_str() ) ;
1764
1765 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1766
1767 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1768 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1769 real dt_dJ = dt_dc * dc_dJ ;
1770
1771 real logrho1 = std::log10(rho1) ;
1772 real logrho2 = std::log10(rhon) ;
1773 real logrho3 = std::log10(rhon) ;
1774
1775 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1776
1777 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1778
1779 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1780 }
1781
1787 inline real
1788 Material::drho_piecewise_dJ( const real normJ, const real T, const real normB, const real angleNxB ) const
1789 {
1790 real rhon = this->rho( T ) ;
1791
1792 if ( T > this->constant_property( MaterialProperty::T_crit ) )
1793 {
1794 return 0.0 ;
1795 }
1796
1797 // dependency-routed evaluation with constants-fallback ( O1 policy )
1798 real ec = this->constant_property( MaterialProperty::ec ) ;
1799 real jc = this->jc_eval( T, normB, angleNxB ) ;
1800 real n = this->n_eval( T, normB, angleNxB ) ;
1801
1802 if (normJ < BELFEM_EPSILON) return 0.0;
1803
1804 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1805 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1806
1807 real j1 = jc * std::pow(10.0, 2.5/n) ;
1808 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1809
1810 if (normJ <= j1)
1811 {
1812 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1813 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1814 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1815 }
1816
1817 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1818
1819 if ( normJ > j3 )
1820 {
1821 return 0.0 ;
1822 }
1823
1824 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1825
1826 real logj1 = std::log10(j1) ;
1827 real logj2 = std::log10(j2) ;
1828 real logj3 = std::log10(j3) ;
1829
1830 real a = logj1 - 2.0*logj2 + logj3 ;
1831 real b = logj1 - logj2 ;
1832 real c = std::log10(normJ) - logj1 ;
1833
1834 BELFEM_ERROR( std::abs( a ) > BELFEM_EPSILON,
1835 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1836 ( double ) std::abs( a ), mLabel.c_str() ) ;
1837
1838 cplx d = b*b + a*c ;
1839 BELFEM_ERROR( std::abs( d ) > BELFEM_EPSILON,
1840 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1841 ( double ) std::abs( d ), mLabel.c_str() ) ;
1842
1843 cplx tParam = ( b + std::sqrt( d ) ) / a ;
1844
1845 cplx dt_dc = 0.5 / std::sqrt( d ) ;
1846 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1847 cplx dt_dJ = dt_dc * dc_dJ ;
1848
1849 real logrho1 = std::log10(rho1) ;
1850 real logrho2 = std::log10(rhon) ;
1851 real logrho3 = std::log10(rhon) ;
1852
1853 cplx dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1854
1855 cplx rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1856
1857 cplx drho_dJ = rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1858
1859 BELFEM_ERROR( std::imag( drho_dJ ) < BELFEM_EPSILON,
1860 "Piecewise power law derivative: imaginary part of drho_dJ is non-zero (|a| = %g) for material %s, %g + i * %g",
1861 ( double ) std::imag( drho_dJ ), mLabel.c_str(), std::real( drho_dJ ) , std::imag( drho_dJ ) ) ;
1862
1863 return std::real( drho_dJ ) ;
1864 }
1865
1871 inline real
1872 Material::drho_piecewise_dJ( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
1873 {
1874 real rhon = this->rho( T ) ;
1875
1876 if ( T > this->constant_property( MaterialProperty::T_crit ) )
1877 {
1878 return 0.0 ;
1879 }
1880
1881 // dependency-routed evaluation with constants-fallback ( O1 policy )
1882 real ec = this->constant_property( MaterialProperty::ec ) ;
1883 real jc = this->jc_eval( T, normB, angleNxB ) * ((this->mDefectFunction) (x,y,z,t)) ;
1884 real n = this->n_eval( T, normB, angleNxB ) ;
1885
1886 if (normJ < BELFEM_EPSILON) return 0.0;
1887
1888 BELFEM_ASSERT( jc > 0.0, "Piecewise power law derivative requires jc > 0 (got jc = %g) for material %s", ( double ) jc, mLabel.c_str() ) ;
1889 BELFEM_ASSERT( n > 1.0, "Piecewise power law derivative requires n > 1 (got n = %g) for material %s", ( double ) n, mLabel.c_str() ) ;
1890
1891 real j1 = jc * std::pow(10.0, 2.5/n) ;
1892 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin) ;
1893
1894 if (normJ <= j1)
1895 {
1896 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1897 real drho_PL_dJ = ( ec / std::pow(jc, 2) ) * ( n - 1.0 ) * std::pow( normJ/jc, n - 2.0 );
1898 return drho_PL_dJ/std::pow((1+rhoPL/rhon),2.0) ;
1899 }
1900
1901 real j3 = j1 * std::pow((rhon/rho1), 1.0/(mNff)) ;
1902
1903 if ( normJ > j3 )
1904 {
1905 return 0.0 ;
1906 }
1907
1908 real j2 = j1 * std::pow((rhon/rho1), 1.0/(n-1.0)) ;
1909
1910 real logj1 = std::log10(j1) ;
1911 real logj2 = std::log10(j2) ;
1912 real logj3 = std::log10(j3) ;
1913
1914 real a = logj1 - 2.0*logj2 + logj3 ;
1915 real b = logj1 - logj2 ;
1916 real c = std::log10(normJ) - logj1 ;
1917
1918 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
1919 "Piecewise power law derivative: degenerate Bezier transition (|a| = %g) for material %s",
1920 ( double ) std::abs( a ), mLabel.c_str() ) ;
1921
1922 real tParam = (b + std::pow(b*b + a*c, 0.5)) / a ;
1923
1924 real dt_dc = 0.5 / std::pow(b*b + a*c, 0.5) ;
1925 real dc_dJ = 1.0 / (normJ * std::log(10.0)) ;
1926 real dt_dJ = dt_dc * dc_dJ ;
1927
1928 real logrho1 = std::log10(rho1) ;
1929 real logrho2 = std::log10(rhon) ;
1930 real logrho3 = std::log10(rhon) ;
1931
1932 real dlogrho_dt = -2.0*(1.0-tParam)*logrho1 + 2.0*(1.0-2.0*tParam)*logrho2 + 2.0*tParam*logrho3 ;
1933
1934 real rho = std::pow(10.0, (1.0-tParam)*(1.0-tParam)*logrho1 + 2.0*(1.0-tParam)*tParam*logrho2 + tParam*tParam*logrho3) ;
1935
1936 return rho * std::log(10.0) * dlogrho_dt * dt_dJ ;
1937 }
1938
1939//------------------------------------------------------------------------------
1940// Field derivatives of the HTS laws ( audited 2026-08-13 )
1941//
1942// jc = jc(T,|B|,θ) and n = n(T,|B|,θ) through the JcFunction hooks; at fixed
1943// J, T, θ the unfloored power law p0 = (ec/jc)·(J/jc)^(n−1) has
1944//
1945// ∂p0/∂jc = −n·p0/jc ∂p0/∂n = p0·ln(J/jc)
1946//
1947// so dp0/d|B| = p0·[ −(n/jc)·djc/d|B| + ln(J/jc)·dn/d|B| ].
1948//
1949// Conventions per the audit round ( tmp/ai_exchange/jc_derivative_plumbing.md,
1950// both voices ):
1951// - drho_powerlaw_dB mirrors drho_powerlaw_dJ exactly: parallel factor
1952// ( 1 + rhoPL/rhon )^-2 with the FLOORED rhoPL, unfloored law in the
1953// numerator, J < eps guard.
1954// - drho_piecewise_dB is consistent with rho_piecewise's OWN residual: the
1955// power-law regime returns the RAW dp0 ( rho_piecewise returns raw rhoPL
1956// there, no parallel combination ); the normal regime and T > T_crit are
1957// exactly 0; the Bézier blend is 0 FOR NOW ( staged — the regime
1958// boundaries also move with jc, that derivative is deferred, and the
1959// tangent has a documented jump at j1 ).
1960// - defect overloads modulate jc AND its derivative by the same D(x,y,z,t).
1961//------------------------------------------------------------------------------
1962
1963 inline real
1964 Material::drho_powerlaw_dB( const real normJ, const real T, const real normB, const real angleNxB ) const
1965 {
1966 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
1967
1968 // dependency-routed evaluation with constants-fallback ( O1 policy )
1969 real ec = this->constant_property( MaterialProperty::ec ) ;
1970 real jc = this->jc_eval( T, normB, angleNxB ) ;
1971 real n = this->n_eval( T, normB, angleNxB ) ;
1972
1973 real djcdB = this->djc_eval_dB( T, normB, angleNxB ) ;
1974 real dndB = this->dn_eval_dB( T, normB, angleNxB ) ;
1975
1976 if ( djcdB == 0.0 && dndB == 0.0 ) return 0.0 ;
1977
1978 real rhon = this->rho( T ) ;
1979 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
1980
1981 // unfloored law differentiated, cf. drho_powerlaw_dJ
1982 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
1983
1984 real dp0dB = p0 * ( -( n / jc ) * djcdB
1985 + std::log( normJ / jc ) * dndB ) ;
1986
1987 return dp0dB / std::pow( ( 1.0 + rhoPL / rhon ), 2.0 ) ;
1988 }
1989
1990 inline real
1991 Material::drho_powerlaw_dB( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
1992 {
1993 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
1994
1995 real ec = this->constant_property( MaterialProperty::ec ) ;
1996
1997 // defect modulation applies to jc AND its derivative ( product rule
1998 // with D independent of |B| ); n is not modulated, cf. rho_powerlaw
1999 real tD = ( this->mDefectFunction )( x, y, z, t ) ;
2000 real jc = this->jc_eval( T, normB, angleNxB ) * tD ;
2001 real n = this->n_eval( T, normB, angleNxB ) ;
2002
2003 real djcdB = this->djc_eval_dB( T, normB, angleNxB ) * tD ;
2004 real dndB = this->dn_eval_dB( T, normB, angleNxB ) ;
2005
2006 if ( djcdB == 0.0 && dndB == 0.0 ) return 0.0 ;
2007
2008 real rhon = this->rho( T ) ;
2009 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
2010
2011 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2012
2013 real dp0dB = p0 * ( -( n / jc ) * djcdB
2014 + std::log( normJ / jc ) * dndB ) ;
2015
2016 return dp0dB / std::pow( ( 1.0 + rhoPL / rhon ), 2.0 ) ;
2017 }
2018
2019 inline real
2020 Material::drho_piecewise_dB( const real normJ, const real T, const real normB, const real angleNxB ) const
2021 {
2022 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2023
2024 // above T_crit rho_piecewise returns rhon, independent of jc and n
2025 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2026 {
2027 return 0.0 ;
2028 }
2029
2030 real ec = this->constant_property( MaterialProperty::ec ) ;
2031 real jc = this->jc_eval( T, normB, angleNxB ) ;
2032 real n = this->n_eval( T, normB, angleNxB ) ;
2033
2034 real djcdB = this->djc_eval_dB( T, normB, angleNxB ) ;
2035 real dndB = this->dn_eval_dB( T, normB, angleNxB ) ;
2036
2037 if ( djcdB == 0.0 && dndB == 0.0 ) return 0.0 ;
2038
2039 // outside the power-law regime: normal regime is exactly 0, the
2040 // Bézier blend is staged 0 ( see the block comment above )
2041 real j1 = jc * std::pow( 10.0, 2.5 / n ) ;
2042 if ( normJ > j1 ) return 0.0 ;
2043
2044 // power-law regime: rho_piecewise returns the RAW power law here,
2045 // so the exact derivative carries NO parallel factor
2046 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2047
2048 return p0 * ( -( n / jc ) * djcdB
2049 + std::log( normJ / jc ) * dndB ) ;
2050 }
2051
2052 inline real
2053 Material::drho_piecewise_dB( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
2054 {
2055 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2056
2057 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2058 {
2059 return 0.0 ;
2060 }
2061
2062 real ec = this->constant_property( MaterialProperty::ec ) ;
2063
2064 real tD = ( this->mDefectFunction )( x, y, z, t ) ;
2065 real jc = this->jc_eval( T, normB, angleNxB ) * tD ;
2066 real n = this->n_eval( T, normB, angleNxB ) ;
2067
2068 real djcdB = this->djc_eval_dB( T, normB, angleNxB ) * tD ;
2069 real dndB = this->dn_eval_dB( T, normB, angleNxB ) ;
2070
2071 if ( djcdB == 0.0 && dndB == 0.0 ) return 0.0 ;
2072
2073 real j1 = jc * std::pow( 10.0, 2.5 / n ) ;
2074 if ( normJ > j1 ) return 0.0 ;
2075
2076 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2077
2078 return p0 * ( -( n / jc ) * djcdB
2079 + std::log( normJ / jc ) * dndB ) ;
2080 }
2081
2082//------------------------------------------------------------------------------
2083// Temperature derivatives of the HTS laws ( T-leg, audited 2026-08-13 )
2084//
2085// The quench-feedback tangent: jc(T), n(T) AND rho_n(T) all move. With
2086// a = rho_n, b = p0 = (ec/jc)·(J/jc)^(n−1) and c = parallel(a,b),
2087//
2088// ∂c/∂a = b²/(a+b)² = ( c/a )² ∂c/∂b = a²/(a+b)² = 1/(1+b/a)²
2089//
2090// dp0/dT = p0·( −(n/jc)·djc/dT + ln(J/jc)·dn/dT ), same chain as the |B|
2091// channel. These closed forms replace the retired b = a·c/(c−a)
2092// reconstruction in the old compute_drhodT_hts, which was sign-flipped
2093// ( correct: a·c/(a−c) ) and whose dadT weighting c²/(a−2c)² diverged at
2094// the flux-flow crossover rho_PL = rho_n. Conventions:
2095// - floored rhoPL in the parallel factors, unfloored law differentiated
2096// ( dB/dJ convention );
2097// - NO ( djc==0 && dn==0 ) early-out: the rho_n term is the only correct
2098// T-dependence of a constant-jc material and must survive;
2099// - drho_piecewise_dT follows rho_piecewise's OWN residual branch for
2100// branch; the Bézier blend carries BOTH parts ( 2026-08-14 ): the
2101// frozen-knot partial ( rho_n control points + the explicit rho1
2102// dependence ) AND the knot motion, i.e. dtParam/dT through j1(T),
2103// j2(T), j3(T). The knot motion dominates — with it staged out the
2104// frozen part alone reproduced only 4-33 % of dρ/dT across the blend,
2105// which is why the thermal Newton stalled in flux-flow. It degrades to
2106// the frozen part alone at the degenerate n−1 == mNff transition ( see
2107// the guard at the discriminant ). Jumps at j1 and T_crit are
2108// inherited from the residual itself;
2109// - defect overloads modulate jc AND djc/dT by the same D(x,y,z,t).
2110//------------------------------------------------------------------------------
2111
2112 inline real
2113 Material::drho_powerlaw_dT( const real normJ, const real T, const real normB, const real angleNxB ) const
2114 {
2115 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2116
2117 // dependency-routed evaluation with constants-fallback ( O1 policy )
2118 real ec = this->constant_property( MaterialProperty::ec ) ;
2119 real jc = this->jc_eval( T, normB, angleNxB ) ;
2120 real n = this->n_eval( T, normB, angleNxB ) ;
2121
2122 real djcdT = this->djc_eval_dT( T, normB, angleNxB ) ;
2123 real dndT = this->dn_eval_dT( T, normB, angleNxB ) ;
2124
2125 real rhon = this->rho( T ) ;
2126 real dadT = this->drhodT( T ) ;
2127
2128 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
2129
2130 // unfloored law differentiated, cf. drho_powerlaw_dB
2131 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2132
2133 real dp0dT = p0 * ( -( n / jc ) * djcdT
2134 + std::log( normJ / jc ) * dndT ) ;
2135
2136 return dp0dT / std::pow( ( 1.0 + rhoPL / rhon ), 2.0 )
2137 + dadT * std::pow( rhoPL / ( rhon + rhoPL ), 2.0 ) ;
2138 }
2139
2140 inline real
2141 Material::drho_powerlaw_dT( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
2142 {
2143 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2144
2145 real ec = this->constant_property( MaterialProperty::ec ) ;
2146
2147 // defect modulation applies to jc AND its derivative ( product rule
2148 // with D independent of T ); n is not modulated, cf. rho_powerlaw
2149 real tD = ( this->mDefectFunction )( x, y, z, t ) ;
2150 real jc = this->jc_eval( T, normB, angleNxB ) * tD ;
2151 real n = this->n_eval( T, normB, angleNxB ) ;
2152
2153 real djcdT = this->djc_eval_dT( T, normB, angleNxB ) * tD ;
2154 real dndT = this->dn_eval_dT( T, normB, angleNxB ) ;
2155
2156 real rhon = this->rho( T ) ;
2157 real dadT = this->drhodT( T ) ;
2158
2159 real rhoPL = std::max(( ec / jc ) * std::pow( normJ / jc, n - 1 ), mRhoMin ) ;
2160
2161 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2162
2163 real dp0dT = p0 * ( -( n / jc ) * djcdT
2164 + std::log( normJ / jc ) * dndT ) ;
2165
2166 return dp0dT / std::pow( ( 1.0 + rhoPL / rhon ), 2.0 )
2167 + dadT * std::pow( rhoPL / ( rhon + rhoPL ), 2.0 ) ;
2168 }
2169
2170 inline real
2171 Material::drho_piecewise_dT( const real normJ, const real T, const real normB, const real angleNxB ) const
2172 {
2173 real dadT = this->drhodT( T ) ;
2174
2175 // above T_crit rho_piecewise returns rhon(T) exactly
2176 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2177 {
2178 return dadT ;
2179 }
2180
2181 // dependency-routed evaluation with constants-fallback ( O1 policy )
2182 real ec = this->constant_property( MaterialProperty::ec ) ;
2183 real jc = this->jc_eval( T, normB, angleNxB ) ;
2184 real n = this->n_eval( T, normB, angleNxB ) ;
2185
2186 real djcdT = this->djc_eval_dT( T, normB, angleNxB ) ;
2187 real dndT = this->dn_eval_dT( T, normB, angleNxB ) ;
2188
2189 // power-law regime: rho_piecewise returns the RAW power law here,
2190 // so the exact derivative carries NO parallel factor and NO rho_n
2191 // term ( J < eps is inside this branch: p0 -> 0 for n > 1 )
2192 real j1 = jc * std::pow( 10.0, 2.5 / n ) ;
2193 if ( normJ <= j1 )
2194 {
2195 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2196
2197 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2198
2199 return p0 * ( -( n / jc ) * djcdT
2200 + std::log( normJ / jc ) * dndT ) ;
2201 }
2202
2203 real rhon = this->rho( T ) ;
2204 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin ) ;
2205
2206 // normal regime: rho_piecewise returns rhon(T) exactly
2207 real j3 = j1 * std::pow( rhon / rho1, 1.0 / ( mNff ) ) ;
2208 if ( normJ > j3 )
2209 {
2210 return dadT ;
2211 }
2212
2213 // Bézier blend: FULL derivative ( 2026-08-14 ). The residual's Bézier
2214 // machinery is reproduced verbatim ( cf. rho_piecewise ) to obtain
2215 // tParam and rhoFF, then differentiated in two parts — (i) the
2216 // control points rho1 and rho_n at frozen knots, and (ii) the knot
2217 // motion, i.e. dtParam/dT through j1(T), j2(T), j3(T). Part (ii)
2218 // dominates: it was staged out when the T-leg first landed, and a
2219 // finite-difference check then showed the frozen part alone
2220 // reproduces only 4-33 % of dρ/dT across the blend. With both parts
2221 // the derivative is exact to roundoff.
2222 real j2 = j1 * std::pow( rhon / rho1, 1.0 / ( n - 1.0 ) ) ;
2223
2224 real logj1 = std::log10( j1 ) ;
2225 real logj2 = std::log10( j2 ) ;
2226
2227 real a = logj1 - 2.0 * logj2 + std::log10( j3 ) ;
2228 real b = logj1 - logj2 ;
2229 real c = std::log10( normJ ) - logj1 ;
2230
2231 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
2232 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
2233 ( double ) std::abs( a ), mLabel.c_str() ) ;
2234
2235 real tDisc = b * b + a * c ;
2236 real tSqrt = std::pow( tDisc, 0.5 ) ;
2237 real tParam = ( b + tSqrt ) / a ;
2238
2239 real rhoFF = std::pow( 10.0,
2240 std::log10( rho1 ) * ( 1.0 - tParam ) * ( 1.0 - tParam )
2241 + std::log10( rhon ) * 2.0 * ( 1.0 - tParam ) * tParam
2242 + std::log10( rhon ) * tParam * tParam ) ;
2243
2244 // (i) control points move: rho1 = (ec/jc)·10^{2.5(n−1)/n} and rho_n(T)
2245 real dlnrho1dT = -djcdT / jc + 2.5 * std::log( 10.0 ) / ( n * n ) * dndT ;
2246
2247 real tFrozen = rhoFF * ( ( 2.0 * tParam - tParam * tParam ) * dadT / rhon
2248 + ( 1.0 - tParam ) * ( 1.0 - tParam ) * dlnrho1dT ) ;
2249
2250 // (ii) KNOT MOTION — the dominant part in this regime, and the reason
2251 // the frozen-knot form alone was not usable: j1, j2, j3 all slide as
2252 // jc(T) falls, so an element at fixed J moves DEEPER into the
2253 // transition and tParam itself carries a T-derivative. Measured
2254 // against a finite difference at 77 K REBCO constants, the frozen
2255 // terms above reproduce only 4-33 % of dρ/dT across the blend; with
2256 // this term the agreement is exact to roundoff.
2257 //
2258 // ln j1 = ln jc + 2.5·ln10/n, hence dln(j1)/dT = −dln(rho1)/dT
2259 // exactly; j3 and j2 follow from their definitions, j2 additionally
2260 // through the explicit 1/(n−1) exponent.
2261 real dlnj1dT = -dlnrho1dT ;
2262 real dlnrhondT = dadT / rhon ;
2263 real dlnj3dT = dlnj1dT + ( dlnrhondT - dlnrho1dT ) / mNff ;
2264 real dlnj2dT = dlnj1dT + ( dlnrhondT - dlnrho1dT ) / ( n - 1.0 )
2265 - dndT / ( ( n - 1.0 ) * ( n - 1.0 ) )
2266 * std::log( rhon / rho1 ) ;
2267
2268 // a, b, c are base-10 logs of the knots: d/dT picks up 1/ln10
2269 real tInvLn10 = 1.0 / std::log( 10.0 ) ;
2270 real dA = ( dlnj1dT - 2.0 * dlnj2dT + dlnj3dT ) * tInvLn10 ;
2271 real dB = ( dlnj1dT - dlnj2dT ) * tInvLn10 ;
2272 real dC = -dlnj1dT * tInvLn10 ;
2273
2274 // tSqrt is the Bézier root's discriminant, and it is NOT bounded away
2275 // from zero for every legal material. With R = log10( rho_n / rho1 ),
2276 // p = n−1 and q = mNff, at the upper knot J = j3 it reduces to
2277 //
2278 // b² + ac = R² ( 1/p − 1/q )²
2279 //
2280 // which VANISHES when p == q — n = 4 at the default mNff = 3, legal
2281 // under the n > 1 precondition, and exactly where j2 == j3, i.e. the
2282 // Bézier control points coincide. The existing |a| > eps assert does
2283 // not cover this ( a = −R/q there, nonzero ).
2284 //
2285 // The singularity is REMOVABLE, not a pole: dt/dT carries a 1/tSqrt
2286 // factor, but the weight derivative it multiplies is ∝ ( 1 − t ), and
2287 // at p == q one has a == b, hence 1 − t = −tSqrt/a. The product is
2288 // finite ( checked numerically: it converges as J → j3 ). What the
2289 // guard avoids is therefore an indeterminate 0/0 in floating point,
2290 // not an infinite physical derivative.
2291 //
2292 // The fallback is exact at the endpoint rather than merely safe: at
2293 // tSqrt = 0 the root gives t = 1, where rhoFF = rho_n and tFrozen
2294 // reduces to dadT — the normal-branch derivative — so the tangent
2295 // stays continuous with the J > j3 branch. In the infinitesimal
2296 // neighborhood where the guard bites, a finite knot contribution is
2297 // dropped; that is a documented approximation for a degenerate
2298 // material ( n = 4 ), not a mechanism any REBCO deck reaches.
2299 // Guarded relative to |b|, the natural scale of the root
2300 // ( tSqrt = |b| at c = 0 ).
2301 // Tested on the DISCRIMINANT, not on its root, and with the
2302 // comparison negated: roundoff can push tDisc slightly negative at
2303 // the degenerate point, which makes tSqrt a NaN — and a NaN fails
2304 // every ordinary comparison, so a `tSqrt <= tol` guard would be
2305 // bypassed and the NaN would propagate into the tangent. `!( x > tol )`
2306 // takes the fallback for NaN as well.
2307 if ( ! ( tDisc > BELFEM_EPSILON * ( b * b + 1.0 ) ) )
2308 {
2309 return tFrozen ;
2310 }
2311
2312 real dSdT = ( 2.0 * b * dB + dA * c + a * dC ) / ( 2.0 * tSqrt ) ;
2313 real dtdT = ( ( dB + dSdT ) - tParam * dA ) / a ;
2314
2315 // d(rhoFF)/dt at frozen control points: the weight derivatives
2316 // −2(1−t) on log ρ1 and (2−2t) on log ρ_n collapse to
2317 // 2(1−t)·ln( ρ_n / ρ1 ) once the ln10 factors cancel
2318 return tFrozen
2319 + rhoFF * 2.0 * ( 1.0 - tParam ) * std::log( rhon / rho1 ) * dtdT ;
2320 }
2321
2322 inline real
2323 Material::drho_piecewise_dT( const real normJ, const real T, const real normB, const real angleNxB, const real x, const real y, const real z, const real t ) const
2324 {
2325 real dadT = this->drhodT( T ) ;
2326
2327 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2328 {
2329 return dadT ;
2330 }
2331
2332 real ec = this->constant_property( MaterialProperty::ec ) ;
2333
2334 // defect modulation applies to jc AND its derivative ( product rule
2335 // with D independent of T ); n is not modulated, cf. rho_piecewise
2336 real tD = ( this->mDefectFunction )( x, y, z, t ) ;
2337 real jc = this->jc_eval( T, normB, angleNxB ) * tD ;
2338 real n = this->n_eval( T, normB, angleNxB ) ;
2339
2340 real djcdT = this->djc_eval_dT( T, normB, angleNxB ) * tD ;
2341 real dndT = this->dn_eval_dT( T, normB, angleNxB ) ;
2342
2343 real j1 = jc * std::pow( 10.0, 2.5 / n ) ;
2344 if ( normJ <= j1 )
2345 {
2346 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2347
2348 real p0 = ( ec / jc ) * std::pow( normJ / jc, n - 1 ) ;
2349
2350 return p0 * ( -( n / jc ) * djcdT
2351 + std::log( normJ / jc ) * dndT ) ;
2352 }
2353
2354 real rhon = this->rho( T ) ;
2355 real rho1 = std::max(( ec / jc ) * std::pow( j1 / jc, n - 1 ), mRhoMin ) ;
2356
2357 real j3 = j1 * std::pow( rhon / rho1, 1.0 / ( mNff ) ) ;
2358 if ( normJ > j3 )
2359 {
2360 return dadT ;
2361 }
2362
2363 real j2 = j1 * std::pow( rhon / rho1, 1.0 / ( n - 1.0 ) ) ;
2364
2365 real logj1 = std::log10( j1 ) ;
2366 real logj2 = std::log10( j2 ) ;
2367
2368 real a = logj1 - 2.0 * logj2 + std::log10( j3 ) ;
2369 real b = logj1 - logj2 ;
2370 real c = std::log10( normJ ) - logj1 ;
2371
2372 BELFEM_ASSERT( std::abs( a ) > BELFEM_EPSILON,
2373 "Piecewise power law: degenerate Bezier transition (|a| = %g) for material %s",
2374 ( double ) std::abs( a ), mLabel.c_str() ) ;
2375
2376 real tDisc = b * b + a * c ;
2377 real tSqrt = std::pow( tDisc, 0.5 ) ;
2378 real tParam = ( b + tSqrt ) / a ;
2379
2380 real rhoFF = std::pow( 10.0,
2381 std::log10( rho1 ) * ( 1.0 - tParam ) * ( 1.0 - tParam )
2382 + std::log10( rhon ) * 2.0 * ( 1.0 - tParam ) * tParam
2383 + std::log10( rhon ) * tParam * tParam ) ;
2384
2385 // (i) control points move: rho1 = (ec/jc)·10^{2.5(n−1)/n} and rho_n(T)
2386 real dlnrho1dT = -djcdT / jc + 2.5 * std::log( 10.0 ) / ( n * n ) * dndT ;
2387
2388 real tFrozen = rhoFF * ( ( 2.0 * tParam - tParam * tParam ) * dadT / rhon
2389 + ( 1.0 - tParam ) * ( 1.0 - tParam ) * dlnrho1dT ) ;
2390
2391 // (ii) KNOT MOTION — the dominant part in this regime, and the reason
2392 // the frozen-knot form alone was not usable: j1, j2, j3 all slide as
2393 // jc(T) falls, so an element at fixed J moves DEEPER into the
2394 // transition and tParam itself carries a T-derivative. Measured
2395 // against a finite difference at 77 K REBCO constants, the frozen
2396 // terms above reproduce only 4-33 % of dρ/dT across the blend; with
2397 // this term the agreement is exact to roundoff.
2398 //
2399 // ln j1 = ln jc + 2.5·ln10/n, hence dln(j1)/dT = −dln(rho1)/dT
2400 // exactly; j3 and j2 follow from their definitions, j2 additionally
2401 // through the explicit 1/(n−1) exponent.
2402 real dlnj1dT = -dlnrho1dT ;
2403 real dlnrhondT = dadT / rhon ;
2404 real dlnj3dT = dlnj1dT + ( dlnrhondT - dlnrho1dT ) / mNff ;
2405 real dlnj2dT = dlnj1dT + ( dlnrhondT - dlnrho1dT ) / ( n - 1.0 )
2406 - dndT / ( ( n - 1.0 ) * ( n - 1.0 ) )
2407 * std::log( rhon / rho1 ) ;
2408
2409 // a, b, c are base-10 logs of the knots: d/dT picks up 1/ln10
2410 real tInvLn10 = 1.0 / std::log( 10.0 ) ;
2411 real dA = ( dlnj1dT - 2.0 * dlnj2dT + dlnj3dT ) * tInvLn10 ;
2412 real dB = ( dlnj1dT - dlnj2dT ) * tInvLn10 ;
2413 real dC = -dlnj1dT * tInvLn10 ;
2414
2415 // tSqrt is the Bézier root's discriminant, and it is NOT bounded away
2416 // from zero for every legal material. With R = log10( rho_n / rho1 ),
2417 // p = n−1 and q = mNff, at the upper knot J = j3 it reduces to
2418 //
2419 // b² + ac = R² ( 1/p − 1/q )²
2420 //
2421 // which VANISHES when p == q — n = 4 at the default mNff = 3, legal
2422 // under the n > 1 precondition, and exactly where j2 == j3, i.e. the
2423 // Bézier control points coincide. The existing |a| > eps assert does
2424 // not cover this ( a = −R/q there, nonzero ).
2425 //
2426 // The singularity is REMOVABLE, not a pole: dt/dT carries a 1/tSqrt
2427 // factor, but the weight derivative it multiplies is ∝ ( 1 − t ), and
2428 // at p == q one has a == b, hence 1 − t = −tSqrt/a. The product is
2429 // finite ( checked numerically: it converges as J → j3 ). What the
2430 // guard avoids is therefore an indeterminate 0/0 in floating point,
2431 // not an infinite physical derivative.
2432 //
2433 // The fallback is exact at the endpoint rather than merely safe: at
2434 // tSqrt = 0 the root gives t = 1, where rhoFF = rho_n and tFrozen
2435 // reduces to dadT — the normal-branch derivative — so the tangent
2436 // stays continuous with the J > j3 branch. In the infinitesimal
2437 // neighborhood where the guard bites, a finite knot contribution is
2438 // dropped; that is a documented approximation for a degenerate
2439 // material ( n = 4 ), not a mechanism any REBCO deck reaches.
2440 // Guarded relative to |b|, the natural scale of the root
2441 // ( tSqrt = |b| at c = 0 ).
2442 // Tested on the DISCRIMINANT, not on its root, and with the
2443 // comparison negated: roundoff can push tDisc slightly negative at
2444 // the degenerate point, which makes tSqrt a NaN — and a NaN fails
2445 // every ordinary comparison, so a `tSqrt <= tol` guard would be
2446 // bypassed and the NaN would propagate into the tangent. `!( x > tol )`
2447 // takes the fallback for NaN as well.
2448 if ( ! ( tDisc > BELFEM_EPSILON * ( b * b + 1.0 ) ) )
2449 {
2450 return tFrozen ;
2451 }
2452
2453 real dSdT = ( 2.0 * b * dB + dA * c + a * dC ) / ( 2.0 * tSqrt ) ;
2454 real dtdT = ( ( dB + dSdT ) - tParam * dA ) / a ;
2455
2456 // d(rhoFF)/dt at frozen control points: the weight derivatives
2457 // −2(1−t) on log ρ1 and (2−2t) on log ρ_n collapse to
2458 // 2(1−t)·ln( ρ_n / ρ1 ) once the ln10 factors cancel
2459 return tFrozen
2460 + rhoFF * 2.0 * ( 1.0 - tParam ) * std::log( rhon / rho1 ) * dtdT ;
2461 }
2462
2463//------------------------------------------------------------------------------
2464// The riva law ( 2026-08-27 )
2465//
2466// The same parallel model as rho_powerlaw — the superconducting power-law
2467// channel in parallel with the normal-state channel ( Duron et al. 2004;
2468// Riva 2021, EPFL thesis 8754, Eq. 5.4 ) — but TOTAL over the full range a
2469// measured jc/n table can produce mid-iterate:
2470// - jc_eff ≤ 0 or nonfinite ( dead defect D(x)=0, spline underflow ):
2471// the superconducting channel is gone, the material is fully normal;
2472// - the power-law channel is evaluated in log10 space; past a cap the
2473// parallel combination is ρn to machine precision, so the residual
2474// returns ρn and the tangents return the matching normal-branch values
2475// instead of the raw inf/inf;
2476// - n arrives from n_eval pre-floored at 1 ( ohmic limit ρPL = ec/jc,
2477// J-independent; dn_eval_* are 0 while the floor binds );
2478// - no n > 1 precondition anywhere. ( A PROVABLY bad n source -- a table
2479// or constant whose bound sits at or below 1 -- is refused at setup in
2480// set_resistivity_law; the runtime stays total for what remains. )
2481// Weight algebra as in the powerlaw T-leg: with w = ρn/(ρPL+ρn),
2482// ∂ρ/∂ρPL = w², ∂ρ/∂ρn = (1−w)².
2483// Deviation from Riva Eq. 5.1/5.4: BELFEM keeps its mRhoMin FLOOR semantics
2484// ( zero by default ) instead of Riva's additive 1e-17 regularization — an
2485// additive term would reopen the 2026-08-10 value/tangent desync.
2486//------------------------------------------------------------------------------
2487
2488 inline bool
2489 Material::riva_rho_pl( const real normJ, const real jc, const real n, const real ec, real & rhoPL ) const
2490 {
2491 if ( n <= 1.0 )
2492 {
2493 // ohmic floor: J-independent resistor. A subnormal jc passes the
2494 // caller's positivity gate but overflows ec/jc, and a bad deck
2495 // can set ec <= 0 — an infinite or negative channel resistance
2496 // takes the fully-normal branch
2497 rhoPL = ec / jc ;
2498 return std::isfinite( rhoPL ) && rhoPL >= 0.0 ;
2499 }
2500 else if ( normJ < BELFEM_EPSILON )
2501 {
2502 rhoPL = mRhoMin ;
2503 }
2504 else
2505 {
2506 real lg = std::log10( ec / jc )
2507 + ( n - 1.0 ) * std::log10( normJ / jc ) ;
2508
2509 // past this cap 1/ρPL vanishes to machine precision against any
2510 // physical ρn — the caller takes the fully-normal branch
2511 // negated NaN-aware comparison: a NaN lg ( NaN n from the
2512 // table, negative ec, or inf·0 at J == jc with infinite n )
2513 // must also land here, and NaN fails every ordinary comparison
2514 if ( ! ( lg <= 250.0 ) )
2515 {
2516 return false ;
2517 }
2518 rhoPL = std::max( std::pow( 10.0, lg ), mRhoMin ) ;
2519 }
2520 return true ;
2521 }
2522
2523 inline real
2524 Material::rho_riva( const real normJ, const real T, const real normB, const real angleNxB ) const
2525 {
2526 real rhon = this->rho( T ) ;
2527
2528 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2529 {
2530 return rhon ;
2531 }
2532
2533 real ec = this->constant_property( MaterialProperty::ec ) ;
2534 real jc = this->jc_eval( T, normB, angleNxB ) ;
2535 real n = this->n_eval( T, normB, angleNxB ) ;
2536
2537 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) )
2538 {
2539 return rhon ;
2540 }
2541
2542 real rhoPL ;
2543 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2544 {
2545 return rhon ;
2546 }
2547
2548 // parallel combination, branch-stable against a huge ρPL
2549 return rhoPL > rhon ? rhon / ( 1.0 + rhon / rhoPL )
2550 : rhoPL / ( 1.0 + rhoPL / rhon ) ;
2551 }
2552
2553 inline real
2554 Material::drho_riva_dJ( const real normJ, const real T, const real normB, const real angleNxB ) const
2555 {
2556 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2557 {
2558 return 0.0 ;
2559 }
2560
2561 real ec = this->constant_property( MaterialProperty::ec ) ;
2562 real jc = this->jc_eval( T, normB, angleNxB ) ;
2563 real n = this->n_eval( T, normB, angleNxB ) ;
2564
2565 // fully-normal fallbacks: ρ = ρn there, which is J-independent
2566 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) ) return 0.0 ;
2567 if ( n <= 1.0 ) return 0.0 ;
2568 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2569
2570 real rhoPL ;
2571 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2572 {
2573 return 0.0 ;
2574 }
2575
2576 real rhon = this->rho( T ) ;
2577 real w = rhon / ( rhoPL + rhon ) ;
2578
2579 // grouped so each factor stays bounded: w·ρPL ≤ ρn even when ρPL
2580 // is huge, and w ≤ 1
2581 return ( w * rhoPL ) * ( w * ( n - 1.0 ) / normJ ) ;
2582 }
2583
2584 inline real
2585 Material::drho_riva_dB( const real normJ, const real T, const real normB, const real angleNxB ) const
2586 {
2587 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2588 {
2589 return 0.0 ;
2590 }
2591
2592 real ec = this->constant_property( MaterialProperty::ec ) ;
2593 real jc = this->jc_eval( T, normB, angleNxB ) ;
2594 real n = this->n_eval( T, normB, angleNxB ) ;
2595
2596 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) ) return 0.0 ;
2597
2598 // the J < eps early-out is only valid for n > 1 ( rhoPL -> mRhoMin );
2599 // at the ohmic floor rhoPL = ec/jc is J-independent but still
2600 // B-dependent through jc, so the tangent must follow the residual
2601 if ( n > 1.0 && normJ < BELFEM_EPSILON ) return 0.0 ;
2602
2603 real djcdB = this->djc_eval_dB( T, normB, angleNxB ) ;
2604 real dndB = this->dn_eval_dB( T, normB, angleNxB ) ;
2605
2606 if ( djcdB == 0.0 && dndB == 0.0 ) return 0.0 ;
2607
2608 real rhoPL ;
2609 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2610 {
2611 return 0.0 ;
2612 }
2613
2614 real rhon = this->rho( T ) ;
2615 real w = rhon / ( rhoPL + rhon ) ;
2616
2617 // dρPL/dB = ρPL·( −(n/jc)·djc/dB + ln(J/jc)·dn/dB ); the ln term
2618 // vanishes with dndB = 0 while the n-floor binds
2619 // djc/jc grouped as a ratio ( bounded for every real jc source ) so
2620 // a tiny jc cannot overflow n/jc on its own; weights grouped so each
2621 // factor stays bounded ( w·ρPL ≤ ρn, w ≤ 1 )
2622 return ( w * rhoPL ) * ( w * ( -( n * ( djcdB / jc ) )
2623 + ( normJ < BELFEM_EPSILON ? 0.0
2624 : std::log( normJ ) - std::log( jc ) ) * dndB ) ) ;
2625 }
2626
2627 inline real
2628 Material::drho_riva_dT( const real normJ, const real T, const real normB, const real angleNxB ) const
2629 {
2630 real dadT = this->drhodT( T ) ;
2631
2632 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2633 {
2634 return dadT ;
2635 }
2636
2637 real ec = this->constant_property( MaterialProperty::ec ) ;
2638 real jc = this->jc_eval( T, normB, angleNxB ) ;
2639 real n = this->n_eval( T, normB, angleNxB ) ;
2640
2641 // fully-normal fallbacks: ρ = ρn there, tangent follows it
2642 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) ) return dadT ;
2643
2644 real rhoPL ;
2645 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2646 {
2647 return dadT ;
2648 }
2649
2650 real djcdT = this->djc_eval_dT( T, normB, angleNxB ) ;
2651 real dndT = this->dn_eval_dT( T, normB, angleNxB ) ;
2652
2653 real rhon = this->rho( T ) ;
2654 real w = rhon / ( rhoPL + rhon ) ;
2655 real v = rhoPL / ( rhoPL + rhon ) ;
2656
2657 // dρPL/dT = ρPL·( −(n/jc)·djc/dT + ln(J/jc)·dn/dT ); the ln term is
2658 // suppressed below BELFEM_EPSILON where ρPL is 0 ( n > 1 ) or
2659 // J-independent ( n = 1, dndT = 0 under the floor )
2660 // djc/jc grouped as a ratio and the weights as (w·ρPL)·(w·dln) so
2661 // no factor can overflow on a tiny-but-finite jc
2662 real dlnp = ( normJ < BELFEM_EPSILON && n > 1.0 ) ? 0.0 :
2663 -( n * ( djcdT / jc ) )
2664 + ( normJ < BELFEM_EPSILON ? 0.0
2665 : std::log( normJ ) - std::log( jc ) ) * dndT ;
2666
2667 return ( w * rhoPL ) * ( w * dlnp ) + v * v * dadT ;
2668 }
2669
2670 inline real
2671 Material::rho_riva( const real normJ, const real T, const real normB, const real angleNxB,
2672 const real x, const real y, const real z, const real t ) const
2673 {
2674 real rhon = this->rho( T ) ;
2675
2676 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2677 {
2678 return rhon ;
2679 }
2680
2681 real ec = this->constant_property( MaterialProperty::ec ) ;
2682 real jc = this->jc_eval( T, normB, angleNxB ) * ( ( this->mDefectFunction )( x, y, z, t ) ) ;
2683 real n = this->n_eval( T, normB, angleNxB ) ;
2684
2685 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) )
2686 {
2687 return rhon ;
2688 }
2689
2690 real rhoPL ;
2691 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2692 {
2693 return rhon ;
2694 }
2695
2696 return rhoPL > rhon ? rhon / ( 1.0 + rhon / rhoPL )
2697 : rhoPL / ( 1.0 + rhoPL / rhon ) ;
2698 }
2699
2700 inline real
2701 Material::drho_riva_dJ( const real normJ, const real T, const real normB, const real angleNxB,
2702 const real x, const real y, const real z, const real t ) const
2703 {
2704 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2705 {
2706 return 0.0 ;
2707 }
2708
2709 real ec = this->constant_property( MaterialProperty::ec ) ;
2710 real jc = this->jc_eval( T, normB, angleNxB ) * ( ( this->mDefectFunction )( x, y, z, t ) ) ;
2711 real n = this->n_eval( T, normB, angleNxB ) ;
2712
2713 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) ) return 0.0 ;
2714 if ( n <= 1.0 ) return 0.0 ;
2715 if ( normJ < BELFEM_EPSILON ) return 0.0 ;
2716
2717 real rhoPL ;
2718 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2719 {
2720 return 0.0 ;
2721 }
2722
2723 real rhon = this->rho( T ) ;
2724 real w = rhon / ( rhoPL + rhon ) ;
2725
2726 // grouped so each factor stays bounded: w·ρPL ≤ ρn even when ρPL
2727 // is huge, and w ≤ 1
2728 return ( w * rhoPL ) * ( w * ( n - 1.0 ) / normJ ) ;
2729 }
2730
2731 inline real
2732 Material::drho_riva_dB( const real normJ, const real T, const real normB, const real angleNxB,
2733 const real x, const real y, const real z, const real t ) const
2734 {
2735 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2736 {
2737 return 0.0 ;
2738 }
2739
2740 real ec = this->constant_property( MaterialProperty::ec ) ;
2741 real tD = ( this->mDefectFunction )( x, y, z, t ) ;
2742 real jc = this->jc_eval( T, normB, angleNxB ) * tD ;
2743 real n = this->n_eval( T, normB, angleNxB ) ;
2744
2745 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) ) return 0.0 ;
2746
2747 // J < eps early-out only for n > 1, cf. the non-defect overload
2748 if ( n > 1.0 && normJ < BELFEM_EPSILON ) return 0.0 ;
2749
2750 // defect modulates jc AND its derivative by the same D
2751 real djcdB = this->djc_eval_dB( T, normB, angleNxB ) * tD ;
2752 real dndB = this->dn_eval_dB( T, normB, angleNxB ) ;
2753
2754 if ( djcdB == 0.0 && dndB == 0.0 ) return 0.0 ;
2755
2756 real rhoPL ;
2757 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2758 {
2759 return 0.0 ;
2760 }
2761
2762 real rhon = this->rho( T ) ;
2763 real w = rhon / ( rhoPL + rhon ) ;
2764
2765 // djc/jc grouped as a ratio ( bounded for every real jc source ) so
2766 // a tiny jc cannot overflow n/jc on its own; weights grouped so each
2767 // factor stays bounded ( w·ρPL ≤ ρn, w ≤ 1 )
2768 return ( w * rhoPL ) * ( w * ( -( n * ( djcdB / jc ) )
2769 + ( normJ < BELFEM_EPSILON ? 0.0
2770 : std::log( normJ ) - std::log( jc ) ) * dndB ) ) ;
2771 }
2772
2773 inline real
2774 Material::drho_riva_dT( const real normJ, const real T, const real normB, const real angleNxB,
2775 const real x, const real y, const real z, const real t ) const
2776 {
2777 real dadT = this->drhodT( T ) ;
2778
2779 if ( T > this->constant_property( MaterialProperty::T_crit ) )
2780 {
2781 return dadT ;
2782 }
2783
2784 real ec = this->constant_property( MaterialProperty::ec ) ;
2785 real tD = ( this->mDefectFunction )( x, y, z, t ) ;
2786 real jc = this->jc_eval( T, normB, angleNxB ) * tD ;
2787 real n = this->n_eval( T, normB, angleNxB ) ;
2788
2789 if ( ! ( jc > 0.0 && std::isfinite( jc ) ) ) return dadT ;
2790
2791 real rhoPL ;
2792 if ( ! this->riva_rho_pl( normJ, jc, n, ec, rhoPL ) )
2793 {
2794 return dadT ;
2795 }
2796
2797 // defect modulates jc AND its derivative by the same D
2798 real djcdT = this->djc_eval_dT( T, normB, angleNxB ) * tD ;
2799 real dndT = this->dn_eval_dT( T, normB, angleNxB ) ;
2800
2801 real rhon = this->rho( T ) ;
2802 real w = rhon / ( rhoPL + rhon ) ;
2803 real v = rhoPL / ( rhoPL + rhon ) ;
2804
2805 // djc/jc grouped as a ratio and the weights as (w·ρPL)·(w·dln) so
2806 // no factor can overflow on a tiny-but-finite jc
2807 real dlnp = ( normJ < BELFEM_EPSILON && n > 1.0 ) ? 0.0 :
2808 -( n * ( djcdT / jc ) )
2809 + ( normJ < BELFEM_EPSILON ? 0.0
2810 : std::log( normJ ) - std::log( jc ) ) * dndT ;
2811
2812 return ( w * rhoPL ) * ( w * dlnp ) + v * v * dadT ;
2813 }
2814
2815}
2816#endif //BELFEM_POWERLAWS_HPP
BELFEM_ERROR
#define BELFEM_ERROR(aCheck,...)
Definition assert.hpp:264
BELFEM_ASSERT
#define BELFEM_ASSERT(aCheck,...)
Definition assert.hpp:244
cl_Material.hpp
belfem::Material::is_constant
bool is_constant(const MaterialProperty aProperty) const
Check if a property is constant (temperature-independent).
Definition cl_Material.hpp:1781
belfem::Material::rho_powerlaw
real rho_powerlaw(const real normJ) const
Power law resistivity for HTS (constant jc and n).
Definition powerlaws.hpp:184
belfem::Material::rho
virtual real rho(const real T) const
Electrical resistivity (isotropic).
Definition cl_Material.hpp:2038
belfem::Material::jc_custom
virtual real jc_custom(const real T) const
Definition cl_Material.cpp:574
belfem::Material::n
real n(const real normB, const real angleNxB, const real T) const
Direct evaluation of the n-value as a function of , and .
Definition powerlaws.hpp:976
belfem::Material::n_eval_raw
real n_eval_raw(const real T, const real normB, const real angleNxB) const
Definition powerlaws.hpp:76
belfem::Material::drho_riva_dB
real drho_riva_dB(const real normJ, const real T, const real normB, const real angleNxB) const
dρ/d|B| of rho_riva: w²·dρ_PL/dB ( dρ_n/dB = 0 on this path )
Definition powerlaws.hpp:2585
belfem::Material::mJcFunction
const material::JcFunction * mJcFunction
Definition cl_Material.hpp:414
belfem::Material::drhodT
virtual real drhodT(const real T) const
Definition cl_Material.hpp:2049
belfem::Material::drho_piecewise_dT
real drho_piecewise_dT(const real normJ, const real T, const real normB, const real angleNxB) const
Derivative of piecewise resistivity with respect to T at fixed J, |B|, θ ( T-leg ).
Definition powerlaws.hpp:2171
belfem::Material::jc_eval
real jc_eval(const real T, const real normB, const real angleNxB) const
O1 "full-signature policy" helpers: the assembly path always passes the full ( T, normB,...
Definition powerlaws.hpp:66
belfem::Material::drho_piecewise_dB
real drho_piecewise_dB(const real normJ, const real T, const real normB, const real angleNxB) const
Derivative of piecewise resistivity with respect to |B| at fixed J, T, θ
Definition powerlaws.hpp:2020
belfem::Material::drho_piecewise_dJ
real drho_piecewise_dJ(const real normJ, const real T, const real normB, const real angleNxB) const
Derivative of piecewise resistivity with respect to current density magnitude.
Definition powerlaws.hpp:1788
belfem::Material::drho_riva_dJ
real drho_riva_dJ(const real normJ, const real T, const real normB, const real angleNxB) const
dρ/d|J| of rho_riva: w²·dρ_PL/dJ with w = ρ_n/(ρ_PL+ρ_n)
Definition powerlaws.hpp:2554
belfem::Material::dn_eval_dB
real dn_eval_dB(const real T, const real normB, const real angleNxB) const
d(n)/d|B|, same routing as djc_eval_dB
Definition powerlaws.hpp:110
belfem::Material::drho_powerlaw_dB
real drho_powerlaw_dB(const real normJ, const real T, const real normB, const real angleNxB) const
Derivative of power-law resistivity with respect to |B| at fixed J, T, θ ( jc = jc(T,...
Definition powerlaws.hpp:1964
belfem::Material::rho_piecewise
real rho_piecewise(const real normJ) const
Piecewise resistivity for HTS (constant jc and n).
Definition powerlaws.hpp:412
belfem::Material::djc_eval_dB
real djc_eval_dB(const real T, const real normB, const real angleNxB) const
Definition powerlaws.hpp:102
belfem::Material::riva_rho_pl
bool riva_rho_pl(const real normJ, const real jc, const real n, const real ec, real &rhoPL) const
Definition powerlaws.hpp:2489
belfem::Material::n_eval
real n_eval(const real T, const real normB, const real angleNxB) const
Definition powerlaws.hpp:86
belfem::Material::jc
real jc(const real normB, const real angleNxB, const real T) const
Direct evaluation of the critical current density .
Definition powerlaws.hpp:1012
belfem::Material::djc_eval_dT
real djc_eval_dT(const real T, const real normB, const real angleNxB) const
Definition powerlaws.hpp:126
belfem::Material::dn_eval_dT
real dn_eval_dT(const real T, const real normB, const real angleNxB) const
d(n)/dT, same routing as djc_eval_dT
Definition powerlaws.hpp:134
belfem::Material::drho_powerlaw_dT
real drho_powerlaw_dT(const real normJ, const real T, const real normB, const real angleNxB) const
Derivative of power-law resistivity with respect to T at fixed J, |B|, θ ( T-leg: jc(T),...
Definition powerlaws.hpp:2113
belfem::Material::constant_property
real constant_property(const MaterialProperty aProperty) const
Get constant property value.
Definition cl_Material.hpp:1787
belfem::Material::mNFunction
const material::JcFunction * mNFunction
Definition cl_Material.hpp:415
belfem::Material::drho_riva_dT
real drho_riva_dT(const real normJ, const real T, const real normB, const real angleNxB) const
dρ/dT of rho_riva: w²·dρ_PL/dT + (1−w)²·dρ_n/dT
Definition powerlaws.hpp:2628
belfem::Material::drho_powerlaw_dJ
real drho_powerlaw_dJ(const real normJ, const real T, const real normB, const real angleNxB) const
Derivative of power-law resistivity with respect to current density magnitude.
Definition powerlaws.hpp:1240
belfem::Material::rho_riva
real rho_riva(const real normJ, const real T, const real normB, const real angleNxB) const
Riva-law resistivity: the superconducting power-law channel in parallel with the normal-state channel...
Definition powerlaws.hpp:2524
belfem::Material::n_custom
virtual real n_custom(const real T) const
Definition cl_Material.cpp:581
belfem::material::JcFunction::deval_dB
virtual real deval_dB(const real normB, const real angle, const real T) const
Derivative of eval with respect to the field magnitude.
Definition cl_JcFunction.cpp:52
belfem::material::JcFunction::eval
virtual real eval(const real normB, const real angle) const
Evaluate with field and angle dependence.
Definition cl_JcFunction.cpp:26
belfem::material::JcFunction::depends_on
bool depends_on(const JcParameter aParameter) const
Check if function depends on a parameter.
Definition cl_JcFunction.hpp:281
belfem::material::JcFunction::deval_dT
virtual real deval_dT(const real normB, const real angle, const real T) const
Derivative of eval with respect to temperature.
Definition cl_JcFunction.cpp:72
belfem::material::JcParameter::angleNxB
@ angleNxB
Definition cl_JcFunction.hpp:82
belfem::material::JcParameter::normB
@ normB
Definition cl_JcFunction.hpp:81
belfem::material::JcParameter::T
@ T
Definition cl_JcFunction.hpp:83
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::cplx
std::complex< real > cplx
Definition typedefs.hpp:37
belfem::gTbulk
greal gTbulk
Definition globals.hpp:45
belfem::BELFEM_EPSILON
constexpr real BELFEM_EPSILON
Definition typedefs.hpp:90
belfem::MaterialProperty::ec
@ ec
Definition cl_Material.hpp:185
belfem::MaterialProperty::n
@ n
Definition cl_Material.hpp:187
belfem::MaterialProperty::jc
@ jc
Definition cl_Material.hpp:186
belfem::MaterialProperty::T_crit
@ T_crit
Definition cl_Material.hpp:167
belfem::real
double real
Definition typedefs.hpp:36