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